7.1.3 Transformada z unilateral e inversa con Sympy-Python

Referencia: Tabla de transformadas z: Lathi Tabla 5.1 Transformada z p492. Oppenheim tabla 10.2 p776, Schaum Hsu Tabla 4-1 p170. Propiedades:  Schaum Hsu Tabla 4-2 p173. Lathi Tabla 5.2 p509. Oppenheim Tabla 10.3 p793

Usando los algoritmos presentados para transformadas z a partir de la tabla y propiedades, se incorpora en la instrucción z_transform(f,n,z) incorporada en telg1001.py. Para demostración se desarrollan algunos ejemplos que se desarrollan a partir de la tabla de transformadas z.

Transformada z de impulso unitario δ(t)

>>> import sympy as sym
>>> import telg1001_z as fcnm
>>> n = sym.Symbol('n', real=True)
>>> z = sym.Symbol('z')
>>> d = sym.DiracDelta(n)
>>> f = d
>>> fcnm.z_transform(f,n,z)
1
>>>

Transformada z de impulso unitario desplazado δ(t-2) y δ(t+2)

>>> fcnm.z_transform(f.subs(n,n-2),n,z)
z**(-2)
>>> fcnm.z_transform(f.subs(n,n+2),n,z)
z**2
>>>

Tabla de Transformadas z unilateral e Inversas usando algoritmo en telg1001.py

Para comprobar las transformadas de los algoritmos presentados e integrados en telg1001.py, se crea una tabla con expresiones f[n] de prueba en el dominio n de tiempo discreto.

Al resultado se aplica la inversa de F[z] que muestra que la expresión es invertible usando la transformada z unilateral.

El resultado del algoritmo se muestra a continuación y luego el algoritmo.

z_transform :  0
  f[n]: DiracDelta(n)
  F[z]: (1, True, True)
  f[n]: (DiracDelta(n), True, True)

z_transform :  1
  f[n]: DiracDelta(n - 2)
  F[z]: (z**(-2), Abs(z) > 0, True)
  f[n]: (DiracDelta(n - 2), True, True)

z_transform :  2
  f[n]: DiracDelta(n + 2)
  F[z]: (z**2, True, True)
  f[n]: (DiracDelta(n + 2), True, True)

z_transform :  3
  f[n]: 3*DiracDelta(n + 2)
  F[z]: (3*z**2, True, True)
  f[n]: (3*DiracDelta(n + 2), True, True)

z_transform :  4
  f[n]: -3*DiracDelta(n + 2)
  F[z]: (-3*z**2, True, True)
  f[n]: (-3*DiracDelta(n + 2), True, True)

z_transform :  5
  f[n]: Heaviside(n)
  F[z]: (z/(z - 1), Abs(z) > 1, True)
  f[n]: (Heaviside(n), Abs(z) > 1, True)

z_transform :  6
  f[n]: Heaviside(n - 2)
  F[z]: (1/(z*(z - 1)), Abs(z) > 1, True)
  f[n]: [None, [1/(z*(z - 1))]]

z_transform :  7
  f[n]: Heaviside(n + 2)
  F[z]: (z**3/(z - 1), Abs(z) > 1, True)
  f[n]: (Heaviside(n + 2), Abs(z) > 1, True)

z_transform :  8
  f[n]: 3*Heaviside(n + 2)
  F[z]: (3*z**3/(z - 1), Abs(z) > 1, True)
  f[n]: (3*Heaviside(n + 2), Abs(z) > 1, True)

z_transform :  9
  f[n]: n*Heaviside(n)
  F[z]: (z/(z - 1)**2, Abs(z) > 1, True)
  f[n]: (n*Heaviside(n), Abs(z) > 1, True)

z_transform :  10
  f[n]: n**2*Heaviside(n)
  F[z]: (z*(z + 1)/(z - 1)**3, Abs(z) > 1, True)
  f[n]: (n**2*Heaviside(n), Abs(z) > 1, True)

z_transform :  11
  f[n]: n**3*Heaviside(n)
  F[z]: (z*(z**2 + 4*z + 1)/(z - 1)**4, Abs(z) > 1, True)
  f[n]: (n**3*Heaviside(n), Abs(z) > 1, True)

z_transform :  12
  f[n]: 0.5**n*Heaviside(n)
  F[z]: (z/(z - 0.5), Abs(z) > 0.5, True)
  f[n]: (0.5**n*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  13
  f[n]: 0.5**n*Heaviside(n)
  F[z]: (z/(z - 0.5), Abs(z) > 0.5, True)
  f[n]: (0.5**n*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  14
  f[n]: Heaviside(n)/4**n
  F[z]: (z/(z - 1/4), Abs(z) > 1/4, True)
  f[n]: (Heaviside(n)/4**n, Abs(z) > 1/4, True)

z_transform :  15
  f[n]: -0.5**n*Heaviside(-n - 1)
  F[z]: (z/(z - 0.5), Abs(z) < 0.5, True) f[n]: (0.5**n*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  16
  f[n]: 0.5**(n - 1)*Heaviside(n - 1)
  F[z]: (1.0/(z - 0.5), Abs(z) > 0.5, True)
  f[n]: [None, [1/(z - 0.5)]]

z_transform :  17
  f[n]: 4**(1 - n)*Heaviside(n - 1)
  F[z]: (1/(z - 1/4), Abs(z) > 1/4, True)
  f[n]: [None, [1/(z - 1/4)]]

z_transform :  18
  f[n]: 0.5**n*n*Heaviside(n)
  F[z]: (0.5*z/(z - 0.5)**2, Abs(z) > 0.5, True)
  f[n]: (0.5**n*n*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  19
  f[n]: -0.5**n*n*Heaviside(-n - 1)
  F[z]: (0.5*z/(z - 0.5)**2, Abs(z) < 0.5, True) f[n]: (0.5**n*n*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  20
  f[n]: 0.5**n*(n + 1)*Heaviside(n)
  F[z]: z/(z - 0.5) + 0.5*z/(z - 0.5)**2
  f[n]: (2*0.5**(n + 1)*(n + 1)*Heaviside(n + 1), Abs(z) > 0.5, True)

z_transform :  21
  f[n]: 0.5**n*n**2*Heaviside(n)
  F[z]: (0.5*z*(z + 0.5)/(z - 0.5)**3, Abs(z) > 0.5, True)
  f[n]: (0.5**n*n**2*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  22
  f[n]: 1.33333333333333*0.5**n*n*(n - 2)*(n - 1)*Heaviside(n)
  F[z]: 1.33333333333333*z/(z - 0.5)**2 - 2*z*(z + 0.5)/(z - 0.5)**3 + 0.666666666666667*z*(z**2 + 2*z + 0.25)/(z - 0.5)**4
  f[n]: (1.33333333333333*0.5**n*n*(n - 2)*(n - 1)*Heaviside(n), Abs(z) > 1, 0)

z_transform :  23
  f[n]: cos(2*n)*Heaviside(n)
  F[z]: (z*(z - cos(2))/(z**2 - 2*z*cos(2) + 1), Abs(z) > 1, True)
  f[n]: (cos(2*n)*Heaviside(n), Abs(z) > 1, True)

z_transform :  24
  f[n]: 0.5**n*cos(2*n)*Heaviside(n)
  F[z]: (z*(z - 0.5*cos(2))/(z**2 - z*cos(2) + 0.25), Abs(z) > 0.5, True)
  f[n]: (0.5**n*cos(2*n)*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  25
  f[n]: sin(2*n)*Heaviside(n)
  F[z]: (z*sin(2)/(z**2 - 2*z*cos(2) + 1), Abs(z) > 1, True)
  f[n]: (sin(2*n)*Heaviside(n), Abs(z) > 1, True)

z_transform :  26
  f[n]: 0.5**n*sin(2*n)*Heaviside(n)
  F[z]: (0.5*z*sin(2)/(z**2 - z*cos(2) + 0.25), Abs(z) > 0.5, True)
  f[n]: (0.5**n*sin(2*n)*Heaviside(n), Abs(z) > 0.5, True)

z_transform :  27
  f[n]: cos(2*n + 0.25)*Heaviside(n)
  F[z]: (z*(z*cos(0.25) - cos(1.75))/(z**2 - 2*z*cos(2) + 1), Abs(z) > 1, True)
  f[n]: (1.0*cos(2*n + 0.25)*Heaviside(n), Abs(z) > 1, True)

z_transform :  28
  f[n]: 3*0.5**n*cos(2*n + 0.25)*Heaviside(n)
  F[z]: (3*z*(z*cos(0.25) - 0.5*cos(1.75))/(z**2 - z*cos(2) + 0.25), Abs(z) > 0.5, True)
  f[n]: (3.0*0.5**n*cos(2*n + 0.25)*Heaviside(n), Abs(z) > 0.5, True)

Algoritmo en Python

Las intrucciones se pueden modificar para observar una expresión en detalle si se activa la variable sym.SYMPY_DEBUG=True. Por la extensión del resultado se recomienda activarla para una expresión.

import sympy as sym
import telg1001 as fcnm
sym.SYMPY_DEBUG = False

# INGRESO
n = sym.Symbol('n')
z = sym.Symbol('z')
u = sym.Heaviside(n)
d = sym.DiracDelta(n)

a0 = sym.Rational(1,4).limit_denominator(100)

tabla = [
           d,
           d.subs(n,n-2),
           d.subs(n,n+2),
           3*d.subs(n,n+2),
           -3*d.subs(n,n+2),
           u,
           u.subs(n,n-2),
           u.subs(n,n+2),
           3*u.subs(n,n+2),
           n*u,
           (n**2)*u,
           (n**3)*u,
           0.5**n*u,((1/2)**n)*u,(a0**n)*u,
           -(0.5**n)*u.subs(n,-n-1),
           (0.5**(n-1))*u.subs(n,n-1),
           (a0**(n-1))*u.subs(n,n-1),
           n*(0.5**n)*u,
           -n*(0.5**n)*u.subs(n,-n-1),
           (n+1)*(0.5**n)*u,
           n**2*(0.5**n)*u,
           n*(n-1)*(n-2)*(0.5**n)*u/((0.5**3)*6),
           sym.cos(2*n)*u,
           (0.5**n)*sym.cos(2*n)*u,
           sym.sin(2*n)*u,
           (0.5**n)*sym.sin(2*n)*u,
           sym.cos(2*n+0.25)*u,
           3*(0.5**n)*sym.cos(2*n+0.25)*u
           ]
# seleccionar revisar
revisar = tabla[:]

# revisa detalle de una transformada y la inversa
if len(revisar)<2:
    sym.SYMPY_DEBUG = True
else:
    sym.SYMPY_DEBUG = False

# PROCEDIMIENTO
for i in range(0,len(revisar),1):
    f = revisar[i]
    # SALIDA
    print('\nz_transform : ', i)
    print('  f[n]:',f)
    Fz = fcnm.z_transform(f,n,z)
    print('  F[z]:',Fz)
    if type(Fz)==tuple:
        fn = fcnm.inverse_z_transform(Fz[0],z,n)
    else:
        fn = fcnm.inverse_z_transform(Fz,z,n)
    print('  f[n]:',fn)

7.1.2 Transformada z – Propiedades con Sympy-Python

Luego de intentar obtener la Transformada z con la Tabla de Pares f[n] y F[z], si del resultado sigue siendo None se busca aplicar las reglas descritas en la tabla de propiedades de la Transformada z.

Referencia: Schaum Hsu Tabla 4-2 p173. Lathi Tabla 5.2 p509. Oppenheim Tabla 10.3 p793


1. Propiedad de multiplicación por n de la transformada z

Presentada en la tabla de propiedades y también conocida como Diferenciación en z, se presenta como:

n x[n] u[n] -z \frac{\delta}{\delta z}X(z) R’ = R

Se aplica sobre la señal x[n]=1

x[n] = n \mu[n]

A partir de la tabla de transformadas se tiene que la expresión para μ[n] es

F[z] = \frac{z}{z-1} = (z)\frac{1}{z-1} \frac{\delta}{\delta z}F[z] = \frac{1}{z-1}-\frac{z}{(z-1)^2}

por lo que para X[z] de obtiene como:

X[z] = -z \frac{\delta}{\delta z} F[z] = -z \Big[ \frac{1}{z-1}-\frac{z}{(z-1)^2}\Big] = -z \Big[ -\frac{1}{(z-1)^2}\Big] X[z] = \frac{z}{(z-1)^2}

Algoritmo en Python

El desarrollo del algoritmo es posterior a la búsqueda en la tabla de pares de transformadas z, por lo que se importa el proceso desde telg1001.py. En la tabla, para comprobar el funcionamiento de la propiedad con los resultados del texto, se usa una tabla básica, que luego puede ser ampliada.

El algoritmo se desarrolla para una expresión básica como la de la tabla. Luego los procedimientos se integran en una sola instrucción tal como z_transform(f,n,z) que usa primero la tabla, luego revisa las propiedades y fórmula de la transformada para mostrar un resultado.

u = sym.Heaviside(n)

#f = u
f = n*u

Resultados con el algoritmo:

 k ; func: 1  ;  n*Heaviside(n)
no hubo expresión similar en tabla de transformadas
buscando aplicar una propiedad n*f[n]
  f_powna ; n**a: n  ;  {a_: 1}

 _z_propiedad n*f[n] -------------
  n**a ; func/n, : n  ;  Heaviside(n)
_z_pairs_table match:
  f:           Heaviside(n)
  z_pair f[n]: Heaviside(n)
  z_pair F[z]: z/(z - 1)
  try,check: True -> True
  _z_propiedad n*f[n] <--> -z*diff(Fz):
   (z/(z - 1)**2, Abs(z) > 1, True)
>>> 

Instrucciones en Python

El coeficiente constante k se separa de la expresión para simplificar el análisis de la expresión en func. La constante se reincorpora al resultado al final.

El procedimiento con la tabla de pares de transformadas  z_pairs_properties se obtiene desde telg1001.py , se usa con el parámetro apply_properties=False para mostrar el funcionamiento del algoritmo para el ejercicio. la instrucción z_pairs_properties hace el llamado a una función que analiza la expresión con las propiedades ya implementadas y se omite su escritura en condiciones normales.

Las expresiones pueden tener dos factores con exponentes: an, na, que corresponden a dos propiedades diferentes. El factor  de interés para el ejercicio es na.  La búsqueda de éstos factores se realiza con las instrucciones sym.pow(n,a) y sym.pow(a,n), en la variable f_powna que se inicializan con sym.S.One.

Cuando se encuentra la expresión con f_powna se modifica expresión func_n =func/n, con la que se vuelve a buscar el resultado en la tabla. El proceso para analizar la expresión puede realizarse de forma recursiva hasta obtener una expresión simple como μ[n] y regresar al estado anterior para completar las operaciones. La forma recursiva requiere implementar el procedimiento con def z_propiedades().

La aplicación de la propiedad de multiplicación para n(x[n]) prepara la expresión para aplicar:

n x[n] \longleftrightarrow -z \frac{\delta}{\delta z}X(z)
Fz = (k*(-z)*sym.factor(sym.diff(Fz[0],z,1)),
                  Fz[1], Fz[2])

que contiene las partes de la transformada como el plano de convergencia y condiciones de aplicación que también deben tomarse en cuenta durante el proceso. Aunque para expresiones se sumas, no se toman en cuenta para el resultado final como en los ejercicios de los textos de referencia.

Esta sección expone cómo se implementa la propiedad de transformada z, el resultado final se encuentra en el archivo.py.

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

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

#f = u
f = n*u
#f = sym.cos(2*n)
#f = sym.sin(2*n)

# PROCEDIMIENTO
f = sym.expand(sym.powsimp(f))
a = sym.Wild('a', exclude=[n])
b = sym.Wild('b', exclude=[n])
y = sym.Wild('y')
g = sym.WildFunction('g', nargs=1)

# separa constantes del término
k, func = f.as_independent(n, as_Add=False)
print(' k ; func:',k,' ; ',func)
Fz = z_pairs_properties(func, n, z, apply_properties=False)
if not(Fz==None):
    Fz = k*Fz
    print(' usando tabla de pares: \n ',Fz)

if Fz==None: # no encontrada en tabla
    print('no hubo expresión similar en tabla de transformadas')
    print('buscando aplicar una propiedad n*f[n]')
    # crear una función para aplicar varias veces las propiedades
    # en las expresiones, en forma recursiva.
    
    # busca factores pow(n,a) or pow(a,n)
    f_powna = sym.S.One ; f_powan = sym.S.One
    factor_Mul = sym.Mul.make_args(func)
    for factor_k in factor_Mul:
        ma_powna = factor_k.match(sym.Pow(n,a))
        ma_powan = factor_k.match(sym.Pow(a,n))
        if ma_powna and not(ma_powna[a]==sym.S.Zero) and ma_powna[a].is_integer:
            f_powna = f_powna*factor_k
            print('  f_powna ; n**a:',factor_k,' ; ',ma_powna)
        if ma_powan:
            f_powan = f_powan*factor_k
            print('  f_powan ; a**n:',factor_k,' ; ',ma_powan)

    # otras formas de expresión a revisar
    ma_un = func.match(sym.Heaviside(y))
    ma_gn = func.match(g)
    ma_gu = func.match(g*sym.Heaviside(y))

    # z_propiedad de multipliación n*f[n] <--> -z*dF(z)/dz --------
    if not(f_powna==sym.S.One): # n**a factor encontrado
        ma_powna = f_powna.match(sym.Pow(n,a))
        func_n = func/n # aplicar luego de forma recursiva
        print('\n _z_propiedad n*f[n] -------------')
        print('  n**a ; func/n, :',n,' ; ',func_n)
        Fz = z_pairs_properties(func_n,n,z,
                                apply_properties=False)
        if not(Fz==None):
            Fz = (k*(-z)*sym.factor(sym.diff(Fz[0],z,1)),
                  Fz[1], Fz[2])
            print('  _z_propiedad n*f[n] <--> -z*diff(Fz):\n  ',Fz)

if Fz==None: # no encontrada en tabla, tampoco con la propiedad
    print('propiedad implementada no es suficiente')
    print('implementar otras propiedades.')

Un procedimiento semejante se implementa para la transformada z inversa.

2. Propiedad de multiplicación por n de la transformada z inversa

Se aplica a partir de F(z), donde la expresión será tipo polinomio, por lo que se prepara la expresión para conocer el grado del numerador y denominador, separar su signo, constante de multiplicación, coeficientes, etc.

Para el resultado del ejercicio anterior f(n) = nμ(n) se obtuvo la expresión:

F[z] = \frac{z}{(z-1)^2}
n x[n] u[n] -z \frac{\delta}{\delta z}X(z) R’ = R

el sentido inverso de la propiedad, se aplica con un integral sobre la expresión con denominador Q=(z-a)b. y a partir de la respuesta F[z] con al menos el grado del numerador P de 1.

El integral se desarrolla con la expresión despejada de la derivada y con un valor F(0)=0 para determinar la constante del integral.

X[z] = \int{\Big[\frac{z}{(z-1)^2}\Big]\Big[\frac{-1}{z}\Big]} + C X[z] = \frac{1}{(z-1)} + C

usando el valor de para la constante F(0)=0

0 = \frac{1}{(0-1)} + C C = 1 X[z] = \frac{1}{(z-1)} +1 = \frac{1+(z-1)}{(z-1)} X[z] = \frac{z}{(z-1)}

de la tabla la transformada z inversa es μ(n), y al aplicar al resultado la propiedad se tiene:

x[n] = n \mu [n]

resultados con el algoritmo:

 Parametros de P/Q -------------
 P: z
 Q: (z - 1)**2
 P_signo ;k :  1  ;  1
 P_grado, P_ceros: 1  ;  {0: 1}
 Q_grado, Q_polos (real):  2  ;  [1, 1]
 ma_P1 (a*z+ b) :  {a_: 1, b_: 0}
 ma_Q1 (z-a)**b :  {a_: 1, b_: 2}
 ma_Q2 a*z**2+ b*z + r**2:  None
 Fz :  z/(z - 1)**2

 _z_propiedad multiplicación nf[n] <--> -z*diff(F[z])-- 
 Func = integrate(factor(Fz)/z,z):
  z/(z - 1)
_z_pairs_table match:
  k ; F:      1 z/(z - 1)
  z_pair F[z]: z/(z - 1)
         ma_z: {}
  z_pair f[n]: Heaviside(n)
  try,check: True -> True
  _z_propiedad multiplicación nf[n]:
  (n*Heaviside(n), Abs(z) > 1, True)
>>> 

Instrucciones en Python

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

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

# multiplicación nf[n] <--> -z*diff(F[z])
#F = z/(z - 1) #f = u
F = z/(z - 1)**2 #f = n*u
#F = z*(z + 1)/(z - 1)**3 # (n**2)*u

# PROCEDIMIENTO
# separa constantes del término
fn = z_pairs_prop_inverse(F, z, n, apply_properties=False)
if not(fn==None):
    print(' usando tabla de pares: \n ',fn)

if fn==None: # no encontrada en tabla
    Fz = sym.simplify(F) #sym.simplify(F)
    a = sym.Wild('a', exclude=[n,z])
    b = sym.Wild('b', exclude=[n,z])
    r = sym.Wild('r', exclude=[n,z])
    y = sym.Wild('y')
    g = sym.WildFunction('g', nargs=1)
    # analiza como polinomio F[z]:
    # P_signo, constante, F[z] pares
    [P,Q] = F.as_numer_denom()
    # P
    P = sym.Poly(P,z)
    ma_P1 = P.match(a*z+ b)
    P_zeros  = sym.roots(P)
    P_degree = sym.degree(P,z)
    P_leadcoef = sym.LC(P)
    k = sym.Abs(P_leadcoef)
    P_sign = P_leadcoef/k
    P = P/P_leadcoef
    # Q
    Q = sym.factor(Q)
    ma_Q1 = Q.match((z-a)**b)
    ma_Q2 = Q.match(a*z**2+ b*z + r**2)
    Q_poles = sym.real_roots(Q)
    Q_degree = sym.degree(Q,z)
    # separa constante
    Fz = sym.factor(Fz*P_sign/k)
    print('\n Parametros de P/Q -------------')
    print(' P:',P)
    print(' Q:',Q)
    print(' P_signo ;k : ',P_sign,' ; ',k)
    print(' P_grado, P_ceros:',P_degree,' ; ', P_zeros)
    print(' Q_grado, Q_polos (real): ',Q_degree,' ; ',Q_poles)
    print(' ma_P1 (a*z+ b) : ',ma_P1)
    print(' ma_Q1 (z-a)**b : ',ma_Q1)
    print(' ma_Q2 a*z**2+ b*z + r**2: ',ma_Q2)
    print(' Fz : ',Fz)
    
    # _z_property nf[n] <--> -z*diff(F[z])
    if ma_Q1 and ma_Q1[a]==1 and P_degree>0:
        FuncI = sym.factor(Fz/(-z),z)
        Func = sym.integrate(FuncI,z)
        F0 = 0 # para Constante de integral
        C = -Func.subs(z,0)+F0 
        FunC = sym.factor(Func+C)
        print('\n _z_propiedad multiplicación nf[n] <--> -z*diff(F[z])-- ')
        print(' Func = integrate(factor(Fz)/z,z):\n ',FunC)
        fn = z_pairs_prop_inverse(FunC, z, n, apply_properties=True)
        if not(fn==None):
            fn = (P_sign*k*n*fn[0],fn[1],fn[2])
            print('  _z_propiedad multiplicación nf[n]:\n ',fn)

if fn==None: # no encontrada en tabla, tampoco con la propiedad
    print('propiedad implementada no es suficiente')
    print('implementar otras propiedades.')

Si la propiedad debe aplicarse más de una vez, será necesario convertir el bloque a una función para hacer llamadas recursivas a si misma, como es el caso de n2*u

F = z*(z + 1)/(z - 1)**3 # (n**2)*u

o al final de la sección de la propiedad, usar apply_properties=True, para usar el algoritmo implementado en telg1001_z .

7.1.1 Transformada z – Pares f[n] y F[z] con match() de Sympy-Python

Para facilitar el uso de la Transformada z se usa la Tabla de Pares f[n] y F[z], donde se busca una expresión semejante que permite cambiar del dominio_n al dominio_z. Los Pares se usan en conjunto con la tabla de propiedades de la Transformada z con lo que se extiende las posibilidades de uso para la tabla.

Se propone primero usar la tabla de pares con Sympy para buscar las expresiones semejantes, usando la instrucción ‘f.match()‘, donde f es la función en el dominio_n y F es la función en el domino_z. En la siguiente sección se incorpora la tabla de propiedades, para en conjunto realizar una búsqueda más completa de las transformadas z o de sus inversas.

Referencia: Schaum Hsu Tabla 4-1 p170. Lathi Tabla 5.1 Transformada z p492. Oppenheim tabla 10.2 p776


Ejercicio 1. transformada z de f[n] = cos(2n) con Sympy y la instrucción f.match()

Revise si existe una expresión para f[n] en una pequeña lista de pares de transformadas:

f[n] = \cos [2n]

Algoritmo en Python

El algoritmo inicia con las variables a usar para cada dominio y la expresión de Heaviside como en los ejemplos de las unidades anteriores.

En la tabla se requiere unas variables tipo «comodín» (Wild) para buscar la semejanza de las expresiones, se usan como sym.Wild('a', exclude=[n]) de la que se excluye la variable de tiempo discreto n.

Se incluye una pequeña lista como tabla de pares de transformadas como un ejemplo básico, que luego puede ser ampliada. Los pares se incluyen como tuplas con las expresiones para cada domino.

Para encontrar semejanzas, se recorre cada par de transformadas, comparando las expresiones con f.match(par_nz[0]), que para revisar la expresión del domino n, solo toma la primera parte de la tupla. Si existe coincidencia, se crea un diccionario que indica los valores comodín que hacen que las expresiones sean iguales.

Resultado del algoritmo

 se encontró expresión similar:
  f:            cos(2*n)
  z_pares f[n]: cos(n*a_)
   similar con: {a_: 2}
  z_pares F[z]: (z**2 - z*cos(a_))/(z**2 - 2*z*cos(a_) + 1)
  Fz:         : (z**2 - z*cos(2))/(z**2 - 2*z*cos(2) + 1)
>>> 

Instrucciones en Python

# transformada z con f.match de Sympy
# expresiones similares o semejantes f.match
import sympy as sym

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

#f = u
#f = 5*u
f = sym.cos(2*n)
#f = sym.sin(2*n)

# PROCEDIMIENTO
# para revisar semejanza de expresiones fn y n_dom
a = sym.Wild('a', exclude=[n])
b = sym.Wild('b', exclude=[n])

# tabla de pares [(n_dom, z_dom)]
z_pares = [
    # impulso unitario d[n], DiracDelta
    (DiracDelta(n),
     S.One,
     S.true, S.Zero, dco),
    (DiracDelta(a*n),
     S.One,
     Abs(a)>0, S.Zero, dco),
    # escalon unitario u[n], Heaviside
    (Heaviside(n),
     z/(z-1),
     S.true, Abs(z) > 1, dco),
    # cos[n], sin[n] ,trigonometricas
    (cos(a*n),
     (z*(z-cos(a)))/(z**2-(2*cos(a))*z+1),
     Abs(a)>0, Abs(z) > 1, dco),
    (sin(a*n),
     (sin(a)*z+0)/(z**2-(2*cos(a))*z+1),
     Abs(a)>0, Abs(z) > 1, dco),
    ]

Fz = None; f_pares = None # sin similar
z_pares_len = len(z_pares) ; i=0
while i<z_pares_len and Fz==None:
    par_nz = z_pares[i]
    n_dom = par_nz[0]
    z_dom = par_nz[1]
    similar = f.match(n_dom)
    # entrega diccionario de expresion similar
    # si el diccionario es vacio, es coincidente
    if similar or similar=={}:
        f_pares = par_nz
        f_args = similar
        Fz_ = z_dom
        Fz  = z_dom.xreplace(similar)
    i = i+1 # siguiente par

# SALIDA
if not(Fz==None):
    print(' se encontró expresión similar:')
    print('  f:           ',f)
    print('  z_pares f[n]:',f_pares[0])
    print('   similar con:',f_args)
    print('  z_pares F[z]:',f_pares[1])
    print('  Fz:         :',Fz)
else:
    print(' NO se encontró una expresión similar...')

El concepto de expresiones similares se prueba con:

f[n] = \mu [n]

donde se observa que a respuesta es un diccionario vacío

 se encontró expresión similar:
  f:            Heaviside(n)
  z_pares f[n]: Heaviside(n)
   similar con: {}
  z_pares F[z]: z/(z - 1)
  Fz:         : z/(z - 1)
>>> 

mientras que para una función que no se encuentra en la tabla, el resultado debe ser Fz=None

f[n] = \sin [2n]
 NO se encontró una expresión similar...
>>> type(Fz)
<class 'NoneType'>
>>> 

El concepto básico del algoritmo se extiende para las tablas de Pares de transformadas que se adjuntan en un archivo de funciones telg1001.py.


Ejercicio 2. Transformada z inversa de F[z] = z/(z-1) con Sympy y la instrucción f.match()

Revisar si existe una expresión para F[z] en una pequeña lista de pares de transformadas:

F[z] = \frac{z}{z-1}

Algoritmo en Python

Semejante al ejercicio 1 para f[n] a F[z], el algoritmo inicia con las variables a usar para cada dominio y de F[z].

Para encontrar semejanzas, se recorre cada par de transformadas, comparando las expresiones con F.match(par_nz[1]), que para revisar la expresión del domino z, solo toma la segunda parte de la tupla. Si existe coincidencia, se crea un diccionario que indica los valores comodín que hacen que las expresiones sean iguales.

Resultado del algoritmo

 se encontró expresión similar:
  F:            z/(z - 1)
  z_pares F[z]: z/(z - 1)
   similar con: {}
  z_pares f[n]: Heaviside(n)
  fn:         : Heaviside(n)
>>> 

Instrucciones en Python

# transformada z inversa con F.match de Sympy
# expresiones similares o semejantes F.match
import sympy as sym

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

F = z/(z-1)
#F = (z**2 - z*sym.cos(2))/(z**2 - 2*z*sym.cos(2) + 1)

# PROCEDIMIENTO
# para revisar semejanza de expresiones fn y n_dom
a = sym.Wild('a', exclude=[n])
b = sym.Wild('b', exclude=[n])

# tabla de pares [(n_dom, z_dom)]
z_pares = [
    # impulso unitario d[n], DiracDelta
    (DiracDelta(n),
     S.One,
     S.true, S.Zero, dco),
    (DiracDelta(a*n),
     S.One,
     Abs(a)>0, S.Zero, dco),
    # escalon unitario u[n], Heaviside
    (Heaviside(n),
     z/(z-1),
     S.true, Abs(z) > 1, dco),
    # cos[n], sin[n] ,trigonometricas
    (cos(a*n),
     (z*(z-cos(a)))/(z**2-(2*cos(a))*z+1),
     Abs(a)>0, Abs(z) > 1, dco),
    (sin(a*n),
     (sin(a)*z+0)/(z**2-(2*cos(a))*z+1),
     Abs(a)>0, Abs(z) > 1, dco),
    ]

fn = None; F_pares = None # sin similar
z_pares_len = len(z_pares) ; i=0
while i<z_pares_len and fn==None:
    par_nz = z_pares[i]
    n_dom = par_nz[0]
    z_dom = par_nz[1]
    similar = F.match(z_dom)
    # entrega diccionario de expresion similar
    # si el diccionario es vacio, es coincidente
    if similar or similar=={}:
        F_pares = par_nz
        F_args = similar
        fn_ = n_dom
        fn  = n_dom.xreplace(similar)
    i = i+1 # siguiente par

# SALIDA
if not(fn==None):
    print(' se encontró expresión similar:')
    print('  F:           ',F)
    print('  z_pares F[z]:',F_pares[1])
    print('   similar con:',F_args)
    print('  z_pares f[n]:',F_pares[0])
    print('  fn:         :',fn)
else:
    print(' NO se encontró una expresión similar...')

otros ejercicios realizados con

F[z] =\frac{z^2 - zcos(2)}{z^2 - 2zcos(2) + 1}

resultado del algoritmo

 se encontró expresión similar:
  F:            (z**2 - z*cos(2))/(z**2 - 2*z*cos(2) + 1)
  z_pares F[z]: (z**2 - z*cos(a_))/(z**2 - 2*z*cos(a_) + 1)
   similar con: {a_: 2}
  z_pares f[n]: cos(n*a_)
  fn:         : cos(2*n)
>>> 

Referencias: Sympy: match(pattern, old=False) https://docs.sympy.org/latest/modules/core.html?highlight=match#sympy.core.basic.Basic.match

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', integer=True, 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:

suma de convolución

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.

diagrama de sistema transformada z


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:

suma de convolución

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]

diagrama sistema transformada z respuesta impulso

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]

grafica respuesta a impulso


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.