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