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 desplazar 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 raíces del denominador obtenidas con Sympy-Python:
>>> sym.roots(Pz)
{1/2: 1, 0: 2}
>>> sym.factor(Pz)
z**2*(2*z - 1)/2
>>> sym.factor(Qz)
(4*z - 1)*(9*z**2 - 9*z - 2)/36
>>> sym.roots(Qz)
{1/4: 1, 1/2 - sqrt(17)/6: 1, 1/2 + sqrt(17)/6: 1}
>>>
Observación: existe una raíz con distancia al origen superior al radio=1, por lo que habría un termino creciente no acotado.
que se observa también como:
El modelo se puede plantear como:
H[z] = \frac{z^2\Big(z-\frac{1}{2}\Big)}{\Big(z-\frac{1}{4}\Big)\Big(z^2-z-\frac{2}{9}\Big)}para usar fracciones parciales modificadas, se multiplica cada lado por 1/z:
\frac{H[z]}{z} = \frac{z\Big(z-\frac{1}{2}\Big)}{\Big(z-\frac{1}{4}\Big)\Big(z^2-z-\frac{2}{9}\Big)}Se puede plantear un modelo de respuesta por cada raíz como:
\frac{H[z]}{z} = \frac{k1}{z-\frac{1}{4}}+\frac{k_2 z +k_3}{z^2-z-\frac{2}{9}}usando el método de Heaviside,
k_1 = \frac{z\Big(z-\frac{1}{2}\Big)}{\cancel{\Big(z-\frac{1}{4}\Big)}\Big(z^2-z-\frac{2}{9}\Big)}\Big|_{z=\frac{1}{4}} k_1 = \frac{\frac{1}{4}\Big(\frac{1}{4}-\frac{1}{2}\Big)}{\Big(\frac{1}{4}^2-\frac{1}{4}-\frac{2}{9}\Big)} k_1 = -\frac{1}{16}\frac{1}{\frac{(4)(9)-(16)(9)-(2)(16)(4)}{(16)(4)(9)}} =\frac{36}{236}con lo que ahora en la expresión se convierte en:
\frac{H[z]}{z} = \frac{z\Big(z-\frac{1}{2}\Big)}{\Big(z-\frac{1}{4}\Big)\Big(z^2-z-\frac{2}{9}\Big)} =\frac{\frac{36}{236}}{z-\frac{1}{4}}+\frac{k_2 z +K_3}{z^2-z-\frac{2}{9}} =Usando el método de los factores cuadráticos, se multiplica ambos lados por z y z→∞
\frac{H[z]}{z} z = \frac{z^2\Big(z-\frac{1}{2}\Big)}{\Big(z-\frac{1}{4}\Big)\Big(z^2-z-\frac{2}{9}\Big)} =\frac{\frac{36}{236}z}{z-\frac{1}{4}}+z\frac{k_2 z +K_3}{z^2-z-\frac{2}{9}} =se divide numerador y denominador del lado izquierdo para 1/z3 y el lado derecho el primer termino 1/z y el segundo termino 1/z2
\frac{\Big(1-\frac{1}{2}\frac{1}{z}\Big)}{\Big(1-\frac{1}{4}\frac{1}{z}\Big)\Big(1-\frac{1}{z}-\frac{2}{9}\frac{1}{z^2}\Big)} =\frac{\frac{36}{236}}{1-\frac{1}{4}\frac{1}{z}}+\frac{k_2 +K_3\frac{1}{z}}{1-\frac{1}{z}-\frac{2}{9}\frac{1}{z}}y cuando z→∞
\frac{\Big(1-\frac{1}{2}(0)\Big)}{\Big(1-\frac{1}{4}(0)\Big)\Big(1-(0)-\frac{2}{9}(0)\Big)} =\frac{\frac{36}{236}}{1-\frac{1}{4}(0)}+\frac{k_2 +K_3(0)}{1-(0)-\frac{2}{9}(0)} 1 =\frac{36}{236}+k_2 k_2 = 1-\frac{36}{236}=\frac{200}{236} = \frac{50}{59}con lo que K2=50/59 , para encontrar K3 se usa un valor conveniente de z=0
H[0] = \frac{(0)\Big((0)-\frac{1}{2}\Big)}{\Big((0)-\frac{1}{4}\Big)\Big((0)^2-(0)-\frac{2}{9}\Big)} =\frac{\frac{36}{236}}{(0)-\frac{1}{4}}+\frac{\frac{50}{59}(0) +K_3}{(0)^2-(0)-\frac{2}{9}} = 0 =\frac{\frac{36}{236}}{-\frac{1}{4}}+\frac{K_3}{-\frac{2}{9}} = K_3 = \frac{2}{9}\frac{\frac{36}{236}}{-\frac{1}{4}} = -4\frac{2}{9}\frac{36}{236} = -\frac{8}{59} \frac{H[z]}{z} = \frac{36}{236}\frac{1}{z-\frac{1}{4}}+\frac{50/59 z -8/59}{z^2-z-\frac{2}{9}} = \frac{H[z]}{z} = \frac{36}{236}\frac{1}{z-\frac{1}{4}}+\frac{2}{59}\frac{25 z -4}{z^2-z-\frac{2}{9}}y restaurando en fracciones parciales al multiplicar por z cada lado
H[z] = \frac{9}{59}\frac{z}{z-\frac{1}{4}}+\frac{2}{59}\frac{z(25 z -4)}{z^2-z-\frac{2}{9}}los valores no se ajustan al modelo planteado en el enunciado, y existe un polo con radio>1, por lo que el sistema es creciente,y otro polo en posición r<0.
h[n] = a \alpha^n \mu [n] + b \beta^n \mu[n] + c \gamma^n \mu [n]obtenga entonces los valores pertinentes.
| a= 9/59 = 0.1525 | b= … creciente sin cota, con polo mayor a 1 | c=… con polo negativo, lo que es un término con signo alternado. |
| α= 1/4 = 0.25 | β = 1/2 + np.sqrt(17)/6 = 1.1871842709362768 | γ = 1/2 – np.sqrt(17)/6 = -0.18718427093627676 |
Resultados iniciales con el algoritmo
Se encuentra que hay un término que no se puede usar para encontrar la transformada z.
Hz:
2
3 z
z - --
2
-------------------
2
3 5*z z 1
z - ---- + -- + --
4 36 18
Hz en fracciones parciales
50*z*(z - 4/25) 9*z
--------------- + ------------
/ 2 2\ 59*(z - 1/4)
59*|z - z - -|
\ 9/
Hz en factores
2
z *(z - 0.5)
---------------------------------------
/ 2 \
(z - 0.25)*\z - z - 0.222222222222222/
{Q_polos:veces}: {1/4: 1, 1/2 - sqrt(17)/6: 1, 1/2 + sqrt(17)/6: 1}
{P_ceros:veces}: {1/2: 1, 0: 2}
parametros cuadraticos:
termino: 50*z*(z - 4/25)/(59*(z**2 - z - 2/9))
r : 120.72788283210394
gamma : None
beta : None
theta : None
estabilidad asintótica en z:
circ1_dentro : 2
circ1_repetidos : 0
circ1_sobre : 0
circ1_fuera : 1
unicos : 3
repetidos : 0
asintota : inestable
h[n]:
-n
9*4 *Heaviside(n)
------------------
59
revisar terminos sin transformada de tabla:
50*z*(z - 4/25)/(59*(z**2 - z - 2/9))
señal discreta h[n]
n : [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
h[n]: [1.52542373e-01 3.81355932e-02 9.53389831e-03 2.38347458e-03
5.95868644e-04 1.48967161e-04 3.72417903e-05 9.31044756e-06
2.32761189e-06 5.81902973e-07]
Instrucciones en Python
Desarrollo con las expresiones iniciales, sin corrección sobre los parámetros que se indican en el enunciado.
Nota: cuando se produzca el siguiente error con Numpy para evaluar una expresión con exponente negativo,
Traceback (most recent call last):
File "D:\MATG1052Ejemplos\Transformadaz\ejercicio03.py", line 93, in
fi = f_n(ki)
File "", line 2, in _lambdifygenerated
return (9/59)*4**(-n)*Heaviside(n, 1/2)
ValueError: Integers to negative integer powers are not allowed.
proceda actualizando los valores a evaluar como tipo real (dtype float), tan solo usando en la línea de ki con lo siguiente:
ki = np.arange(0,muestras_fn,1.0)
quedando las instrucciones de la siguiente forma, que si evalúa valores para realizar gráficas.
# Transformada z- Fracciones parciales # http://blog.espol.edu.ec/telg1001/lti-dt-transformada-z-xz-fracciones-parciales-con-python/ import numpy as np import sympy as sym import matplotlib.pyplot as plt import telg1001 as fcnm #sym.SYMPY_DEBUG=True # INGRESO z = sym.Symbol('z') n = sym.Symbol('n', real=True) # coeficientes como racional en dominio 'ZZ' enteros a0 = sym.Rational(1,2) a1 = sym.Rational(5,4) a2 = sym.Rational(1,36) a3 = sym.Rational(1,18) Pz = z**3-a0*z**2 Qz = z**3-a1*z**2+a2*z+a3 F = Pz/Qz # para graficar f_nombre = 'H' # nombre de función[z]: H,X,Y, etc muestras_fn = 10 # muestras para f[n] # PROCEDIMIENTO Fz = fcnm.apart_z(F) Fz_factor = sym.factor(F.evalf()) Fz_factor = fcnm._round_float_is_int(Fz_factor) # polos y ceros de Hz [P,Q] = F.as_numer_denom() P = sym.poly(P,z) Q = sym.poly(Q,z) P_ceros = sym.roots(P) Q_polos = sym.roots(Q) estable_z = fcnm.estabilidad_asintotica_z(Q_polos) # Inversa de transformada z fn = 0*n ; f_noeval = 0*n ; Qz2_term =[] term_sum = sym.Add.make_args(Fz) for term_k in term_sum: term_kn = fcnm.inverse_z_transform(term_k,z,n) if type(term_kn)==tuple: fn = fn + term_kn[0] elif term_kn is not None: fn = fn + term_kn elif term_kn is None: f_noeval = f_noeval + term_k Qz2 = fcnm.Q_cuad_z_parametros(term_k) if Qz2: Qz2_term.append(Qz2) fn = fn.collect(sym.Heaviside(n)) fn = fn.collect(sym.DiracDelta(n)) # SALIDA print('\n '+f_nombre+'z:') sym.pprint(F) print('\n '+f_nombre+'z en fracciones parciales') sym.pprint(Fz) print('\n '+f_nombre+'z en factores') sym.pprint(Fz_factor) print('\n {Q_polos:veces}:',Q_polos) print(' {P_ceros:veces}:',P_ceros) if len(Qz2_term)>0: print('\nparametros cuadraticos: ') for i in range(0,len(Qz2_term),1): for unterm in Qz2_term[i]: print(' termino:',unterm) fcnm.print_resultado_dict(Qz2_term[i][unterm]) print('\nestabilidad asintótica en z:') fcnm.print_resultado_dict(estable_z) print('\n '+f_nombre.lower()+'[n]:') sym.pprint(fn) if not f_noeval==sym.S.Zero: print('revisar terminos sin transformada de tabla:') print(f_noeval) # # GRAFICA ----------- fig_ROC = fcnm.graficar_Fz_polos(Fz_factor,Q_polos,P_ceros, muestras=101,f_nombre=f_nombre) fig_Fz = fcnm.graficar_Fs(Fz_factor,Q_polos,P_ceros, muestras=101, f_nombre=f_nombre) # graficar f[n] ------- f_n = sym.lambdify(n,fn.expand(),modules=fcnm.equivalentes) ki = np.arange(0,muestras_fn,1.0) fi = f_n(ki) print('\nseñal discreta '+f_nombre.lower()+'[n]') print('n :',ki) print(f_nombre.lower()+'[n]:',fi) # graficar f[n] fig_fn, grafxn = plt.subplots() plt.axvline(0,color='grey') plt.stem(ki,fi) plt.grid() plt.xlabel('n') plt.ylabel(f_nombre.lower()+'[n]') etiqueta = r''+f_nombre.lower()+'[n]= $'+str(sym.latex(fn))+'$' plt.title(etiqueta) plt.show()