7.1.6 LTI DT – H[z] Diagrama de bloques con «1/z» – ejercicios

Referencia: Lathi 5.4 p519. Oppenheim 10.8 p783, Hsu Ejercicio 4.34 p201

Dada la similitud entre sistemas LTIC y LTID, las convenciones para diagramas de bloques  y las reglas de interconexión son idénticas a los sistemas contínuos. Los diagramas mostrados se realizan con el programa Xcos de SciLab.

Una función de transferencia general se expresa como:

H(z) = \frac{b_0 z^N +b_1 z^{N-1} + \text{ ... } + b_{N-1} z + b_N}{z^N + a_1 z^{N-1} +\text{ ... } + a_{N-1}z +a_N}

Ejemplo 1. H(z)

Referencia : Lathi 5.8a p522

H(z) =\frac{2}{z+5}

La función de transferencia es de primer orden (N=1), por lo que solo se usará un retraso para el diagrama. Los coeficientes de retroalimentación de retraso y adelanto son:

a1= 5 y b0 = 0, b1=2

La imagen presenta dos formas de realizar los diagramas

Ejemplo 2. H(z)

Referencia : Lathi 5.8b p522

H(z) =\frac{4z+28}{z+1}

la función de transferencia también es de primer orden (N=1), los coeficientes de retraso y adelanto son:

a1= 1 y b0 = 4, b1=28

La imagen presenta dos formas de realizar el diagrama de bloques

Ejemplo 3. H(z)

Referencia : Lathi 5.8c p522

H(z) =\frac{z}{z+7}

la función de transferencia también es de primer orden (N=1), los coeficientes de retraso y adelanto son:

a1= 7 y b0 = 1, b1=0

Ejemplo 4. H(z)

Referencia : Lathi 5.8d p522

H(z) =\frac{4z+28}{z^2+6z+5}

la función de transferencia también es de segundo orden (N=2), los coeficientes de retraso y adelanto son:

a1= 6, a2= 5 y b0 = 0, b1=4, b2=28

Ejemplo 5. Sistema LTI D

Referencia: Ejemplo Openheim 10.28 p784

Considere el sistema LTI causal descrito mediante:

y[n] - \frac{1}{4} y[n-1] = x[n]

con función del sistema:

H(z) =\frac{1}{1- \frac{1}{4} z^{-1} }

Aqui z-1 es la función del sistema con retraso unitario. El diagrama de bloques en la figura contiene un lazo de retroalimentación.

y[n] -\frac{1}{4} z^{-1} y[n] = x[n] y[n] = x[n] + \frac{1}{4} z^{-1} y[n]


Ejemplo 6. Función H(z)

Referencia: Ejemplo Openheim 10.29 p785

Considere un sistema LTI causal con función del sistema:

H(z) = \frac{ 1 - 2 z^{-1} }{ 1 - \frac{1}{4} z^{-1} }

La función se puede separar en bloque de denominador y numerador, semejante a los polinomios Q(E)  y P(E).

= \Bigg[\frac{1}{1 - \frac{1}{4}z^{-1}} \Bigg] \Bigg[ 1 - 2z^{-1} \Bigg]

Se observa en el diagrama (a) que z-1 se encuentra duplicado pues ambos bloques toman la misma señal para aplicarle un atraso, cada uno le da a una ganancia diferente. Por lo que un solo bloque z-1 puede realizar la misma operación.

En el mundo de los circuito digitales y componentes, la situación se interpreta como un componente duplicado, lo que tiene implicaciones de costos de implementación, recuerde el tema cuando construye un circuito digital usando puertas lógicas y registros.

7.1.5 Transformada z Inversa con Sympy-Python

La Transformada z  inversa, usa la Tabla de Pares f[n] y F[z] y los algoritmos para buscar una expresión semejante F.match() en conjunto con la tabla de propiedades y los algoritmos desarrollados en la sección anterior. El resultado se integra desde telg1001.py y se muestra en los siguientes ejemplos:

Ejemplo 1.  Transformada z inversa de una fracción X[1/z]

Referencia: Oppenheim Ejemplo 10.13 p761, Lathi Ejemplo 5.1 p490

Continuando con el ejercicio presentado para convertir x[n] a X[z], se tiene que:

X[z] = \frac{1}{1-a z^{-1}}\text{ , } |z|>|a| X[z] = \frac{z}{z-a}\text{ , } \Big|z|>|a|

siendo la expresión de ingreso al algoritmo:

F = z/(z - 2)

el resultado es:

f[n] : (2**n*Heaviside(n), Abs(z) > 2, True)
>>> 

Ejemplo 1a. Instrucciones Python

# transformada z propiedades con Sympy
# aplicar luego de buscar en tabla de pares
import sympy as sym
import telg1001 as fcnm

# INGRESO
n = sym.Symbol('n', real=True)
z = sym.Symbol('z')
u = sym.Heaviside(n)

# (a**n)*f[n] <--> F(z/a)
F = z/(z - 2)  #f = (2**n)*u

# PROCEDIMIENTO
fn = fcnm.inverse_z_transform(F,z,n)

# SALIDA
print('f[n] :',fn)

El algoritmo presentado puede ser tomado como punto de partida para los siguientes ejercicios. Realizado para explicación conceptual.

Ejemplo 2. Transformada z inversa usando expansión en fracciones parciales

Referencia: Lathi 5.3 a

Encuentre la transformada z inversa para:

F[z] = \frac{8z-19}{(z-2)(z-3)}

el ingreso para el algoritmo del ejercicio anterior es:

F = (8*z-19)/((z-2)*(z-3))

con el resultado siguiente:

f[n] : (9*2**n + 10*3**n)*Heaviside(n)/6 - 19*DiracDelta(n)/6

activando la variable sym.SYMPY_DEBUG=True, se observa el proceso de fracción parcial modificada, la búsqueda en la tabla, y aplicación de una propiedad de la transformada z.

inverse_z_transform(F,z,n): (8*z - 19)/((z - 3)*(z - 2))

apart_z(F) : (8*z - 19)/((z - 3)*(z - 2))
_simplify_z: (8*z - 19)/((z - 3)*(z - 2))
  P_leadcoef: 8  ; P_zc: 1  ; P: z - 19/8
  _simplify_z k_sign: 1  ; k: 8  ; P_zc: 1  ; Q_zc: 1
    F[z]: (z - 19/8)/((z - 3)*(z - 2))
_simplify_z: (z - 19/8)/(z*(z - 3)*(z - 2))
  P_leadcoef: 8  ; P_zc: 1  ; P: z - 19/8
  Q_leadcoef: 8  ; Q_zc: z  ; Q: (z - 3)*(z - 2)
  _simplify_z k_sign: 1  ; k: 1  ; P_zc: 1  ; Q_zc: z
    F[z]: (z - 19/8)/((z - 3)*(z - 2))
  P: (8*z - 19)/8 
  Q: (z - 3)*(z - 2) 
  P_degree: 1  ; P_zeros: {19/8: 1} 
  Q_degree: 3  ; Q_poles: {3: 1, 2: 1, 0: 1}
  apart(F/z): 3/(16*(z - 2)) + 5/(24*(z - 3)) - 19/(48*z)
  ma_Qk: {a_: 48, b_: 0, c_: 1}  ; ma_Qk2: {b_: 48, a_: 0, c_: 0}
_simplify_z: -19/(48*z)
  P_leadcoef: -19  ; P_zc: 1  ; P: 1
  Q_leadcoef: 48  ; Q_zc: z  ; Q: 1
  _simplify_z k_sign: -1  ; k: 19/48  ; P_zc: 1  ; Q_zc: z
    F[z]: 1
  ma_Qk: {a_: 16, b_: -32, c_: 1}  ; ma_Qk2: {b_: 16, c_: -32, a_: 0}
_simplify_z: 3/(16*(z - 2))
  P_leadcoef: 3  ; P_zc: 1  ; P: 1
  Q_leadcoef: 16  ; Q_zc: 1  ; Q: z - 2
  _simplify_z k_sign: 1  ; k: 3/16  ; P_zc: 1  ; Q_zc: 1
    F[z]: 1/(z - 2)
  ma_Qk: {a_: 24, b_: -72, c_: 1}  ; ma_Qk2: {b_: 24, c_: -72, a_: 0}
_simplify_z: 5/(24*(z - 3))
  P_leadcoef: 5  ; P_zc: 1  ; P: 1
  Q_leadcoef: 24  ; Q_zc: 1  ; Q: z - 3
  _simplify_z k_sign: 1  ; k: 5/24  ; P_zc: 1  ; Q_zc: 1
    F[z]: 1/(z - 3)

 apart z*(Fz/z)  : 3*z/(2*(z - 2)) + 5*z/(3*(z - 3)) - 19/6 

 apart_z(F) as partial fractions: (-19/6, 3*z/(2*(z - 2)), 5*z/(3*(z - 3)))
  term_k[z]: -19/6
[0] _z_pairs_prop_inverse ----
_simplify_z: -19/6
  P_leadcoef: -19  ; P_zc: 1  ; P: 1
  Q_leadcoef: 6  ; Q_zc: 1  ; Q: 1
  _simplify_z k_sign: -1  ; k: 19/6  ; P_zc: 1  ; Q_zc: 1
    F[z]: 1
 _z_pairs_table match: 
  F: 1
  z_pair F[z]: 1  ; ma_z: {}
  z_pair f[n]: DiracDelta(n)
  try,check  : True -> True
[0]   pair f[n]: (-19*DiracDelta(n)/6, True, True)
  term_k[n]: (-19*DiracDelta(n)/6, True, True)
  term_k[z]: 3*z/(2*(z - 2))
[0] _z_pairs_prop_inverse ----
_simplify_z: 3*z/(2*(z - 2))
  P_leadcoef: 3  ; P_zc: 1  ; P: z
  Q_leadcoef: 2  ; Q_zc: 1  ; Q: z - 2
  _simplify_z k_sign: 1  ; k: 3/2  ; P_zc: 1  ; Q_zc: 1
    F[z]: z/(z - 2)

[1] z_properties_inverse_z_transform
  P: z 
  Q: z - 2 
  P_degree: 1  ; P_zeros: {0: 1} 
  Q_degree: 1  ; Q_poles: {2: 1}
  ma_P1 (a*z+ b):    {a_: 1, b_: 0}
  ma_Q1 (a*z-b)**c : {a_: 1, b_: 2, c_: 1}

 _z_property multiply (a**n)*f[n] <--> F(z/a) 
 k_sign: 1 ; k: 3/2  ; k_b: 2
 F[z]: z/(z - 1)
[2] _z_pairs_prop_inverse ----
_simplify_z: z/(z - 1)
    F[z]: z/(z - 1)
 _z_pairs_table match: 
  F: z/(z - 1)
  z_pair F[z]: z/(z - 1)  ; ma_z: {}
  z_pair f[n]: Heaviside(n)
  try,check  : True -> True
[2]   pair f[n]: (Heaviside(n), Abs(z) > 1, True)

  _z_property multiply (a**n)*f[n]:
  (3*2**n*Heaviside(n)/2, Abs(z) > 2, True)
[1]  _z_properties_inverse f[n]: (3*2**n*Heaviside(n)/2, Abs(z) > 2, True)
[0]   prop f[n]: (3*2**n*Heaviside(n)/2, Abs(z) > 2, True)
  term_k[n]: (3*2**n*Heaviside(n)/2, Abs(z) > 2, True)
  term_k[z]: 5*z/(3*(z - 3))
[0] _z_pairs_prop_inverse ----
_simplify_z: 5*z/(3*(z - 3))
  P_leadcoef: 5  ; P_zc: 1  ; P: z
  Q_leadcoef: 3  ; Q_zc: 1  ; Q: z - 3
  _simplify_z k_sign: 1  ; k: 5/3  ; P_zc: 1  ; Q_zc: 1
    F[z]: z/(z - 3)

[1] z_properties_inverse_z_transform
  P: z 
  Q: z - 3 
  P_degree: 1  ; P_zeros: {0: 1} 
  Q_degree: 1  ; Q_poles: {3: 1}
  ma_P1 (a*z+ b):    {a_: 1, b_: 0}
  ma_Q1 (a*z-b)**c : {a_: 1, b_: 3, c_: 1}

 _z_property multiply (a**n)*f[n] <--> F(z/a) 
 k_sign: 1 ; k: 5/3  ; k_b: 3
 F[z]: z/(z - 1)
[2] _z_pairs_prop_inverse ----
_simplify_z: z/(z - 1)
    F[z]: z/(z - 1)
 _z_pairs_table match: 
  F: z/(z - 1)
  z_pair F[z]: z/(z - 1)  ; ma_z: {}
  z_pair f[n]: Heaviside(n)
  try,check  : True -> True
[2]   pair f[n]: (Heaviside(n), Abs(z) > 1, True)

  _z_property multiply (a**n)*f[n]:
  (5*3**n*Heaviside(n)/3, Abs(z) > 3, True)
[1]  _z_properties_inverse f[n]: (5*3**n*Heaviside(n)/3, Abs(z) > 3, True)
[0]   prop f[n]: (5*3**n*Heaviside(n)/3, Abs(z) > 3, True)
  term_k[n]: (5*3**n*Heaviside(n)/3, Abs(z) > 3, True)

 f[n]: 3*2**n*Heaviside(n)/2 + 5*3**n*Heaviside(n)/3 - 19*DiracDelta(n)/6
 f[n]: (9*2**n + 10*3**n)*Heaviside(n)/6 - 19*DiracDelta(n)/6
f[n] : (9*2**n + 10*3**n)*Heaviside(n)/6 - 19*DiracDelta(n)/6
>>> 

7.1.4 Transformada z Inversa – Conceptos

Referencia: Lathi 5.1-1 p495, Oppenheim 10.3 p757

Muchas de las transformadas X[z] de interés práctico son fracciones de polinomios en z, que pueden ser expresadas como una suma de fracciones parciales cuya inversa puede ser revisada en una tabla de transformadas z.

El método de fracciones parciales funciona dado que para cada x[n] transformable con n≥0 le corresponde una X[z] única definida para |z|>a y viceversa. Para lo mostrado, se supone que se usa la fracción parcial mas simple.

Ejemplo 1a.  Transformada z inversa de una fracción X[1/z]

Referencia: Oppenheim Ejemplo 10.13 p761, Lathi Ejemplo 5.1 p490

Continuando con el ejercicio presentado para convertir x[n] a X[z], se tiene que:

X[z] = \frac{1}{1-a z^{-1}}\text{ , } |z|>|a|

la expresión se puede expandir en series se potencias y obtener los primeros términos de la secuencia de x[n]. Los términos son útiles si se se requieren los m primeros términos de la secuencia, es decir n se encuentra en [0,m-1]

X[z]: 
  z   
------
-a + z

 F[z]:
         2    3    4    5    6                 
    a   a    a    a    a    a     /1          \
1 + - + -- + -- + -- + -- + -- + O|--; z -> oo|
    z    2    3    4    5    6    | 7         |
        z    z    z    z    z     \z          /
Términos x[n] entre[0,7]
[1, a, a**2, a**3, a**4, a**5, a**6]

 x[n] con a= 1/2
[1.0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625] 

Al final de la serie F[z] se muestra un término con el orden del error.

De los términos mostrados, se observa que X[n] = an

Usando por ejemplo a=1/2, se obtiene como resultado de x[n],

Transformada Z Ej01 X[z] x[n]

>>> xi
[1.0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625]

Ejemplo 1a. Instrucciones Python

# transformada z inversa de X[z]
import sympy as sym

# INGRESO
z = sym.symbols('z',)
n = sym.symbols('n', positive=True)
a = sym.symbols('a')
u = sym.Heaviside(n)

Xz = 1/(1-a*z**(-1))

# valor a como racional en dominio 'ZZ' enteros
a_k = sym.Rational(1/2).limit_denominator(100)

m = 7   # Términos a graficar

# PROCEDIMIENTO
# Series de potencia z**(-1)
Fz = sym.series(Xz,z**(-1), n=m)

# Terminos de X[n]
xn = []
termino = Fz.args
for i in range(0,m,1):
    xn.append(termino[i]*(z**i))
 
# SALIDA
print('X[z]: ')
sym.pprint(Xz)
print('\n F[z]:')
sym.pprint(Fz)
print('\n Términos x[n] entre [0,'+str(m)+']')
print(xn)

# GRAFICA valores ---------------
import numpy as np
import matplotlib.pyplot as plt
# Terminos de X[n]
xi = [] ; ki=[]
for i in range(0,m,1):
    valor = xn[i].subs({z:1,a:a_k})
    xi.append(float(valor))
    ki.append(i)

print('\n x[n] con a=',a_k)
print(xi)

# grafica entrada x[n]
plt.axvline(0,color='grey')
plt.stem(ki,xi,label='x[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.title(r'x[n]= $'+str(sym.latex(xn))+'$')
plt.grid()
plt.show()

Ejemplo 1b.  Transformada z inversa de una fracción dada con X[z]

Referencia: Oppenheim Ejemplo 10.13 p761, Lathi Ejemplo 5.1 p490

Si la expresión se usa de la forma simplificada, que es la más común para tratar las expresiones con Sympy, el algoritmo anterior para las series no realiza la misma respuesta.

X[z] = \frac{z}{z-a}\text{ , } \Big|z|>|a| X[z] = \frac{1}{1-a z^{-1}}\text{ , } \Big|z|>|a|

Para este caso se multiplica el numerador y denominador por 1/z, obteniendo la expresión del ejemplo anterior. Para generar la serie se realiza el, cambio de variable Z=1/z para generar la serie sobre Z. Luego se restablece en la expresión antes de mostrar los términos, y se obtiene el mismo resultado.

De los términos mostrados, se observa que X[n] = an

Ejemplo 1b. Instrucciones en Python

# transformada z inversa de X[z]
# supone que la expresión son fracciones parciales
import sympy as sym

# INGRESO
z = sym.symbols('z')
n = sym.symbols('n', integer=True, positive=True)
a = sym.symbols('a')

Xz = z/(z-a)

# valor a como racional en dominio 'ZZ' enteros
a_k = sym.Rational(1/2).limit_denominator(100)

m = 7   # Términos a graficar

# PROCEDIMIENTO
# expresión Xz con z para aplicar serie
# separa numerador y denominador
[Pz,Qz] = Xz.as_numer_denom()
Pz = (Pz*(1/z)).expand(1/z)
Qz = (Qz*(1/z)).expand(1/z)

# cambia Z por 1/z
Z = sym.symbols('Z')
PZ = Pz.subs(1/z,Z)
QZ = Qz.subs(1/z,Z)
XZ = PZ/QZ

# Series de potencia de Z
FZ = sym.series(XZ,Z, n=m)
Fz = FZ.subs(Z,1/z) # restituye 1/z

# Terminos de X[n]
xn = []
termino = Fz.args
for i in range(0,m,1):
    xn.append(termino[i]*(z**i))

# SALIDA
print('X[z]: ')
sym.pprint(Xz)
print('\n F[z]:')
sym.pprint(Fz)
print('Términos x[n] entre[0,'+str(m)+']')
print(xn)

# GRAFICA valores ---------------
import numpy as np
import matplotlib.pyplot as plt
# Terminos de X[n]
xi = [] ; ki=[]
for i in range(0,m,1):
    valor = xn[i].subs({z:1,a:a_k})
    xi.append(float(valor))
    ki.append(i)

print('\n x[n] con a=',a_k)
print(xi)

# grafica entrada x[n]
plt.axvline(0,color='grey')
plt.stem(ki,xi,label='x[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.title(r'x[n]= $'+str(sym.latex(xn))+'$')
plt.grid()

plt.show()

El algoritmo presentado puede ser tomado como punto de partida para los siguientes ejercicios. Realizado para explicación conceptual.

7.1 Transformada z – Sumatoria con Sympy-Python

Referencia: Lathi 5.1 p488, Oppenheim 10.1 p741

Se define como X[z] a la transformada z de una señal x[n] como

X[z] = \sum_{n=-\infty}^{\infty} x[n] z^{-n}

Por razones semejantes a las descritas en la unidad 4 para transformadas de Laplace, es conveniente considerar la transformada z unilateral. Dado que muchos de los sistemas y señales son causales, en la práctica se considera que las señales inician en n=0 (señal causal), la definición de la transformada z unilateral es la misma que la bilateral. excepto por los limites de la sumatoria que son [0,∞]

X[z] = \sum_{n=0}^{\infty} x[n] z^{-n}

Ejemplo 1. Transformada z de exponencial causal

Referencia: Lathi ejemplo 5.1 p490, Oppenheim ejemplo 10.1 p743

Encuentre la Transformada z y la correspondiente ROC para la señal x[n]

x[n] = a ^{n} \mu[n]

Transformada Z Ej01x[n]

1.1 Desarrollo analítico

por definición

X[z] = \sum_{n=0}^{\infty} a ^{n} \mu[n] z^{-n}

dado que μ[n]=1 para todo n≥0,

X[z] = \sum_{n=0}^{\infty} \Big( a^n \Big) \Big(\frac{1}{z} \Big)^{n} = \sum_{n=0}^{\infty} \Big(\frac{a}{z} \Big)^{n} =1 + \Big(\frac{a}{z} \Big)^{1}+\Big(\frac{a}{z} \Big)^{2}+\Big(\frac{a}{z} \Big)^{3}+\text{...}

revisando la progresión geométrica,

1+x+x^2+x^3+\text{...} = \frac{1}{1+x} \text{ , }|x|<1

y aplicando en la expresión del problema, se tiene:

X[z] = \frac{1}{1-\frac{a}{z}}\text{ , } \Big|\frac{a}{z}\Big| <1 X[z] = \frac{z}{z-a}\text{ , } |z|>|a|

observe que X[z] existe solo si |z|>|a|. Para |z|<|a| la sumatoria no converge y tiende al infinito. Por lo que, la Región de Convergencia ROC de X[z] es la región sombreada de un círculo de radio |a| centrado en el origen en el plano z.

Transformada Z Ej01 ROC
Luego se puede mostrar que la transformada z de la señal -an u[-(n+1)] también es z/(z-a). Sin embargo la región de convergencia en éste caso es |z|<|a|. Se observa que la transformada z inversa no es única, sin embargo,al restringir la transformada inversa a causal, entonces la transformada inversa es única como la mostrada.

Transformada Z Ej01 Xz

1.2 Transformada z con Sympy-Python

En Sympy de Python existe la función sym.summation(), que se usa para generar la expresión del resultado para la transformada z.

f[n]: 
 n
a 

 F[z] desde sumatoria:
    1     |a|     
(-------, |-| < 1)
   a      |z|     
 - - + 1          
   z              

 F[z] simplificada
  z   
------
-a + z

 ROC: 
|a|    
|-| < 1 
|z| 

{Q_polos:veces}: {a: 1} 
{P_ceros:veces}: {0: 1} 
>>> 

Instrucciones en Python

# transformada z de x[n]u[n]
import sympy as sym

# INGRESO
z = sym.symbols('z')
n = sym.symbols('n', positive=True)
a = sym.symbols('a')
u = sym.Heaviside(n)

fn = (a**n)*u

# valor a como racional en dominio 'ZZ' enteros
a_k = sym.Rational(1/2).limit_denominator(100)
m   = 7        # Términos a graficar
muestras = 101 # dominio z

# PROCEDIMIENTO
fnz = fn*(z**(-n)) # f(n,z) para sumatoria
# sumatoria transformada z
Fz_sum = sym.summation(fnz,(n,0,sym.oo))
Fz_eq  = Fz_sum.args[0]  # primera ecuacion e intervalo
Fz = Fz_eq[0].simplify() # solo expresion

ROC = Fz_eq[1]  # condicion ROC

# polos y ceros de Fz
[P,Q] = Fz.as_numer_denom()
P = sym.poly(P,z)
Q = sym.poly(Q,z)
P_ceros = sym.roots(P)
Q_polos = sym.roots(Q)

# SALIDA
print('f[n]: ')
sym.pprint(fn)
print('\n F[z] desde sumatoria:')
sym.pprint(Fz_eq)
print('\n F[z] simplificada')
sym.pprint(Fz)
print('\n ROC: ')
sym.pprint(ROC)
print('\n {Q_polos:veces}:',Q_polos)
print(' {P_ceros:veces}:',P_ceros)

La parte gráfica se desarrolla reutilizando algunas funciones de los algoritmos telg1001.py

Ser requiere dar valor a la constante ‘a‘ y se sustituye en las funciones con a_k=1/2. Para actualizar los polos, se usa la función buscapolos() realizada en la unidad 4.3.2, añadiendo en la función la sustitución en F[z] de variable independiente a F(s) para que sea transparente el cálculo. Lo mismo se aplica a la grafica de la transformada F[z] y los resultados son los mostrados.

# GRAFICA  -----------
import numpy as np
import matplotlib.pyplot as plt
import telg1001 as fcnm

def graficar_Fz_polos(Fz,Q_polos={},P_ceros={},
                z_a=1,z_b=0,muestras=101,f_nombre='F',
                solopolos=False):
    ''' polos y ceros plano z imaginario
    '''
    fig_zROC, graf_ROC = plt.subplots()
    # limite con radio 1
    radio1 = plt.Circle((0,0),1,color='lightsalmon',
                        fill=True)
    radio2 = plt.Circle((0,0),1,linestyle='dashed',
                        color='red',fill=False)
    graf_ROC.add_patch(radio1)
    for unpolo in Q_polos.keys():
        [r_real,r_imag] = unpolo.as_real_imag()
        unpolo_radio = np.abs(unpolo)
        unpolo_ROC = plt.Circle((0,0),unpolo_radio,
                          color='lightgreen',fill=True)
        graf_ROC.add_patch(unpolo_ROC)
    graf_ROC.add_patch(radio2) # borde r=1
    graf_ROC.axis('equal')
    # marcas de r=1 y polos
    for unpolo in Q_polos.keys():
        x_polo = sym.re(unpolo)
        y_polo = sym.im(unpolo)
        etiqueta = 'polo: '+str(unpolo)
        graf_ROC.scatter(x_polo,y_polo,marker='x',
                        color='red',label = etiqueta)
        etiqueta = "("+str(float(x_polo)) + ','
        etiqueta = etiqueta + str(float(y_polo))+")"
        plt.annotate(etiqueta,(x_polo,y_polo), rotation=45)
    # marcas de ceros
    for uncero in P_ceros.keys():
        x_cero = sym.re(uncero)
        y_cero = sym.im(uncero)
        etiqueta = 'cero: '+str(uncero)
        graf_ROC.scatter(x_cero,y_cero,marker='o',
                        color='blue',label = etiqueta)
        etiqueta = "("+str(float(x_cero)) + ','
        etiqueta = etiqueta + str(float(y_cero))+")"
        plt.annotate(etiqueta,(x_cero,y_cero), rotation=45)
    # limita radio 1
    graf_ROC.plot(1,0,'o',color='red',
                 label ='radio:'+str(1))
    graf_ROC.axhline(0,color='grey')
    graf_ROC.axvline(0,color='grey')
    graf_ROC.grid()
    graf_ROC.legend()
    graf_ROC.set_xlabel('Re[z]')
    graf_ROC.set_ylabel('Imag[z]')
    untitulo = r'ROC '+f_nombre+'[z]=$'
    untitulo = untitulo+str(sym.latex(Fz))+'$'
    graf_ROC.set_title(untitulo)
    return(fig_zROC)

# a tiene valor a_k
fn = fn.subs(a,a_k) 
Fz = Fz.subs(a,a_k)
# polos y ceros de Fz
polosceros = fcnm.busca_polosceros(Fz)
Q_polos = polosceros['Q_polos']
P_ceros = polosceros['P_ceros']
# estima intervalo para z
z_a = list(Q_polos.keys())
z_a.append(1) # comparar con radio 1
z_a = 2*int(max(z_a))

fig_ROC = graficar_Fz_polos(Fz,Q_polos,P_ceros,
                      muestras=101,f_nombre='X')

fig_Fz = fcnm.graficar_Fs(Fz,Q_polos,P_ceros,
                     -z_a,z_a,muestras=101,
                     f_nombre='X')

# para graficar f[n]
f_n = sym.lambdify(n,fn)
ki  = np.arange(0,m,1)
fi  = f_n(ki)

# entrada x[n]
fig_fn, grafxn = plt.subplots()
plt.axvline(0,color='grey')
plt.stem(ki,fi)
plt.grid()
plt.xlabel('n')
plt.ylabel('x[n]')
etiqueta = r'x[n]= $'+str(sym.latex(fn))+'$'
plt.title(etiqueta)

plt.show()

Referencia: The z-transform.dynamics-and-control. https://dynamics-and-control.readthedocs.io/en/latest/1_Dynamics/9_Sampled_systems/The%20z%20transform.html

Convolución de Sumas – Tabla

Referencia: Lathi Tabla 3.1 p285

\gamma_1 \neq \gamma_2
No x1[n] x2[n] x1[n]⊗x2[n] = x2[n]⊗x1[n]
1 δ[n-k] x[n] x[n-k]
2 \gamma^{n} \mu[n] μ[n] \frac{1-\gamma^{n+1}}{1-\gamma} \mu[n]
3 μ[n] μ[n] (n+1) μ[n]
4 \gamma_1^{n} \mu[n] \gamma_2^{n} \mu[n] \frac{\gamma_1^{n+1} - \gamma_2^{n+1}}{\gamma_1 - \gamma_2} \mu[n]
5 \gamma_1^{n} \mu[n] \gamma_2^{n} \mu[-(n+1) ] \frac{\gamma_1}{\gamma_2 - \gamma_1} \gamma_1^{n} \mu[n] +
+ \frac{\gamma_2}{\gamma_2 - \gamma_1} \gamma_2^{n} \mu[-(n+1)]
|\gamma_2| > |\gamma_1|
6 n\gamma_1^{n} \mu[n] \gamma_2^{n} \mu[n] \frac{\gamma_1 \gamma_2}{(\gamma_1 - \gamma_2)^2} \Big[ \gamma_2^{n} - \gamma_1^{n} + \frac{\gamma_1 - \gamma_2}{\gamma_2}n \gamma_1^n \Big] \mu [n]
\gamma_1\neq \gamma_2
7 n μ[n] n μ[n] \frac{1}{6} n (n-1) (n+1) \mu [n]
8 \gamma^{n} \mu[n] \gamma^{n} \mu[n] (n+1) \gamma^{n} \mu[n]
9 \gamma^{n} \mu[n] n \mu[n] \Big[ \frac{\gamma(\gamma^{n}-1)+n(1-\gamma)}{(1-\gamma)^2} \Big] \mu[n]
10 |\gamma_1|^{n} \cos (\beta n + \theta) \mu [k] |\gamma_2|^{n}\mu [n] \frac{1}{R} \Big[ |\gamma_1|^{n+1} \cos [\beta (n+1) +\theta -\phi] -\gamma_2 ^{n+1} \cos (\theta - \phi) \Big] \mu[n]
R=\Big[|\gamma_1|^2 + \gamma_2^2 -2|\gamma_1|\gamma_2 \cos(\beta) \Big]^{\frac{1}{2}}

\phi = \tan ^{-1} \Big[ \frac{|\gamma_1| \sin(\beta)}{|\gamma_1| \cos (\beta) -\gamma_2} \Big]

11 \mu [n] n\mu [n] \frac{n(n+1)}{2}\mu [n]

6.6 LTI DT – Respuesta Total

Referencia:  Lathi 3.8-3 p297

La respuesta total de un sistema LTID se puede expresar como la suma respuesta a entrada cero y respuesta a estado cero.

\text{respuesta total} = \sum_{j=1}^{N} c_j \gamma_j^{n} + x[n] \circledast h[n]

respuesta total = componente entrada cero + componente estado cero

Es el concepto aplicado al modelo de sistemas contínuos, cambiando a x[n], h[n] y y[n]


Para el sistema LTID desarrollado en las secciones anteriores y descrito por la ecuación,

y[n+2] - 0.6 y[n+1] - 0.16 y[n] = 5x[n+2]

dadas las condiciones iniciales y[-1]=0, y[-2]=25/4, ante una entrada
x[n]=4-nμ[n], se han determinado los dos componentes, ambos para n≥0:

\text{respuesta total} = y [n] = \frac{1}{5} (-0.2)^n + \frac{4}{5} (0.8)^n + + \Big[-1.26 (0.25)^{n} + 0.444 (-0.2)^{n} +5.81(0.8)^{n}\Big]

al simplificar las expresiones se tiene la respuesta expresada como

y [n] = \Big(\frac{1}{5}+ 0.444 \Big) (-0.2)^n + \Big(\frac{4}{5} +5.81 \Big) (0.8)^n -1.26 (0.25)^{n} y [n] = 0.644 (-0.2)^n + 6.61 (0.8)^n -1.26 (0.25)^{n} ; n \geq 0

Para el resultado se ha asumido que es un sistema invariante en el tiempo LTI, entonces la respuesta a δ[n-m] se puede expresar como h[n-m].


Solución clásica a ecuaciones lineales de diferencias

considera la respuesta total como la suma de una respuesta natural y una respuesta a componente forzados o de entrada.

\text{respuesta total} = y_c[n] +y_0[n]

se expresa como,

Q[E](y_c[n] + y_{\phi}[n]) = P[E] x[n]

dado que yc[n] es el resultado de los modos característicos,

Q[E]y_c[n] = 0 Q[E] y_{\phi}[n] = P[E] x[n]

La respuesta natural es una combinación lineal de los modos característicos. Las constantes arbitraras se determinan de las condiciones auxiliares dadas como y[0], y[1], … y[n-1], o por condiciones iniciales y[-1], y[-2],…, y[-N].

Las respuesta forzadas yΦ[n] satisface la ecuación anterior y por definición contiene solamente los términos que «nomodos»

6.5 LTI DT – Convolución de sumas – Python

Referencia: Lathi 3.8 p282, Oppenheim 2.1.2 p77 pdf108, Hsu 2.6.C p62

Una señal discreta de entrada x[n] se representa como una superposición de versiones escaladas de un conjunto de impulsos unitarios desplazados δ[n-k], cada uno con valor diferente de cero en un solo punto en el tiempo, especificado por el valor de k.

La respuesta de un sistema lineal y[n] a x[n] es la superposición de las respuestas escaladas del sistema a cada uno de estos impulsos desplazados.

Si usamos hk[n] como la respuesta del sistema lineal al impulso unitario desplazado por δ[n-k]. La respuesta y[n] del sistema lineal a la entrada x[n] en la ecuación será la combinación lineal ponderada de las respuesta básicas.

y[n] = \sum_{k=-\infty}^{+\infty} x[k]h_{k}[n]

Si se conoce la respuesta de un sistema lineal al conjunto de impulsos unitarios desplazados, podemos construir una respuesta a una entrada arbitraria.

Si el sistema lineal también es invariante en el tiempo, entonces estas respuestas a impulsos unitarios desplazados en el tiempo son todas las versiones desplazadas en el tiempo unas de otras.

h[n] es la salida del sistema LTI cuando δ[n] es la entrada. Entonces para un sistema LTI la ecuación se vuelve.

y[n] = \sum_{k=-\infty}^{+\infty} x[k]h[n-k]

El resultado se conoce como la «convolución de suma» o «suma de superposición» y la operación miembro derecho de la ecuación se llama convolución de las secuencias x[n] y h[n] que se representa de manera simbólica como:

y[n] = x[n]*h[n]

Ejemplo

Referencia: Openheim Ejemplo 2.4 p85 pdf116

Considere una entrada x[n] y una respuesta al impulso unitario h[n] dada por:

x[n] = (u[n]-u[n-5])
h[n] = αn (u[n]-u[n-7])

con α>1.  Para el ejercicio, α=1.5.

Nota: Para el algoritmo, si α es un entero, por ejemplo 2, usar α=2.0, para que la operación potencia se realice con números reales, como entero se puede saturar y dar error.


Al aplicar la suma de convolución obtiene:

Para aplicar el algoritmo se requiere definir u[n], por ser parte de las funciones x[n] y h[n]. Dado que la operación requiere valores fuera del rango muestreado para n, la sección suma convolución utiliza las funciones en lugar de los vectores xi, hi.

La función está definida en un intervalo simétrico, por lo que el rango de trabajo [a,b] se mantiene de la forma [-b,b] en las instrucciones.

# Respuesta a impulsos - forma discreta
# Ejemplo Oppenheim Ejemplo 2.3 p83/pdf111
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Rango [a,b], simétrico a 0
b = 15 ; a = -b

alfa = 1.5
u = lambda n: np.piecewise(n,n>=0,[1,0])

x = lambda n: u(n)-u(n-5)
h = lambda n: (alfa**n)*(u(n)-u(n-7))


# PROCEDIMIENTO
ni = np.arange(a,b+1,1)
xi = x(ni)
hi = h(ni)

# Suma de Convolucion x[n]*h[n]
muestras = len(xi)
yi = np.zeros(muestras, dtype=float)
for i in range(0,muestras):
    suma = 0
    for k in range(0,muestras):
        suma = suma + x(ni[k])*h(ni[i]-ni[k])
    yi[i] = suma

# yi = np.convolve(xi,hi,'same')

# SALIDA - GRAFICA
plt.figure(1)
plt.suptitle('Suma de Convolución x[n]*h[n]')

plt.subplot(311)
plt.stem(ni,xi, linefmt='b--',
         markerfmt='bo',basefmt='k-')
plt.ylabel('x[n]')

plt.subplot(312)
plt.stem(ni,hi, linefmt='b--',
         markerfmt='ro',basefmt='k-')
plt.ylabel('h[n]')

plt.subplot(313)
plt.stem(ni,yi, linefmt='g-.',
         markerfmt='mo', basefmt='k-')

plt.ylabel('x[n]*h[n]')
plt.xlabel('n')

plt.show()

La suma convolución se encuentra también disponible con Numpy en np.convolve(), la  sección de suma convolución se puede reemplazar y obtener los mismos resultados. Considere que para éste caso se usan los vectores xi y hi.

yi = np.convolve(xi,hi,'same')

el algoritmo se puede aplicar a otros ejercicios para comprobar los resultados.

Continúe con el ejercicio 2.5 del libro.

6.4.1 LTI DT – Respuesta estado cero. Ejercicio con Python

Referencia: Lathi Ejemplo 3.21 p286

Determine la respuesta a estado cero de un sistema LTID descrito por la ecuación de diferencias mostrada,

y[n+2] - 0.6 y[n+1] - 0.16 y[n] = 5x[n+2]

Si la entrada es x[n] = 4-n u[n]


Para realizar el diagrama, se desplaza la ecuación en dos unidades

y[n] - 0.6 y[n-1] - 0.16 y[n-2] = 5x[n] y[n] = 0.6 y[n-1] + 0.16 y[n-2] + 5x[n]

las condiciones iniciales es estado cero son todas cero.


La entrada puede expresada como:

x[n] = 4^{-n} \mu [n] = \Big( \frac{1}{4} \Big)^n \mu [n] x[n] = 0.25^{n} \mu [n]

la respuesta al impulso del sistema se encontró en el ejemplo de la sección anterior:

h[n] = [ (-0.2)^n +4 (0.8)^n ] \mu [n]

Por lo que:

y[n] = x[n] \circledast h[n] = 0.25^{n} \mu [n] \circledast \Big[ (-0.2)^n +4 (0.8)^n \Big] \mu [n] = 0.25^{n} \mu [n] \circledast (-0.2)^n \mu [n] + 0.25^{n} \mu [n] \circledast 4 (0.8)^n \mu [n]

usando la tabla de convolución de sumas para tener:

y[n] = \Bigg[ \frac{0.25^{n+1}-(-0.2)^{n+1}}{0.25-(-0.2)} + 4 \frac{0.25^{n+1} - 0.8^{n+1}}{0.25-0.8}\Bigg] \mu [n] = \Bigg[2.22 \Big[ 0.25^{n+1} - (-0.2)^{n+1}\Big] + - 7.27 \Big[ 0.25^{n+1} - (0.8)^{n+1}\Big]\Bigg] \mu [n] = \Big[ -5.05 (0.25)^{n+1} - 2.22(-0.2)^{n+1} + 7.27 (0.8)^{n+1} \Big] \mu [n]

recordando que ϒn+1 = ϒ(ϒ)n se puede expresar,

y[n] = \Big[-1.26 (0.25)^{n} + 0.444 (-0.2)^{n} +5.81(0.8)^{n}\Big] \mu [n]

Desarrollo numérico con Python

La respuesta a estado cero se encuentra en forma numérica con la expresión:

y[n] = x[n] \circledast h[n]

que en Numpy es np.convolve() y facilita el cálculo de la convolucion de sumas. El tema de convolución con Python se desarrolla en la siguiente sección.

La respuesta del algoritmo se presenta en la gráfica:

En el algoritmo se compara la solución numérica yi[n] con la solución analítica yia[n], obteniendo un error del orden 10-3 dado por el truncamiento de dígitos en ambos métodos.

 error numerico vs analítico: 
 maximo de |error|:  0.006000000000000227
 errado: 
[-0.          0.         -0.          0.         -0.          0.
 -0.          0.         -0.          0.         -0.         -0.
 -0.          0.          0.         -0.006      -0.0058     -0.00509
 -0.0041445  -0.00334173 -0.00267831 -0.0021442  -0.00171569 -0.00137264
 -0.00109813 -0.00087851 -0.00070281 -0.00056225 -0.0004498  -0.00035984
 -0.00028787]

El algoritmo continúa como un anexo a los algoritmos de respuesta a entrada cero e impulso.

Instrucciones en Python

# Sistema LTID. Respuesta entrada cero, 
# Respuesta a impulso
# resultado con raices Reales
# ejercicio lathi ejemplo 3.121 p286
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
# coeficientes E con grado descendente
Q = [1., -0.6, -0.16]
P = [5.,  0,    0   ]
# condiciones iniciales con n ascendente -2,-1, 0
inicial = [25/4, 0.] 

# PROCEDIMIENTO
# raices
gamma = np.roots(Q)

# revisa raices numeros complejos
revisaImag = np.iscomplex(gamma)
escomplejo = np.sum(revisaImag)
if escomplejo == 0:
    
    # coeficientes de ecuacion
    m = len(Q)-1
    A = np.zeros(shape=(m,m),dtype=float)
    for i in range(0,m,1):
        for j in range(0,m,1):
            A[i,j] = gamma[j]**(-m+i)
    ci = np.linalg.solve(A,inicial)

    # ecuacion yn
    n = sym.Symbol('n')
    y0 = 0
    for i in range(0,m,1):
        y0 = y0 + ci[i]*(gamma[i]**n)

    # Respuesta impulso
    ni = np.arange(-m,m,1)
    hi = np.zeros(m, dtype=float)
    xi = np.zeros(2*m, dtype=float)
    xi[m] = 1 # impulso en n=0

    # h[n] iterativo
    p_n = len(P)
    for i in range(0,m,1):
        for k in range(m,2*m,1):
            hi[i] = hi[i] - Q[k-m]*hi[m-k+1]
        for k in range(0,p_n,1):
            hi[i] = hi[i] + P[k]*xi[m+i]
    Bh = hi[:m]

    # coeficientes de ecuacion
    m = len(Q)-1
    Ah = np.zeros(shape=(m,m),dtype=float)
    for i in range(0,m,1):
        for j in range(0,m,1):
            Ah[i,j] = gamma[j]**(i)
    ch = np.linalg.solve(Ah,Bh)

    # ecuacion hn
    hn = 0
    for i in range(0,m,1):
        hn = hn + ch[i]*(gamma[i]**n)

    # Respuesta de estado cero ----------------
    # Rango [a,b], simetrico a 0
    b = 15 ; a = -b
    
    u = lambda n: np.piecewise(n,[n>=0,n<0],[1,0])
    
    x = lambda n: (0.25)**(n)*u(n)
    hL = sym.lambdify(n,hn,'numpy')
    h = lambda n: hL(n)*u(n)
    
    # PROCEDIMIENTO
    ni = np.arange(a,b+1,1)
    xi = x(ni)
    hi = h(ni)
    
    # Suma de Convolucion x[n]*h[n]
    yi = np.convolve(xi,hi,'same')
    
    # compara respuesta numerica con resultado analitico
    ya = lambda n: (-1.26*(0.25**n)+0.444*((-0.2)**n)+5.81*(0.8**n))*u(n)
    yia = ya(ni)

    errado = yia-yi
    erradomax = np.max(np.abs(errado))

        
# SALIDA
print('respuesta entrada cero: ')
print('raices: ', gamma)
if escomplejo == 0:
    print('Matriz: ')
    print(A)
    print('Ci:     ', ci)
    print('yn:  ')
    sym.pprint(y0)
    
    print('\n respuesta impulso: ')
    print('hi:',hi)
    print('Matriz: ')
    print(Ah)
    print('Ch: ',ch)
    print('n>=0, hn:  ')
    sym.pprint(hn)

    print('\n error numerico vs analitico: ')
    print(' maximo de abs(error): ', erradomax)
    print(' errado: ')
    print(errado)

    # SALIDA - GRAFICA
    plt.figure(1)
    plt.suptitle('Suma de Convolucion x[n]*h[n]')

    plt.subplot(311)
    plt.stem(ni,xi, linefmt='b--',
             markerfmt='bo',basefmt='k-')
    plt.ylabel('x[n]')

    plt.subplot(312)
    plt.stem(ni,hi, linefmt='b--',
             markerfmt='ro',basefmt='k-')
    plt.ylabel('h[n]')

    plt.subplot(313)
    plt.stem(ni,yi, linefmt='g-.',
             markerfmt='mo', basefmt='k-')
    plt.stem(ni,yia, linefmt='y-.',
             markerfmt='mD', basefmt='k-')

    plt.ylabel('x[n]*h[n]')
    plt.xlabel('n')

    plt.show()
    
else:
    print(' existen raices con números complejos.')
    print(' usar algoritmo de la sección correspondiente.')

6.4 LTI DT – Respuesta estado cero – Concepto

Referencia: Lathi 3.8 p280

Es la respuesta de un sistema cuando el estado es cero. se incia con que la entrada x[n] es la suma de componentes δ[n]

x[n] = x[0] \delta[n] + x[1] \delta[n-1] + x[2] \delta[n-2] + \text{...} + x[-1] \delta[n+1] + x[-2] \delta[n+2] + \text{...} x[n] = \sum_{m=-\infty}^{\infty} x[m] \delta[n-m]

En un sistema lineal, conociendo su respuesta impulso, responde a cualquier entrada como la suma de sus componentes.

y[n] = \sum_{m=-\infty}^{\infty} x[m] h[n-m] = x[n] \circledast h[n]

Propiedades de la sumatoria de convolución

Conmutativa

x_1[n] \circledast x_2[n] = x_2[n] \circledast x_1[n]

Distributiva

x_1[n] \circledast (x_2[n]+x_3[n])= (x_1[n] \circledast x_2[n]) + (x_1[n] \circledast x_3[n])

Asociativa

x_1[n] \circledast (x_2[n] \circledast x_3[n]) = (x_1[n] \circledast x_2[n]) \circledast x_3[n]

Desplazamiento

Si

x_1[n] \circledast x_2[n] = c[n]

entonces

x_1[n-m] \circledast x_2[n-p] = c[n-m-p]

Convolución de un impulso

x[n] \circledast \delta[n] = x[n]

Ancho o intervalo

Si x1[n] y x2[n] tienen anchos finitos W1 y W2, el ancho de la convolución entre ellos es W1+W2.

 

6.3.1 LTI DT – Respuesta impulso h[n]. Ejercicio con Python

Referencia: Lathi Ejemplo 3.17, 3.18  y 3.19 p277-278

Determine la respuesta al impulso unitario de un sistema LTID descrito por la ecuación de diferencias mostrada,

y[n+2] - 0.6 y[n+1] - 0.16 y[n] = 5x[n+2]

Para el diagrama se usa la ecuación desplazada en dos unidades,

y[n] - 0.6 y[n-1] - 0.16 y[n-2] = 5x[n] y[n] = 0.6 y[n-1] + 0.16 y[n-2] + 5x[n]

En la entrada se usa un impulso unitario x[n] = δ[n]

Desarrollo Analítico

Se requieren al menos dos valores iniciales en el ejercicio, por lo que se emplea la ecuación en forma iterativa, recordando que x[n] = δ[n], siendo entonces y[n] = h[n]

y[n] - 0.6 y[n-1] - 0.16 y[n-2] = 5x[n] h[n] - 0.6 h[n-1] - 0.16 h[n-2] = 5\delta [n]

aplicando para n=0

h[0] - 0.6 h[-1] - 0.16 h[-2] = 5\delta[0] h[0] - 0.6 (0) - 0.16 (0) = 5(1) h[0] = 5

luego se aplica para n = 1

h[1] - 0.6 h[0] - 0.16 h[-1] = 5 \delta[1] h[1] - 0.6 (5) - 0.16 (0) = 5(0) h[1] = 3

El ejercicio es continuación del primer ejercicio de respuesta a entrada cero, por lo que simplificamos el desarrollo como

(E^2 - 0.6 E - 0.16) y[n] = 5 E^2 x[n] \gamma^2 - 0.6 \gamma - 0.16 = 0

con raíces características γ1=-0.2 y γ2=0.8,

con los modos característicos se tiene que

y_c[n] = c_1 (-0.2)^n + c_2 (0.8)^n

que se convierte a:

h[n] = [c_1 (-0.2)^n + c_2 (0.8)^n ] \mu [n]

Para determinar c1 y c2 se recurre a encontrar dos valores de forma iterativa, con n=0 y n=1.

h[0] = [c_1 (-0.2)^0 + c_2 (0.8)^0 ] \mu [0] h[0] = c_1 + c_2 h[1] = [c_1 (-0.2)^1 + c_2 (0.8)^1 ] \mu [1] h[1] = -0.2 c_1 + 0.8 c_2

los valores de h[0]=5 y h[1]=3 se encontraron de forma iterativa

5 = c_1 + c_2 3 = -0.2 c_1 + 0.8 c_2

resolviendo con Python:

>>> import numpy as np
>>> A = [[ 1., 1.],
         [-0.2, 0.8]]
>>> B = [5.,3.]
>>> np.linalg.solve(A,B)
array([1., 4.])
>>> 

encontrando que c1=1 y c2=4

teniendo como resultado final:

h[n] = [ (-0.2)^n +4 (0.8)^n ] \mu [n]


Algoritmo en Python

La respuesta para el procedimiento anterior realizada con Python es:

 respuesta impulso: 
hi: [5. 3.]
Matriz: 
[[ 1.   1. ]
 [ 0.8 -0.2]]
Ch:  [4. 1.]
n>=0, hn:  
        n          n
1.0*-0.2  + 4.0*0.8 
>>>

se desarrolla el algoritmo a partir de entrada cero, continuando con el cálculo iterativo de h[n] para los valores iniciales, luego determina los valores de los coeficientes de la ecuación h[n]

# Sistema LTID. Respuesta impulso
# QE con raices Reales NO repetidas Lathi ejemplo 3.13 pdf271
# Lathi ejemplo 3.18 y 3.19 pdf278,
# blog.espol.edu.ec/telg1001
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
# coeficientes E con grado descendente
QE = [1., -0.6, -0.16]
PE = [5.,  0,    0   ]
# condiciones iniciales ascendente ...,y[-2],y[-1]
inicial = [25/4, 0.]

tolera = 1e-6  # casi_cero
muestras = 10  # para grafica

# PROCEDIMIENTO
# Respuesta a ENTRADA CERO
# raices, revisa numeros complejos
gamma = np.roots(QE)
revisaImag = np.iscomplex(gamma)
escomplejo = np.sum(revisaImag)

# coeficientes de ecuacion
m_q = len(QE)-1
Ac = np.zeros(shape=(m_q,m_q),dtype=float)

# revisa si parte compleja <tolera o casi_cero 
if escomplejo>0:
    for i in range(0,m_q,1):
        valorimag = np.imag(gamma[i])
        if np.abs(valorimag)<tolera:
            gamma[i] = float(np.real(gamma[i]))
    sumaj = np.sum(np.abs(np.imag(gamma)))
    if sumaj <tolera:
        print(sumaj)
        gamma = np.real(gamma)
        escomplejo = 0

# revisa repetidos
unicoscuenta = np.unique(gamma,return_counts=True)
repetidas = np.sum(unicoscuenta[1]-1)

# Respuesta impulso h[n]
ki = np.arange(-m_q,m_q,1)
hi = np.zeros(m_q, dtype=float)
xi = np.zeros(2*m_q, dtype=float)
xi[m_q] = 1 # impulso en n=0
Ah = np.zeros(shape=(m_q,m_q),dtype=float)

# h[n] iterativo
p_n = len(PE)
for i in range(0,m_q,1):
    for k in range(m_q,2*m_q,1):
        hi[i] = hi[i] - QE[k-m_q]*hi[m_q-k+1]
    for k in range(0,p_n,1):
        hi[i] = hi[i] + PE[k]*xi[m_q+i]
Bh = np.copy(hi[:m_q])

# coeficientes de ecuacion
for i in range(0,m_q,1):
    for j in range(0,m_q,1):
        Ah[i,j] = gamma[j]**(i)
ch = np.linalg.solve(Ah,Bh)

# ecuacion hn
n = sym.Symbol('n')
hn = 0*n
for i in range(0,m_q,1):
    hn = hn + ch[i]*(gamma[i]**n)

# SALIDA
if escomplejo == 0: 
    print('\n respuesta impulso: ')
    print('Bh:',Bh)
    print('Matriz: ')
    print(Ah)
    print('Ch: ',ch)
    print('n>=0, hn:  ')
    sym.pprint(hn)
else:
    print(' existen raices con números complejos.')
    print(' usar algoritmo de la sección correspondiente.')
    

# grafica datos
ki = np.arange(-m_q,muestras,1)
hi = np.zeros(muestras+m_q)

if escomplejo == 0: 
    # evaluación de h[n]
    h_n = sym.lambdify(n,hn)
    hi[m_q:] = h_n(ki[m_q:])

    # grafica h[n]
    plt.stem(ki,hi,label='h[n]',
             markerfmt ='C1o',
             linefmt='C2--')
    plt.legend()
    plt.grid()
    plt.ylabel('h[n]')
    plt.xlabel('ki')
    plt.title('h[n] = '+str(hn))
    plt.show()

Comentario: Aunque es relativamente simple determinar la respuesta al impulso h[n] usando el método descrito, el desarrollo del tema se realiza mejor en la unidad con Transformada z. Lathi p280.