Algoritmos y funciones: telg1001.py

Resumen de algoritmos y funciones del curso señales y sistemas. Descargue una copia del archivo como telg1001.py en el directorio de trabajo, importe las funciones y use llamando con fcnm.funcion().

import telg1001 as fcnm

respuesta = fcnm.funcion(parametros)

El contenido de telg1001.py se actualiza cuando al desarrollar un ejercicio se encuentra que se puede mejorar una función manteniendo compatibilidad con lo realizado anteriormente. Se incorpora la referencia al ejercicio como comentario.

# Transformadas de Laplace, funciones
# http://blog.espol.edu.ec/telg1001/
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
s = sym.Symbol('s')
t = sym.Symbol('t',real=True)

# Transformada de Laplace para f(t) con Sympy-Python
# http://blog.espol.edu.ec/telg1001/transformada-de-laplace-para-ft-con-sympy-python/
def separa_constante(termino):
    ''' separa constante antes de usar
        sym.laplace_transform(term_suma,t,s)
        para incorporarla luego de la transformada
        inconveniente revisado en version 1.11.1
    '''
    constante = 1
    if termino.is_Mul:
        for term_mul in termino.args:
            if not(term_mul.has(t)):
                constante = term_mul
        termino = termino/constante
    return([termino,constante])

...

4.4.5 LTI Laplace – Circuitos electrónicos Activos «Op Amp»

Referencia: Lathi 4.4-1 p382,  Oppenheim 11.50 p896

Se amplian los conceptos de circuitos pasivos analizados con transformadas e Laplace, aplicando a circuitos activos. Se obtienen los circuitos equivalentes o modelos matemáticos y se repiten los procedimientos anteriores.

El elemento activo más conocido es el amplificador operacional  (op amp) que tienen ganancia «muy alta». El  voltaje de salida v2 = – A v1, donde A es del orden de 105 o 106. Un factor importante es que la impedancia de entrada es muy alta del orden de 1012Ω y la impedancia de salida es muy baja (50-100Ω)

La configuración de la ganancia se establece con los resistores Ra y Rb y la forma de conectar las entradas y salidas.

K = 1+\frac{R_a}{R_b} v_2 = K v_1 v_2 = (R_b + R_a) i_o = R_b i_o + R_a i_o v_2 = v_s = Ra i_o = R_a i_o \frac{v_2}{v_1} =\frac{R_b+R_a}{R_a} = 1+\frac{R_b}{R_a} = K

1. Amplificador Operacional en el dominio s

Referencia: Lathi 4.6-5 p399

Dada la alta impedancia del op amp, la corriente de retroalimentación  I(s) fluye solo por los resistores. El voltaje de entrada es cero o muy pequeño dada la ganancia muy grande del op amp. Dadas estas simplificaciones, se aproxima con mucha precisión que:

Y(s) = - I(s) Z_f(s) I(s) = \frac{X(s)}{Z(s)} Y(s) = -\frac{Z_f(s)}{Z(s)} H(s) = -\frac{Z_f(s)}{Z(s)}

1,1 Amplificador Operacional como multiplicador escalar

H(s) = -\frac{R_f}{R}

1,2 Amplificador Operacional como Integrador

Referencia: Oppenheim 11.52 p898

H(s) = \Big(-\frac{1}{RC}\Big) \frac{1}{s}

1,3 Amplificador Operacional como Sumador

Y(s) = - \Big[\frac{R_f}{R_1}X_1(s)+\frac{R_f}{R_2}X_2(s)+\frac{R_f}{R_3}X_3(s) \Big]

Observe que las ganancias del sumador son siempre negativas, hay una inversión de signo en cada señal de entrada.

Y(s) = K_1 X_1(s)+K_2 X_2(s)+K_3 X_3(s)


Ejemplo 1. Implementación con Op-Amp

Referencia: Lathi 4.25 p401

Realizar la implementación de un sistema dado por la función de transferencia H(s)

H(s) = \frac{2s+5}{s^2+4s+10}

El diagrama de bloques de la función de transferencia H(s) es,

Se agrupan algunos elementos como sumadores y sus factores de multiplicación. Para referencia se etiqueta cada punto como señal W(s) en cada punto donde el orden del exponente de ‘s’ es diferente.

Se considera la inversión de signo de la señal de entrada por la configuración del amplificador operacional y el factor K de cada rama a usar.

Se identifica el tipo de op amp a usar y se establecen los valores de los resistores en múltiplos de 10KΩ y los capacitores en el orden de 10 μF.


Ejemplo 1. Desarrollo analítico

Para revisar el comportamiento del circuito, en caso de implementarlo con OpAmps, el resultado de la función de transferencia para el impulso usando el algoritmo de la sección LTI CT Laplace – Ejercicio resuelto para Y(s)=H(s)X(s) con Sympy-Python

 H(s) = P(s)/Q(s):
   2*s + 5   
-------------
 2           
s  + 4*s + 10
 H(s) en factores:
   2*s + 5   
-------------
 2           
s  + 4*s + 10

H(s) parámetros cuadraticos:
(2*s + 5)/(s**2 + 4*s + 10) : 
 {'A': 2.0, 'B': 5.0, 'a': 2.0, 'c': 10.0,
  'r': 2.041241452319315, 'b': 2.449489742783178,
  'theta': -0.20135792079033082}

 h(t) :
/  ___  -2*t    /  ___  \                       \             
|\/ 6 *e    *sin\\/ 6 *t/      -2*t    /  ___  \|             
|------------------------ + 2*e    *cos\\/ 6 *t/|*Heaviside(t)
\           6                                   /             

polosceros:
Q_polos : {-2 - sqrt(6)*I: 1, -2 + sqrt(6)*I: 1}
P_ceros : {-5/2: 1}

Estabilidad de H(s):
 n_polos_real : 0
 n_polos_imag : 2
 enRHP : 0
 unicos : 0
 repetidos : 0
 asintota : estable

 X(s): 
1

Respuesta entrada cero ZIR H(s) y condiciones iniciales
term_cero : 0
ZIR :
0
yt_ZIR :
0

 ZSR respuesta estado cero:
ZSR :
   2*s + 5   
-------------
 2           
s  + 4*s + 10

 ZSR_Qs2 :
 (2*s + 5)/(s**2 + 4*s + 10) :
  {'A': 2.0, 'B': 5.0, 'a': 2.0, 'c': 10.0,
   'r': 2.041241452319315, 'b': 2.449489742783178,
   'theta': -0.20135792079033082}
yt_ZSR :
/  ___  -2*t    /  ___  \                       \             
|\/ 6 *e    *sin\\/ 6 *t/      -2*t    /  ___  \|             
|------------------------ + 2*e    *cos\\/ 6 *t/|*Heaviside(t)
\           6                                   /             

 Y(s)_total = ZIR + ZSR:
   2*s + 5   
-------------
 2           
s  + 4*s + 10

 y(t)_total = ZIR + ZSR:
/  ___  -2*t    /  ___  \                       \             
|\/ 6 *e    *sin\\/ 6 *t/      -2*t    /  ___  \|             
|------------------------ + 2*e    *cos\\/ 6 *t/|*Heaviside(t)
\           6                                   /             
>>>

donde la gráfica de polos muestra que se encuentran todos del lado izquierdo del plano

OpAmpEj01 H(s) Polos

También se muestra la respuesta al impulso h(t) del circuito

Las respuestas fueron obtenidas al usar como bloque de entrada,

# H(s) respuesta impulso
Ps = 2*s+5
Qs = s**2 + 4*s + 10
Hs = Ps/Qs

# X(s) Señal de entrada
xt = d

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [0, 0] # estado cero no se usan

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras = 101  # 51 resolucion grafica

s1Eva2012TI_T1 LTI CT Polos y ceros de H(s) con respuesta DC conocida

Ejercicio: 1Eva2012TI_T1 LTI CT Polos y ceros de H(s) con respuesta DC conocida

Se revisan los polos y ceros para escribir la función de transferencia,

H(s) = \frac{(s-ceros_0)(s-ceros_1)}{(s-polo_0)(s-polo_1)} H(s) = \frac{(s-(-1.5j))(s-(1.5j))}{(s-(-1+0.5j))(s-(-1-0.5j))}

H(s) = \frac{s^2+2.25}{s^2+2s+1.25}

tomando el modelo proporcionado de la ecuación, se compara para encontrar los valores de las variables

H(s) = \frac{k (s^2 + b_1 s + b_2)}{s^2 + a_1 s + a_2}

k = 1, b1=0, b2=2.25, a1=2, a2=1.25

Separando en fracciones parciales, dado que el grado M del numerador y grado N denominador son iguales, ganancia el coeficiente del término de mayor grado del numerador. Lo mismo que se obtiene al dividir el numerador y denominador por s2 y hacer s→∞

H(s) = \frac{1+\frac{2.25}{s^2}}{1+2\frac{1}{s}+\frac{1.25}{s^2}} \Big|_{s=\infty} =1 H(s) = 1+\frac{k_1s+k_2}{s^2+2s+1.25} H(s) = \frac{s^2+2s+1.25+k_1s+k_2}{s^2+2s+1.25} H(s) = \frac{s^2+(2+k_1)s+(1.25+k_2)}{s^2+2s+1.25}

con lo que 2+k1=0, entonces k1=-2,
también 1.25+k2=2.25, entonces K2 = 1

con lo que se obtiene,

H(s) = 1+\frac{-2s +1}{s^2 + 2 s + 1.25}

los polos se encuentran sobre el lado izquierdo del plano s, por lo que el sistema es asintoticamente estable

Observe que para H(s) el grado del polinomio el numerador P de H(s) es igual al grado del denominador Q. por lo que M=N. por lo que se considera un función impropia (Lathi 4.3b p339, 4.2 p359). Se genera una respuesta impulso por el valor de la constante 1 que se obtiene al usar fracciones parciales.

Para la transformada inversa, para el término cuadrático, en la tabla de transformadas de Laplace se identifica el caso 12c:

\frac{As+B}{s^2+2as+c} donde A=-2, B=1,a=1,c=1.25

El valor b = \sqrt{c-a^2} se calcula como:

b = \sqrt{1.25-(1)^2} = \frac{1}{2} r = \sqrt{\frac{A^2 c +B^2 -2ABa}{c-a^2}} r = \sqrt{\frac{(-2)^2 (1.25) +(1)^2 -2(-2)(1)(1)}{1.25-(-2)^2}} =\sqrt{\frac{10}{1/4}} = 6.3245 \theta = \tan ^{-1} \Big( \frac{Aa-B}{A\sqrt{c-a^2}}\Big) \theta = \tan ^{-1} \Big( \frac{(-2)(1)-(1)}{(-2)\sqrt{1.25-(1)^2}}\Big)\tan ^{-1}(\frac{-3}{-2(1/2)} = 1.2490

h(t) = δ(t) + re-atcos (bt+θ) μ(t)

h(t) = δ(t) + 6.3245e-tcos (t/2+1.2490) μ(t)

Algoritmo en Python

Para el bloque de ingreso se debe mantener los coeficientes como números enteros o racionales. Recordando que los algoritmos de Sympy para transformadas de Laplace (hasta versión 1.11) se encuentran escritos para el dominio de los enteros y racionales (‘ZZ’), por lo que los datos se ingresan como:

k1 = sym.Rational(2.25).limit_denominator(100)
k2 = sym.Rational(1.25).limit_denominator(100)
Hs = (s**2+k1)/(s**2+2*s+k2)

con el resultado mostrado,

 H(s) = P(s)/Q(s)
    2   9   
   s  + -   
        4   
------------
 2         5
s  + 2*s + -
           4

 H(s) fracciones parciales:
   4*(2*s - 1)      
- -------------- + 1
     2              
  4*s  + 8*s + 5    

 h(t) :
                   -t    /t\                   -t    /t\             
DiracDelta(t) + 6*e  *sin|-|*Heaviside(t) - 2*e  *cos|-|*Heaviside(t)
                         \2/                         \2/             

 {polos,veces}:  {-1 - I/2: 1, -1 + I/2: 1}
 polos reales:  0  complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  0  unicos: 0
 asintoticamente:  estable

 parametros Q cuadraticos s: 
{'A': -2.0, 'B': 1.0, 'a': 1.0, 'c': 1.25,
 'r': 6.324555320336759, 'b': 0.5,
 'theta': 1.2490457723982544}

Con gráficas de H(s) y polos

s1Eva2012TI_T1 polos H(s)

gráfica de h(t)

s1Eva2012TI_T1 h(t)

Instrucciones en Python

Usando los bloques desarrollados en la Unidad 4 Sistemas LTI – Laplace  y las funciones resumidas como telg1001.py que pueden ser usados en cada pregunta.

# Y(s) Respuesta total con entada cero y estado cero
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import matplotlib.pyplot as plt
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)
d = sym.DiracDelta(t)
u = sym.Heaviside(t)

# H(s) y estabilidad
k1 = sym.Rational(2.25).limit_denominator(100)
k2 = sym.Rational(1.25).limit_denominator(100)
Hs = (s**2+k1)/(s**2+2*s+k2)
#Hs = (s**2+2.25)/(s**2+2*s+1.25)

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [0, 0] # estado cero no se usan

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 5
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs = fcnm.apart_exp(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_parametros(Hs_fc)

polosceros = fcnm.busca_polosceros(Hs)
Q_polos = polosceros['Q_polos']
P_ceros = polosceros['P_ceros']

estable = fcnm.estabilidad_asintotica(Q_polos)

# H(t) respuesta al impulso
ht = 0*s
term_suma = sym.Add.make_args(Hs)
for term_k in term_suma:
    ht_k = sym.inverse_laplace_transform(term_k,s,t)
    # simplifica log(exp()) ej: e**(-2s)/(s**2)
    if ht_k.has(sym.log):
        ht_k = sym.simplify(ht_k,inverse=True)
    ht  = ht + ht_k
lista_escalon = ht.atoms(sym.Heaviside)
ht = sym.expand(ht,t) # terminos suma
ht = sym.collect(ht,lista_escalon)

# SALIDA
print(' H(s) = P(s)/Q(s):')
sym.pprint(Hs)
print(' H(s) en factores:')
sym.pprint(Hs_fc)
if len(Hs_Qs2)>0:
    print('\nH(s) parámetros cuadraticos:')
    fcnm.print_resultado_dict(Hs_Qs2)

print('\n h(t) :')
sym.pprint(ht)

print('\npolosceros:')
fcnm.print_resultado_dict(polosceros)

print('\nEstabilidad de H(s):')
for k in estable:
    print('',k,':',estable[k])

# Graficas polos, H(s), con polos h(t) --------
muestras_H = 101
figura_s  = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,muestras=muestras_H,f_nombre='H')
figura_ht = fcnm.graficar_ft(ht,t_a,t_b,muestras,f_nombre='h')
plt.show()

.

s2Eva2012TI_T3 LTI DT causal, coeficientes de respuesta impulso h[n]

Ejercicio: 2Eva2012TI_T3 LTI DT causal, coeficientes de respuesta impulso h[n]

La ecuación de diferencias descrita se usa para crear el diagrama de bloques.

y[n] - \frac{5}{4} y[n-1] + \frac{1}{36} y[n-2] + \frac{1}{18} y[n-3] = x[n] - \frac{1}{2} x[n-1]

literal a. Respuesta impulso h[n]

La forma para usar transformadas z requiere que se desplar tres unidades.

y[n+3] - \frac{5}{4} y[n+2] + \frac{1}{36} y[n+1] + \frac{1}{18} y[n] = x[n+3] - \frac{1}{2} x[n+2] z^3 Y[z] - \frac{5}{4} z^2 Y[z] + \frac{1}{36} z Y[z] + \frac{1}{18} Y[z] = z^3X[z] - \frac{1}{2} z^2 X[z] (z^3 - \frac{5}{4} z^2 + \frac{1}{36} z + \frac{1}{18}) Y[z] = (z^3 - \frac{1}{2} z^2) X[z] H[z] = \frac{Y[z]}{X[z]} = \frac{z^3 - \frac{1}{2} z^2}{z^3 - \frac{5}{4} z^2 + \frac{1}{36} z + \frac{1}{18}}

las raices del denominador obtenidas con el algoritmo son:

 polos: {-0.1871842709362: 1, 0.2500: 1, 1.187184270936: 1}

 polos Re: [1/2 - sqrt(17)/6, 1/4, 1/2 + sqrt(17)/6]
>>> 

Observación: existe una raíz con distancia al origen superior al radio=1, por lo que habría un termino creciente no acotado.

El modelo se puede plantear como:

H[z] = \frac{z^2(z - \frac{1}{2})}{(z+0.1871842709362)(z-0.25)(z- 1.187184270936)}

usando fraciones parciales modificadas:

\frac{H[z]}{z} = \frac{z(z - \frac{1}{2})}{(z+0.1871842709362)(z-0.25)(z- 1.187184270936)}

Se puede plantear un modelo de respuesta por cada raiz como:

\frac{H[z]}{z} = \frac{k_1}{(z-0.1871842709362)}+\frac{k_2}{(z-0.25)} +\frac{k_3}{(z- 1.187184270936)}

usando el método de Heaviside,

k_1= \frac{z(z - \frac{1}{2})}{\cancel{(z+0.1871842709362)}(z-0.25)(z- 1.187184270936) }\Big|_{z=-0.1871842709362} k_1= \frac{-0.1871842709362(-0.1871842709362 - \frac{1}{2})}{(-0.1871842709362-0.25)(-0.1871842709362- 1.187184270936) } = 0.2140 k_2= \frac{z(z - \frac{1}{2})}{(z+0.1871842709362)\cancel{(z-0.25)}(z- 1.187184270936) }\Big|_{z=0.25} k_2= \frac{0.25(0.25 - \frac{1}{2})}{(0.25+0.1871842709362)(0.25- 1.187184270936)} = .0.1525 k_3= \frac{z(z - \frac{1}{2})}{(z+0.1871842709362)(z-0.25)\cancel{(z- 1.187184270936) }}\Big|_{z=1.187184270936} k_3= \frac{1.187184270936(1.187184270936 - \frac{1}{2})}{(1.187184270936+0.1871842709362)(1.187184270936-0.25)} = 0.6333

reemplazando en la expresión planteada,

\frac{H[z]}{z} = \frac{0.2140}{(z-0.1871842709362)}+\frac{0.1525}{(z-0.25)} +\frac{0.6333}{(z- 1.187184270936)}

y restaurando en fracciones parciales al multiplicar por z cada lado

H[z] = \frac{0.2140z}{(z-0.1871842709362)}+\frac{0.1525z}{(z-0.25)} +\frac{0.6333z}{(z- 1.187184270936)}

lso valores se ajustan al modelo planteado en el enunciado

h[n] = a \alpha^n \mu [n] + b \beta^n \mu[n] + c \gamma^n \mu [n] h[n] = 0.2140(0.1871842709362)^n \mu [n] + 0.1525 (0.25)^n \mu[n] + 0.6333 (1.187184270936)^n \mu [n]

obtenga entonces los valores pertinentes.

a= 0.2140 b= 0.1525 c=0.6333
α= 0.1871842709362 β=0.25 γ=1.187184270936

 

 

 

s2Eva2011TI_T2 LTI DT Determinar H[z] desde bloques

Ejercicio: 2Eva2011TI_T2 LTI DT Determinar H[z] desde bloques

literal a. expresión de la función de transferencia

El diagrama de bloques del enunciado se reordena de la siguiente forma:

El nuevo diagrama muesta que el sistema tiene dos subcomponentes en paralelo.

H[z] = -\frac{(11/2)z +7}{z^2-1z-2} +\frac{z}{2} -\frac{9}{2}

como existen varios componentes, se pueden tratar de forma separada.

H[z] = H_1[z] +\frac{z}{2} -\frac{9}{2} H_1[z] = -\frac{(11/2)z +7}{z^2-z-2}

las raíces del denominador son

polos: {2: 1, -1: 1}
H_1[z] = -\frac{(11/2)z +7}{(z-2)(z+1)}

usando fracciones parciales modificadas

\frac{H_1[z]}{z} = -\frac{(11/2)z +7}{z(z-2)(z+1)} = \frac{k_1}{z}+\frac{k_2}{z-2}+\frac{k_3}{z+1}

usando el algoritmo se tiene:

\frac{H_1[z]}{z} = \frac{3.5}{z}-\frac{3}{z-2}-\frac{0.5}{z+1}

restaurando a fraciones parciales

H_1[z] = 3.5-3\frac{z}{z-2}-0.5\frac{z}{z+1}

la función de tranferencia H[z] completa como componentes en paralelo es:

H[z] = 3.5-3\frac{z}{z-2}-0.5\frac{z}{z+1} +\frac{z}{2} -\frac{9}{2}

usando el algoritmo de la sección X[z] Fracciones parciales modificadas con Python para H1[z]

con entrada:

Pz = -(11/2*z+7)
Qz = z**2-z-2

se obtiene:

 Hz:
-5.5*z - 7
----------
 2        
z  - z - 2

 Hz/z:
 -5.5*z - 7.0 
--------------
  / 2        \
z*\z  - z - 2/

 Hz/z.apart:
   0.5        1.5       3.5
- ----- - ----------- + ---
  z + 1   0.5*z - 1.0    z 

 Hz = (Hz/z)*z
  0.5*z      1.5*z         
- ----- - ----------- + 3.5
  z + 1   0.5*z - 1.0      

 polos:     {2: 1, -1: 1}

 polos Re:  [-1, 2]

>>> 

resultado que se completan con los términos de los otros componentes.

literal b. Ecuación de diferencias

Se realiza la conversión por la suma de cada componente (en paralelo):

Para y1[n]:

Y_1[z] = H_1[z]X[z] = \Big[-\frac{(11/2)z +7}{z^2-z-2}\Big] X[z] Y_1[z] [z^2-z-2]= -(11/2)zX[z] -7X[z] y[n+2] - y[n+1] -2y[n] = -\frac{11}{2} x[n+1] -7x[n]

Para y2[n]:

Y_2[z] = H_2[z]X[z] = \frac{z}{2}X[z] y[n] = \frac{1}{2}x[n+1]

Para y3[n]:

Y_3[z] = H_3[z]X[z] = -\frac{9}{2} X[z] y[n] = -\frac{9}{2} x[n]

se suman las expresiones obtenidas de Y1[z] + Y2[z]+Y3[z]

y[n+2] - y[n+1] -2y[n] +y[n] + y[n]= -\frac{11}{2} x[n+1] -7x[n] +\frac{1}{2}x[n+1]-\frac{9}{2} x[n]

la ecuación de diferencias simplificada es;

y[n+2] - y[n+1] = -5 x[n+1] -\frac{23}{2} x[n]

El sistema global se puede reescribir nuevamente en z como

z^2Y[z] - zY[z] = -5 zX[z] -\frac{23}{2} X[z] z(z-1)Y[z] = (-5 z -\frac{23}{2}) X[z]

se tiene que el sistema tienen polos en 0 y 1, que se encuentran dentro del radio 1 del plano z.

literal d. Respuesta al impulso H[z]

H[z] se obtiene a partir de la última ecuación

H[z] = \frac{Y[z]}{H[z]} = -\frac{5 z +\frac{23}{2}}{z(z-1)}

aplicando fracciones parciales modificadas:

\frac{H[z]}{z} = -\frac{5 z +\frac{23}{2})}{z^2(z-1)} \frac{H[z]}{z} = 16.5 \frac{1}{z} +11.5\frac{1}{z^2} -16.5 \frac{1}{z-1}

restaurando fracciones parciales al multiplicar por z

H[z] = 16.5 +11.5 \frac{1}{z} -16.5 \frac{z}{z-1}

usando la tabla de transformada z se tiene:

h[n] = 16.5\delta [n] +11.5 \mu [n-1] - 16.5 \mu[n]

usando el algoritmo:

 Hz:
-5*z - 11.5
-----------
 z*(z - 1) 

 Hz/z:
-5.0*z - 11.5
-------------
    / 2    \ 
  z*\z  - z/ 

 Hz/z.apart:
   16.5   16.5   11.5
- ----- + ---- + ----
  z - 1    z       2 
                  z  

 Hz = (Hz/z)*z
  16.5*z          11.5
- ------ + 16.5 + ----
  z - 1            z  

 polos:     {1: 1, 0: 1}

 polos Re:  [0, 1]
>>> 

Instrucciones en Python

# Transformada z- Fracciones parciales
# Polos únicos, repetidos y complejos
# Lathi Ejercicio 5.3c pd495
# blog.espol.edu.ec/telg1001
import sympy as sym
import numpy as np

# INGRESO
z = sym.Symbol('z')

Pz = -(5*z+23/2)
Qz = z*(z-1)

# PROCEDIMIENTO
P = Pz.as_poly(z)
Q = Qz.as_poly(z)
Xz = Pz/Qz

# Raices Denominador
Q_raiz = sym.roots(Q)
Q_raizRe = sym.real_roots(Q)
rcompleja = len(Q_raiz)-len(Q_raizRe)

# Raices reales, únicas y repetidas
if (rcompleja<=0 and len(Q_raizRe)>0):
    # fracciones parciales modificadas
    Xzz = (P/Q)/z
    Xzm = Xzz.apart()
    # fracciones parciales restaurada
    terminos = Xzm.args
    Xzp = 0*z
    for untermino in terminos:
        Xzp = Xzp + untermino*z

def parametro_cuadratico(untermino):
    unparametro ={}
    # revisa denominador cuadratico
    [numerador,denominador] = (untermino).as_numer_denom()
    gradoD = 0
    coeficientesD = denominador
    gradoN = 0
    coeficientesN = numerador
    if not(denominador.is_constant()):
        denominador = denominador.as_poly()
        gradoD = denominador.degree()
        coeficientesD = denominador.coeffs()
    if not(numerador.is_constant()):
        numerador = numerador.as_poly()
        gradoN = numerador.degree()
        coeficientesN = numerador.coeffs()
    if gradoD == 2 and gradoN==2:
        a = float(coeficientesD[1])/2
        gamma2 = float(coeficientesD[2])
        gamma = np.sqrt(gamma2)
        A = float(coeficientesN[0])
        B = float(coeficientesN[1])
        rN = (A**2)*gamma2 + B**2 - 2*A*a*B
        rD = gamma2 - a**2
        r = np.sqrt(rN/rD)
        beta = np.arccos(-a/gamma)
        thetaN = A*a-B
        thetaD = A*np.sqrt(gamma2-a**2)
        theta = np.arctan(thetaN/thetaD)
        unparametro = {'r':r,
                       'gamma':gamma,
                       'beta':beta,
                       'theta':theta}
    return (unparametro)

# Terminos cuadraticos
parametros = [] # denominador cuadratico
if (rcompleja>0 and len(Q_raizRe)>0):
    # fracciones parciales modificadas
    Xzz = (P/Q)/z
    Xzm = Xzz.apart()
    # fracciones parciales restaurada
    terminos = Xzm.args
    Xzp = 0*z
    for untermino in terminos:
        Xzp = Xzp + untermino*z
        unparam = parametro_cuadratico(untermino*z)
        if len(unparam)>0:
            parametros.append(unparam)
if (rcompleja>0 and len(Q_raizRe)==0):
    Xzp = P/Q
    parametros.append(parametro_cuadratico(P/Q))

# SALIDA
print('\n Hz:')
sym.pprint(Xz)
if (len(Q_raizRe)>0):
    print('\n Hz/z:')
    sym.pprint(Xzz)
    print('\n Hz/z.apart:')
    sym.pprint(Xzm)
    print('\n Hz = (Hz/z)*z')
    sym.pprint(Xzp)
    print('\n polos:    ', Q_raiz)
print('\n polos Re: ', Q_raizRe)
if len(parametros)>0:
    for unparam in parametros:
        print('parametros cuadraticos: ')
        for cadauno in unparam.keys():
            print(cadauno,'\t',unparam[cadauno])