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.
# TELG1001 - Señales y Sistemas ver 2017/05/29# Resumen de problemas resueltos en clases /ESPOL/FCNM-FIEC/# https://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# https://blog.espol.edu.ec/telg1001/transformada-de-laplace-para-ft-con-sympy-python/defsepara_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:
ifnot(term_mul.has(t)):
constante = term_mul
termino = termino/constante
return([termino,constante])
...
Para interpretar mejor la operación se presenta una animación gráfica de h(t-τ) y y(t), para diferentes valores de t en el intervalo de observación [t_a,t_b] se incorpora la función graf_animada_xh_y().
Otros ejemplos de gif animado de gráficas con matplotlib, se pueden revisar en: Una partícula en movimiento circular del curso CCPG1001 Fundamentos de Programación.
# GRAFICA CON ANIMACION ------------defgraf_animada_xh_y(xt,ht,yt,t_a,t_b,
muestras=101,y_nombre='y',
reprod_x = 4,retardo = 200,
archivo_nombre = ''):
'''grafica animada convolucionx(t) y h(t)
en dos subgráficas con Parametros de animación trama/foto
y_nombre = 'ZSR' # o y nombre de resultado convolución
reprod_x = 4 # velocidad de reproducción
retardo = 200 # milisegundos entre tramas
archivo_nombre = '' # crea gif animado si hay nombre
'''# grafica evaluación numerica
x_t = sym.lambdify(t,xt,modules=equivalentes)
h_t = sym.lambdify(t,ht,modules=equivalentes)
y_t = sym.lambdify(t,yt,modules=equivalentes)
ti = np.linspace(t_a,t_b,muestras)
xi = x_t(ti)
hi = h_t(ti)
yi = y_t(ti)
import matplotlib.animation as animation
# h(t-tau) para cada t
ht_tau = []
for tau inrange(0,muestras,reprod_x):
ht_tau.append(h_t(ti[tau]-ti))
tramas = len(ht_tau) # tramas creadas# figura con dos sub-graficas
fig_anim = plt.figure()
graf_a1 = fig_anim.add_subplot(211)
graf_a2 = fig_anim.add_subplot(212)
# grafico superior
x_linea, = graf_a1.plot(ti,xi, color='blue',
label=r'$x(\tau)$')
h_linea, = graf_a1.plot(ti,hi,color='magenta',
linestyle='dashed',
label=r'$h(\tau)$')
htau_linea, = graf_a1.plot(ti,ht_tau[0],
color='magenta',
label=r'$h(t-\tau)$')
punto1, = graf_a1.plot(0,0, color='magenta',marker=6)
# grafico inferior
color_y = 'green'if y_nombre=='ZSR':
color_y ='dodgerblue'
y_linea, = graf_a2.plot(ti,yi, color=color_y,
label=y_nombre+'(t)')
punto2, = graf_a2.plot(0,0, color=color_y,marker=6)
y_sombra, = graf_a2.plot(ti,yi, color=color_y)
y_sombra.set_visible(False) # Para fijar leyend()# Configura gráfica
titulo = r''+y_nombre+'(t)= x(t)$\circledast$h(t)'
graf_a1.set_title(titulo)
ymax1 = np.max([np.max(xi),np.max(hi)])*1.11
ymin1 = np.min([np.min(xi),np.min(hi)])-0.1*ymax1
graf_a1.set_xlim([t_a,t_b])
graf_a1.set_ylim([ymin1,ymax1])
graf_a1.set_xlabel(r'$\tau$')
graf_a1.legend()
graf_a1.grid()
ymax2 = np.max(yi)*1.1
ymin2 = np.min(yi)-0.1*ymax2
graf_a2.set_xlim([t_a,t_b])
graf_a2.set_ylim([ymin2,ymax2])
graf_a2.set_xlabel('t')
graf_a2.legend()
graf_a2.grid()
# cuadros de texto en gráfico
txt_x = (t_b+t_a)/2
txt_y = ymax1*(1-0.09)
txt_tau = graf_a1.text(txt_x,txt_y,'t='+str(t_a),
horizontalalignment='center')
deftrama_actualizar(i,ti,ht_tau):
# actualiza cada linea
htau_linea.set_xdata(ti)
htau_linea.set_ydata(ht_tau[i])
hasta = i*reprod_x
porusar = (muestras-reprod_x*(i+1))
if porusar>=reprod_x: # en intervalo
y_linea.set_xdata(ti[0:hasta])
y_linea.set_ydata(yi[0:hasta])
punto1.set_xdata([ti[hasta]])
punto1.set_ydata([0])
punto2.set_xdata([ti[hasta]])
punto2.set_ydata([0])
else: # insuficientes en intervalo
y_linea.set_xdata(ti)
y_linea.set_ydata(yi)
punto1.set_xdata([ti[-1]])
punto1.set_ydata([0])
punto2.set_xdata([ti[-1]])
punto2.set_ydata([0])
# actualiza texto
t_trama = np.around(ti[i*reprod_x],4)
txt_tau.set_text('t= '+str(t_trama))
return(htau_linea,y_linea,punto1,punto2,txt_tau)
deftrama_limpiar(): # Limpia Trama anterior
htau_linea.set_ydata(np.ma.array(ti, mask=True))
y_linea.set_ydata(np.ma.array(ti, mask=True))
punto1.set_ydata(np.ma.array(ti, mask=True))
punto2.set_ydata(np.ma.array(ti, mask=True))
txt_tau.set_text('')
return(htau_linea,y_linea,punto1,punto2,txt_tau)
i = np.arange(0,tramas,1) # Trama contador
ani = animation.FuncAnimation(fig_anim,trama_actualizar,i ,
fargs = (ti,ht_tau),
init_func = trama_limpiar,
interval = retardo,
blit=True)
# Guarda archivo GIF animado o videoiflen(archivo_nombre)>0:
ani.save(archivo_nombre+'_animado.gif',
writer='imagemagick')
#ani.save(archivo_nombre+'_video.mp4')
plt.draw()
#plt.show()return(ani)
# grafica animada de convolución
n_archivo = ''# sin crear archivo gif animado n_archivo = 'convolucion01' # requiere 'imagemagick'
figura_animada = graf_animada_xh_y(x,h,y,t_a,t_b,
muestras, reprod_x = 4,
archivo_nombre = n_archivo)
plt.show()
Para obtener un gif animado se debe asignar un 'archivo_nombre' para identificar el ejercicio; si se mantiene el nombre vacio '', solamente se crea la gráfica. El resultado se almacena en el archivo_nombre_animado.gif del directorio de trabajo,
La función también se incorpora a telg1001.py para su uso posterior en los ejercicios
La multiplicación entre impulso δ(t) o escalón unitario μ(t) aparece al desarrollar en integral de convolución cuando ambas señales son de tipo causal.
Las expresiones resultantes pueden incluir algunas de éstas operaciones como se muestra en el ejemplo 3 de una convolución entre señales tipo rectangular y rampa de carácter causal.
Con respuestas al usar el algoritmo con Sympy que integran varias multiplicaciones de escalón unitario.
Como Sympy '1,11,1' revisada en Enero del 2023, aún no incluye este tipo de operaciones, se realizan dos funciones: simplifica_impulso() y simplifica_escalon(). Se desarrolla como un ejercicio de análisis de expresiones con Sympy.
>>> sym.__version__
'1.11.1'
Simplificar multiplicación de impulso unitario δ(t)
Las operaciones encontradas en los ejercicios del curso dan una idea básica de por dónde empezar. También se probó las instrucciones doit(), evalf() sin cambios.
Como punto de partida de plantea encontrar todos los puntos de la expresión donde intervienen los impulsos. La función busca_impulso() permite revisar si existiendo impulsos en diferentes tiempos, el resultado deberá ser cero.
x = d
h = d.subs(t,t-1)
xh = x*h
Para facilitar el análisis se realizan las gráficas de las dos funcines que se multiplican y la operación resultante del algoritmo:
otro caso a considerar es:
x = sym.sin(t)
h = d.subs(t,t-2)
xh = x*h
Instrucciones con Python
El algoritmo se desarrolla como funciones, para ser incorporadas a telg1001.py
# Simplificar multiplicacion de impulso unitario d(t)*d(t-1)# https://blog.espol.edu.ec/telg1001/simplificar-multiplicacion-impulso-o-escalon-unitario-sympy/import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t), respuesta impulso h(t)# x = d# h = d.subs(t,t-1)# h = d # h= d**2# x = 2# h = d# h = sym.pi*d.subs(t,t-2)
x = sym.sin(t)
h = d.subs(t,t-2)
# h = 5*d.subs(t,t-2)**2# x = d.subs(t,t-1)# h = u.subs(t,t-1)# h = d*u.subs(t,t-2)+ u*u.subs(t,t-2)# h = d*d.subs(t,t-1)*u.subs(t,t+2)+3# grafica intervalo [t_a,t_b] plano simétrico
t_b = 4 ; t_a = -t_b
muestras = 81
# PROCEDIMIENTOdefsimplifica_impulso(ft):
''' simplifica d**2, d(t-1)*d
'''defsimplifica_d_d(ft):
'''un termino de suma d**2, d(t+1)*d,
'''
respuesta = ft
if ft.has(sym.DiracDelta):# tiene impulsos
impulso_en = fcnm.busca_impulso(ft)
iflen(impulso_en)==0: # anulado por d(t-a)*d(t-b)
respuesta = 0*t
eliflen(impulso_en)>0: # tiene impulsos
respuesta = 1
factor_mul = sym.Mul.make_args(ft)
for factor_k in factor_mul:
ifnot(factor_k.has(sym.DiracDelta)):
ifnot(factor_k.has(sym.Heaviside)):
termino = factor_k.subs(t,impulso_en[0])
else: # tiene escalón
termino = factor_k.subs(t,impulso_en[0])
if termino == 1/2: #coinciden d,u
termino = 1
respuesta = respuesta*termino
else: # factor con impulsoif factor_k.is_Pow: # tiene exponente
respuesta = respuesta*factor_k.args[0]
else: # termino sin exponente
respuesta = respuesta*factor_k
return(respuesta)
# revisa terminos suma
respuesta = 0*t
ft = sym.expand(ft,t)
term_suma = sym.Add.make_args(ft)
for term_k in term_suma:
respuesta = respuesta + simplifica_d_d(term_k)
return(respuesta)
xh = x*h
ft = simplifica_impulso(xh)
# SALIDAprint('xh inicial:')
sym.pprint(xh)
print('h simplificado:')
sym.pprint(ft)
# grafica
figura = fcnm.graficar_xh_y(x,h,ft,t_a,t_b,muestras,y_nombre='xh')
plt.show()
Simplificar multiplicación de escalón unitario μ(t)
Las operaciones encontradas en los ejercicios del curso dan una idea básica de por dónde empezar
Como punto de partida de plantea encontrar todos los puntos de la expresión donde intervienen el escalón unitario. La función busca_escalon() permite revisar la ubicación y el sentido de desarrollo de cada escalón unitario, con estos datos se puede determinar si la función se anula con la otra o la parte que permanece.
x = u
h = u.subs(t,t-1)
xh = x*h
Con una grafica se puede comprobar el escalón que se superpone en la multiplicación de x(t)h(t).
otro ejercicio a considerar es cuando los escalones unitarios se pueden complementar o anular. En el siguiente ejemplo x(t) h(t) se anulan entre si.
x = u.subs(t,t-1)
h = u.subs(t,-t-1)
xh = x*h
mientras en en otro caso, las señales tienen una región donde pueden coexistir
x = u.subs(t,-t+1)
h = u.subs(t,t-0)
Instrucciones con Python
# Simplificar multiplicacion de escalon unitario u(t)*u(t-1)# https://blog.espol.edu.ec/telg1001/simplificar-multiplicacion-impulso-o-escalon-unitario-sympy/import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t), respuesta impulso h(t)# x = u# h = u.subs(t,t-1)# x = u.subs(t,t-1)# h = u.subs(t,-t-1)
x = u.subs(t,-t+1)
h = u.subs(t,t-0)
#x = 3*u.subs(t,-t+2)#h = 2*u.subs(t,-t+3)# x = u.subs(t,t+1)*u.subs(t,-t+2)# h = 1# h = d.subs(t,-t+2)#*u.subs(t,t-3)# grafica intervalo [t_a,t_b] plano simétrico
t_b = 4 ; t_a = -t_b
muestras = 301
# PROCEDIMIENTOdefsimplifica_escalon(ft):
''' simplifica multiplicaciones
Heaviside(t-a)*Heaviside(t-b) en f(t)
'''defsimplifica_u_u(ft):
'''solo dos pares de u(t-1)*u(t-2),
sin coeficientes
'''
donde_u = fcnm.busca_escalon(ft)
donde_u = np.array(donde_u)
# direccion donde_[:,1],# lugar donde[:,0]# analiza multiplicación
resultado = ft # todo igualif donde_u[0,1]*donde_u[1,1] > 0: # direccion igual
k = 0
if donde_u[0,1]>0: # hacia derecha
k = np.argmax(donde_u[:,0])
k_signo = 1
else: # hacia izquierda
k = np.argmin(donde_u[:,0])
k_signo = -1
ubica = donde_u[k,1]*t-k_signo*donde_u[k][0]
resultado = sym.Heaviside(ubica)
else: # direccion diferenteif donde_u[0][1]>0 and (donde_u[0,0]>donde_u[1,0]):
resultado = 0
if donde_u[0][1]<0 and (donde_u[0,0]<=donde_u[1,0]):
resultado = 0
return(resultado)
defsimplifica_u_term(ft):
''' simplifica un termino de varios
factores que multiplican 2*pi*u*u(t-1)
'''
respuesta = ft
if ft.has(sym.Heaviside): # tiene escalon
escalon_en = fcnm.busca_escalon(ft)
revisa = 1 ; otros = 1 ; cuenta = 0
factor_mul = sym.Mul.make_args(ft)
for factor_k in factor_mul:
if factor_k.has(sym.Heaviside):
if factor_k.is_Pow: # con exponente
revisa = revisa*factor_k.args[0]
cuenta = cuenta + 1
else: # sin exponente
revisa = revisa*factor_k
cuenta = cuenta + 1
if cuenta>1: # simplificar
revisa = simplifica_u_u(revisa)
cuenta = len(fcnm.busca_escalon(revisa))
else: # factor sin Heaviside
otros = otros*factor_k
respuesta = otros*revisa
return(respuesta)
# revisa terminos suma
respuesta = 0*t
ft = sym.expand(ft,t)
if ft.has(sym.DiracDelta): # tiene impulsos
ft = fcnm.simplifica_impulso(ft)
term_suma = sym.Add.make_args(ft)
for term_k in term_suma:
respuesta = respuesta + simplifica_u_term(term_k)
return(respuesta)
xh = x*h
ft = simplifica_escalon(xh)
# SALIDAprint('busca_impulso(ft)', fcnm.busca_impulso(h))
print('busca_escalon(ft)', fcnm.busca_escalon(h))
print('x*h inicial:')
sym.pprint(xh)
print('simplificado:')
sym.pprint(ft)
# grafica
figura = fcnm.graficar_xh_y(x,h,ft,t_a,t_b,y_nombre='xh')
plt.show()
En los sistemas lineales, la causalidad se muestra con respuestas para t≥0, por lo que se usan términos multiplicados por escalón unitario μ(t). Mientras que la respuesta al impulso genera al menos un término con impulso unitario δ(t) con un coeficiente constante.
Para el desarrollo de los ejercicios, es de gran utilidad disponer de una función que permita buscar el impulso unitario dentro de una expresión de varios términos. Listar las ubicaciones de los impulsos es un resultado que se usa en la generación de gráficas, pues no siempre la muestra en el intervalo de tiempo usado se encuentra exactamente sobre la aplicación del impulso unitario.
Se propone usar dos funciones. una para encontrar dónde ocurre el impulso y otra para buscar donde y el sentido del escalón unitario.
1. Busca impulso unitario, DiracDelta(t) o δ(t) con Sympy
Si f(t) fuese solo δ(t), δ(t-1) o δ(t+1), el valor de desplazamiento de obtiene con el primer argumento de la función mediante ft.args[0] que al evaluarlo en cero y multiplicarlo por (-1) se tendría donde se aplica.
ecuacion = sym.Eq(ft.args[0],0)
donde = sym.solve(ecuacion,t)
Para los casos como δ(-t-1), δ(-t+1) es necesario trabajar con el signo de la variable t. Una forma de evaluar la expresión interiores de los impulsos es obtener la ecuacion de la forma at+b=0 con sym.Eq(ft.args[0],0) y luego despejar el valor de t con la instrucción sym.solve(ecuacion,t).
En expresiones simples de un solo término con coeficientes, se tiene:
# Busca impulso unitario en h(t)# https://blog.espol.edu.ec/telg1001/busca-impulso-unitario-o-escalon-unitario-con-sympy/import sympy as sym
import matplotlib.pyplot as plt
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
a = sym.Symbol('a',real=True)
# busca impulso en h(t)# h = d# h = 2*sym.pi*d.subs(t,t-2)
h = d.subs(t,t-1) + 2*sym.pi*d.subs(t,t-2)
# h = 3*d.subs(t,t+1)# h = 5*sym.exp(-t)*d# h = 4# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 51
# PROCEDIMIENTOdefbusca_impulso(ft):
''' busca en f(t) impulsos sym.DiracDelta
entrega una lista ordenada de resultados
'''defimpulso_donde(ft):
''' busca posicion de impulso en ft simple d, d**2
un solo termino,sin factores
'''
ft = sym.sympify(ft,t) # convierte a sympy si es constante
donde = [] # revisar f(t) sin impulsosif ft.has(sym.DiracDelta):
ifnot(ft.is_Pow): # sin exponentes
ecuacion = sym.Eq(ft.args[0],0)
else: # con exponentes d**2
ecuacion = sym.Eq(ft.args[0].args[0],0)
donde = sym.solve(ecuacion,t)
return(donde)
defimpulso_donde_unterm(ft):
''' revisa impulso en un termino suma de ft
'''
donde = [] # revisar f(t) sin impulsos
factor_mul = sym.Mul.make_args(ft)
for factor_k in factor_mul:
if factor_k.has(sym.DiracDelta):
donde_k = impulso_donde(factor_k)
iflen(donde) == 0: # vacio
donde.extend(donde_k)
iflen(donde)>0: # sin repetirifnot(donde_k[0] in donde):
donde = [] # anulado por d(t-a)*d(t-b)return(donde)
# revisa terminos suma
ft = sym.sympify(ft,t) # convierte a sympy si es constante
ft = sym.expand(ft)
respuesta = []
term_suma = sym.Add.make_args(ft)
for term_k in term_suma:
donde = impulso_donde_unterm(term_k)
ifnot(donde in respuesta) andlen(donde)>0:
respuesta.extend(donde)
iflen(respuesta)>1: # ordena ascendente
respuesta.sort()
respuesta = list(respuesta)
return(respuesta)
donde_d = busca_impulso(h)
# SALIDAprint('h(t):')
sym.pprint(h)
print('busca_impulso(h): ',donde_d)
# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,muestras,f_nombre='h')
plt.show()
1. Busca escalón unitario, Heaviside(t) o μ(t) con Sympy
Para pruebas se tienen algunas expresiones de h(t), recordando que se analiza expresiones de un solo término, sin sumas o mutiplicaciones.
# entrada x(t)
# h = u
h = u.subs(t,t-1)
# h = u.subs(t,t+1)
# h = 5
Semejante al ejemplo anterior, se busca dónde se produce la transición del escalón unitario y el sentido si es de subida o de bajada como el signo de 't'.
\frac{d}{dt}\mu (t) = \delta(t)
Con la función busca_impulso() obtiene donde se produce el impulso.
El sentido del escalón se obtiene al derivar la expresión del argumento. Si el escalón se desarrolla hacia de donde se produce, será positivo, y hacia la izquierda será negativo. Al aplicar escalon_donde(ft) permite determinar el punto de subida y bajada.
Para determinar el sentido del escalón unitario, se procede a tomar la expresion del argumento ft.args[0] y luego derivar con la instrucción sym.diff(ft.args[0]) para obtener el signo_t. con esto se observa si el desarrollo es hacia el lado derecho del plano (RHS). Con los parametros de donde y el sentido del escalón se puede aplicar un criterio de causalidad.
defescalon_donde(ft):
''' ft sin factores o terminos suma
'''
ft = sym.sympify(ft,t) # convierte a sympy si es constante
respuesta = []
if ft.has(sym.Heaviside):
eq_u = sym.Eq(ft.args[0],0)
sentido = sym.diff(eq_u.lhs,t,1)
donde = sym.solve(eq_u,t)[0]
respuesta = [donde,sentido]
return(respuesta)
Instrucciones en Python
# Busca escalon unitario en h(t)# https://blog.espol.edu.ec/telg1001/busca-impulso-unitario-o-escalon-unitario-con-sympy/import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
k = sym.Symbol('k',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada h(t)# h = u
h = u.subs(t,t-1)
# h = u.subs(t,t+1)# h = 5# h = u + d# h = 3*sym.pi*u.subs(t,t-1)+ 2*d# h = sym.exp(-2*t)*u.subs(t,t+1) + sym.exp(-4*t)*u.subs(t,t-3)# h = sym.exp(-2*t)*sym.cos(t) + 4# h = u.subs(t,-t+2)*u.subs(t,-t+1)# h = u.subs(t,-t+1)*u.subs(t,-t+1)# h = 3*u.subs(t,-t+1)*u.subs(t,-t+1)# h = u.subs(t,t-1)-u.subs(t,t-2)# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 51
# PROCEDIMIENTOdefbusca_escalon(ft):
''' busca en f(t) el donde,sentido de escalon unitario
en un termino simple sin sumas, para ubicar
lado del plano de f(t)
'''defescalon_donde(ft):
''' ft sin factores o terminos suma
'''
ft = sym.sympify(ft,t) # convierte a sympy si es constante
respuesta = []
if ft.has(sym.Heaviside):
eq_u = sym.Eq(ft.args[0],0)
sentido = sym.diff(eq_u.lhs,t,1)
donde = sym.solve(eq_u,t)[0]
respuesta = [donde,sentido]
return(respuesta)
defescalon_donde_unterm(ft):
'''revisa termino de factores multiplica
'''
respuesta = []
factor_mul = sym.Mul.make_args(ft)
for factor_k in factor_mul:
if factor_k.has(sym.Heaviside):
ubicado = []
ifnot(factor_k.is_Pow): # sin exponente
ubicado = escalon_donde(factor_k)
else: # con exponentes d**k
ubicado = escalon_donde(factor_k.args[0])
iflen(ubicado)>0:
respuesta.append(ubicado)
return(respuesta)
# revisa terminos suma
ft = sym.sympify(ft,t) # convierte a sympy si es constante
ft = sym.expand(ft)
respuesta = []
term_suma = sym.Add.make_args(ft)
for term_k in term_suma:
donde = escalon_donde_unterm(term_k)
ifnot(donde in respuesta) andlen(donde)>0:
respuesta.extend(donde)
iflen(respuesta)>1: # ordena ascendente
respuesta = np.array(respuesta)
columna = respuesta[:,0]
ordena = np.argsort(columna)
respuesta = respuesta[ordena]
return(respuesta)
donde_u = busca_escalon(h)
# SALIDAprint('h(t):')
sym.pprint(h)
print('busca_impulso:',fcnm.busca_impulso(h))
iflen(donde_u)<=1:
print('busca_escalon:',donde_u)
else:
print('busca_escalon:')
print(busca_escalon(h))
# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,f_nombre='h')
plt.show()
Las funciones se las incorpora a los algoritmos telg1001.py
La respuesta total del sistema integra las respuestas obtenidas para entrada cero y estado cero. Se resume el desarrollo de un ejercicio con todos los componentes para el ejemplo 1 Modelo entrada-salida para circuitos RLC realizada con Sympy-Python y las funciones del curso en telg1001.py.
Para la respuesta se considera que la respuesta total del sistema se puede describir en forma gráfica como el resultado de dos componentes:
Se plantea el ejercicio para encontrar la respuesta a entrada cero ZIR y la respuesta a estado cero ZSR
Es de considerar que para las gráficas de las respuestas se considera que el sistema es observable, aplicable, usable desde t=0, dado que la solución para entrada cero no se ha acotado en los resultados. Imagine sistema con un capacitor entre sus componentes que a la salida que se descarga desde un voltaje de 1 voltio desde t=0, esta situación no necesariamente implica que podríamos conocer si viene descargandose desde 1.5, 3, 9, 12, o 1000 voltios. Observación a considerar con en la interpretación de los ejercicios para valores obtenidos antes de t=0 que es el punto inicial de observacion.
# Respuesta total del sistema# y(t) = ZIR(t) + ZSR(t)# https://blog.espol.edu.ec/telg1001/lti-ct-yszirzsr-respuesta-del-sistema-con-sympy-python/# Revisar causalidad de x(t) y h(t)import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
y = sym.Function('y')
x = sym.Function('x')
h = sym.Function('h')
u = sym.Heaviside(t)
# ecuacion: lado izquierdo = lado derecho# Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1,evaluate=False)
ecuacion = sym.Eq(LHS,RHS)
# condiciones iniciales [y'(t0),y(t0)]
t0 = 0
cond_inicio = [-5,0]
# entrada x(t)
x = 10*sym.exp(-3*t)*u
# grafica intervalo [t_a,t_b]
t_a = 0; t_b = 5
muestras = 201
# PROCEDIMIENTO# Respuesta entrada cero ZIR
sol_ZIR = fcnm.respuesta_ZIR(ecuacion,cond_inicio,t0)
ZIR = sol_ZIR['ZIR']
# Respuesta al impulso h(t)
sol_h = fcnm.respuesta_impulso_h(ecuacion)
h = sol_h['h']
# respuesta a estado cero ZSR
sol_ZSR = fcnm.respuesta_ZSR(x,h)
ZSR = sol_ZSR['ZSR']
xh = sol_ZSR['xh']
# respuesta a y(t) = ZIR(t)+ZSR(t)
y_total = ZIR+ZSR
# revisa si grafica ZSRifnot(sol_ZSR['cond_graf']):
print('revisar acortar x(t) o h(t) para BIBO')
# SALIDAprint('\n ZIR(t):')
sym.pprint(ZIR)
print('\n h(t):')
sym.pprint(h)
print('\n ZSR(t):')
sym.pprint(ZSR)
print('xcausal: ',sol_ZSR['xcausal'])
print('hcausal: ',sol_ZSR['hcausal'])
print('limites de integral:', sol_ZSR['[tau_a,tau_b]'])
print('\n y(t) = ZIR(t)+ZSR(t):')
sym.pprint(y_total)
# GRAFICA
figura_ZIR = fcnm.graficar_ft(ZIR,t_a,t_b,
muestras,'ZIR')
figura_h = fcnm.graficar_ft(h,t_a,t_b,
muestras,'h')
# grafica animada de convolución
n_archivo = ''# sin crear archivo gif animado # n_archivo = 'LTIC_YTotal_ZIR_Ej01' # requiere 'imagemagick'if sol_ZSR['cond_graf']:
fig_ZSR = fcnm.graficar_xh_y(x,h,ZSR,t_a,t_b,
muestras,y_nombre='ZSR')
fig_ytotal = fcnm.graficar_xhy(ZIR,ZSR,y_total,t_a,t_b,
muestras,x_nombre='ZIR',
h_nombre='ZSR',y_nombre='y')
figura_animada = fcnm.graf_animada_xh_y(x,h,ZSR,-t_b,t_b,
muestras, reprod_x = 4,y_nombre='ZSR',
archivo_nombre = n_archivo)
plt.show()
El algoritmo en Sympy incorpora los pasos realizados con el desarrollo analítico de solución a ecuaciones diferenciales lineales.
Existen otros métodos como el de Transformada de Laplace que simplifica el paso de realizar el integral de convolución, pues al cambiar al dominio 's' en lugar de tiempo las operaciones son mayoritariamete sobre polinomios. La siguiente unidad desarrolla el tema.
El ejemplo es continuación del presentado para respuesta a entrada cero, que tiene la ecuación:
(D^2 + 3D +2)y(t) = Dx(t)
Al diagrama desarrollado en la sección Respuesta entrada cero – Diagrama Bloques, se añaden los bloques para formar la señal de entrada x(t). La variable tiempo se indica con el Bloque Reloj, el escalón μ(t) se configura para empezar en cero y el valor -3 del exponente se los indica con un bloque de ganancia antes de eu.
Para diferenciación de los bloques por sección, la entrada se usa el color azul, el sistema en color naranja y la salida en color morado.
Para continuar con la simulación, asegúrese que no existan condiciones iniciales en el sistema, es decir los valores de estado inicial para los bloques 1/s se encuentren en cero.
Los osciloscopios de observación, tienen su propia referencia de tiempo con otro reloj. El osciloscopio de la señal de entrada x(t) se muestra primero, y la salida en la gráfica posterior.
La señal de salida se obtiene de forma semejante a la respuesta a impulso, es decir en el punto donde está la primera derivada. Revisar teoría para aplicar la forma: P(D)[y(t)μ(t)], que en éste ejercicio se aplica una D.
El resultado de y(t), solo como referencia, es el que se obtiene en el osciloscopio ubicado a la derecha abajo del diagrama.
Con lo que se muestra que es posible observar resultados usando algunos bloques en un simulador.
2. La integral de convolución con Numpy.convolve()
El algoritmo requiere la expresión de las funciones x(t) y h(t) para determinar y(t). Se realiza el cálculo de la convolución, que para aproximar a la integral requiere multiplicar por dt.
Note que debe mantener el valor dt pequeño, en el ejercicio se ha definido como 0.01 que corresponden a casi 1000 muestras en el intervalo [-5,5].
Para usar la función np.convolve(), debido a la naturaleza numérica del algoritmo, se requieren muestras en tiempos, simétricos alrededor de cero. Al mostrar la gráfica, se puede prescindir de la sección negativa, sin embargo, para recordarlo se ha dejado la parte negativa del intervalo de tiempo.
# Señales contínuas en convolución# Compare con el algoritmo en forma discreta# Lathi ejemplo 2.6 p166import numpy as np
import matplotlib.pyplot as plt
# INGRESO
u = lambda t: np.heaviside(t,1)
# u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)
# tiempo [a,b] simétrico alrededor de 0
b = 5 ; a = -b
dt = 0.01
# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)
# Integral de Convolucion x[t]*h[t]# corrección de magnitud por dt para en integral
yi = np.convolve(xi,hi,'same')*dt
# SALIDA - GRAFICA
plt.suptitle('Integral de Convolución x(t)*h(t)')
plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()
plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()
3. La integral de convolución con algoritmo numérico
Esta sección desarrolla el algoritmo para la integral de convolución, tomando paso a paso las operaciones para obtener la gráfica del resultado.
El algoritmo se presenta con propósitos didácticos, pues parece ser un poco más lento debido a que se realizan las operaciones de forma básica (interpretada) de Python. La función np.convolve() usa elementos compilados a lenguaje de maquina que se ejecutan de forma más eficiente.
Note que la única sección que cambia es la de convolución, por lo que la gráfica de salida es la misma que la sección anterior:
# Señales contínuas en convolución. Algoritmo numérico.# compare con el algoritmo en forma discreta# Lathi ejemplo 2.6 p166import numpy as np
import matplotlib.pyplot as plt
# INGRESO
u = lambda t: np.heaviside(t,1)
# u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)
# tiempo [a,b] simétrico alrededor de 0
b = 5 ; a = -b
dt = 0.01
# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)
# Integral de Convolución x[t]*h[t]
muestras = len(xi)
yi = np.zeros(muestras, dtype=float)
for i inrange(0,muestras):
suma=0
for k inrange(0,muestras):
suma = suma + x(ti[k])*h(ti[i]-ti[k])
yi[i]= suma*dt
# SALIDA - GRAFICA
plt.suptitle('Integral de Convolución x(t)*h(t)')
plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()
plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()
Se desarrolla el integral de convolución con Sympy con los criterios indicados en el desarrollo analítico. Para determinar los límites del integral se requiere revisar la Causalidad de h(t) descrita con la función es_causal(ft). Los pasos para los resultados se completan con la función desarrollada en la sección anterior para Integral de convolución. ..
Ejemplo 1. Convolución con x(t) causal y h(t) causal
El límite inferior del integral se usa como 0-, implica aunque se escriba solo 0 se pretende evitar la dificultad cuando x(t) tiene un impulso en el origen.
Para iniciar el algoritmo se definen las variables t y tau, simplificando la expresión del escalón unitario o Heaviside con u.
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
u = sym.Heaviside(t)
Se definen las señales de entrada x(t) y respuesta al impulso h(t) usando las expresiones en Sympy. Para seleccionar los límites del integral de convolución se indica si x(t) y h(t) son causales.
# entrada x(t), respuesta impulso h(t)
x = 10*sym.exp(-3*t)*u
h = (2*sym.exp(-2*t)-sym.exp(-t))*u
Se realiza la comprobación de causalidad revisando cada término de la función, que contenga un escalón o un impulso, tomando como referencia la forma en que se construye h(t) con el impulso y los modos caraterísticos.
En el procedimiento se definen los límites del integral de convolución, la expresión dentro del integral xh(t) y con el resultado se aplica la instrucción expand() para obtener términos suma mas sencillos.
En el caso que la respuesta al impulso h(t) no sea causal y la entrada x(t) es causal, se intercambian las funciones para mantener la simplicidad del integral con límite superior t en lugar de infinito. Se aplica la propiedad conmutativa de la convolución.
# intercambia si h(t) no es_causal# con x(t) es_causal por propiedad conmutativa
intercambia = Falseif hcausal==Falseand xcausal==True:
temporal = h
h = x
x = temporal
xcausal = False
hcausal = True
intercambia = True# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
Se añaden las instrucciones al algoritmo del ejercicio para la convolución de la sección anterior y convierte a una función respuesta_ZSR(x,h) .
Con los resultados del algoritmo, queda añadir las instrucciones para mostrar las ecuaciones y las gráficas como los siguientes:
El integral de convolución de dos funciones x(t) y h(t) usa como operador el símbolo ⊗ que representa el integral de la ecuación mostrada.
La función h(t) representa la repuesta del sistema, función de transferencia o respuesta a impulso, que por la forma de obtenerla se considera causal por tener componente μ(t) y δ(t) que aseguraría que la función se desarrolle en el lado derecho del plano.
La función x(t) representa la entrada del sistema y el integral de convolución se utiliza para determinar la respuesta de estado cero de un sistema ZSR.
Si las funciones x(t) y el sistema h(t) son causales, la respuesta también será causal. Teniendo como límites del integral τa = 0 y τb = t
El límite inferior del integral se usa como 0-, implica aunque se escriba solo 0 se pretende evitar la dificultad cuando x(t) tiene un impulso en el origen.
Ejemplo 1. Algoritmo para el integral de convolución con x(t) y h(t) causales con Sympy
Realizar la convolución entre la respuesta al impulso h(t) y la entrada x(t) mostradas:
x(t) = e^{-t} \mu(t) h(t) = e^{-2t} \mu(t)
Para el algoritmo, se define las variables a usar, escalón e impulso unitarios, definiendo las funciones x(t) y y(t) como:
# entrada x(t), respuesta impulso h(t)
x = sym.exp(-t)*u
h = sym.exp(-2*t)*u
Se revisa la causalidad de las señales para actualizar los límites del integral de convolución.
xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)
# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
El integral de convolución se aplica a la multiplicación entre x(t) y h(t), realizando el cambio de variable por τ y (t-τ). Para facilitar la operación del integral, se simplifica la expresión xh como términos de suma
# integral de convolucion
xh = x.subs(t,tau)*h.subs(t,t-tau)
xh = sym.expand(xh,t)
y = sym.integrate(xh,(tau,tau_a,tau_b))
y = sym.expand(y,t)
y = y.subs(u**2,u) # Heaviside**2=Heaviside
El resultado del integral se muestra como,
xcausal: True
hcausal: True
limites de integral: [ 0 , t ]
x(t)*h(t-tau):
-2*t tau
e *e *Heaviside(tau)*Heaviside(t - tau)
y(t):
-t -2*t
e *Heaviside(t) - e *Heaviside(t)
Instrucciones en Python
# Integral de convolucion con Sympy# Revisar causalidad de x(t) y h(t)import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t), respuesta impulso h(t)
x = sym.exp(-t)*u
h = sym.exp(-2*t)*u
# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 201
# PROCEDIMIENTO
xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)
# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
# integral de convolucion x(t)*h(t)
xh = x.subs(t,tau)*h.subs(t,t-tau)
xh = sym.expand(xh,t)
y = sym.integrate(xh,(tau,tau_a,tau_b))
y = sym.expand(y,t)
ifnot(y.has(sym.Integral)):
y = fcnm.simplifica_escalon(y)
#ZSR = ZSR.subs(u**2,u) # Heaviside**2=Heaviside# SALIDAprint('xcausal: ',xcausal)
print('hcausal: ',hcausal)
print('limites de integral: [',tau_a,',',tau_b,']')
print('x(t)*h(t-tau):')
sym.pprint(xh)
print('\n y(t):')
sym.pprint(y)
if hcausal==False:
print('revisar causalidad de h(t)')
if xcausal==False:
print('revisar causalidad de x(t)')
Grafica para el integral de convolución entre x(t) y h(t)
Con los resultados se puede construir la gráfica, para observar en un intervalo simétrico que permita visualizar la operación de convolución.
Instrucciones complementarias para las gráficas, se realiza como una función a ser añadida a telg1001.py por ser reutilizada en los ejercicios.
# GRAFICA -----defgraficar_xh_y(xt,ht,yt,t_a,t_b,
muestras=101,x_nombre='x',
h_nombre='h',y_nombre='y'):
'''dos subgraficas, x(t) y h(t) en superior
h(t) nn inferior
'''# grafica evaluacion numerica
x_t = sym.lambdify(t,xt,modules=equivalentes)
h_t = sym.lambdify(t,ht,modules=equivalentes)
ti = np.linspace(t_a,t_b,muestras)
xi = x_t(ti)
hi = h_t(ti)
y_t = sym.lambdify(t,yt,modules=equivalentes)
yi = y_t(ti)
colorlinea_y = 'green'if y_nombre =='ZSR':
colorlinea_y = 'dodgerblue'
fig_xh_y, graf2 = plt.subplots(2,1)
untitulo = y_nombre+'(t) = $'+ str(sym.latex(yt))+'$'
graf2[0].set_title(untitulo)
graf2[0].plot(ti,xi, color='blue', label='x(t)')
graf2[0].plot(ti,hi, color='magenta', label='h(t)')
graf2[0].legend()
graf2[0].grid()
graf2[1].plot(ti,yi,colorlinea_y,
label = y_nombre+'(t)')
graf2[1].set_xlabel('t')
graf2[1].legend()
graf2[1].grid()
#plt.show()return(fig_xh_y)
figura_xh_y = graficar_xh_y(x,h,y,t_a,t_b,muestras)
plt.show()
Ejercicio 2. Convolución entre h(t) causal y x(t) no causal
Referencia: Opennheim 2.8 p101
Realizar la convolución entre h(t) y x(t) mostradas:
x(t) = e^{2t} \mu (-t) h(t) = \mu (t-3)
La señal x(t) es creciente desde el lado izquierdo del plano, es decir es no causal, también es acotada en cero por el escalón unitario invertido en el tiempo.
Mientras que h(t) presenta una respuesta con retraso de 3 unidades de tiempo con el escalón unitario, se encuentra en el lado derecho del plano, por lo que se la considera causal.
Para el análisis, el gráfico se presenta la animación de la operación de convolución
Para el ejercicio, es importante considerar que en un sistema lineal causal, solo se desarrolla una salida y(t) ante una entrada x(t). La señal de entrada x(t) no se registra a partir de tiempo con valor cero. Observe en la gráfica que la señal de salida y(t) se desarrolla con retardo de tres unidades respecto a la entrada.
Desde luego el sistema responderá en cualquier momento que se presenta la señal, por lo que x(t) podría ser no causal, La referencia de t=0 se considera desde el punto de vista del sistema h(t).
# entrada x(t), respuesta impulso h(t)
x = sym.exp(-2*t)*u
h = u
# grafica intervalo [t_a,t_b] plano simétrico
t_b = 5 ; t_a = -t_b
muestras = int(t_b)*10+1
Esto implica que de darse el caso que exista un h(t) no causal y x(t) causal, por la propiedad conmutativa de la convolución, para el desarrollo del integral se pueden intercambiar las entradas y hacerlo mas simple.
# intercambia si h(t) no es_causal# con x(t) es_causal por propiedad conmutativa
intercambia = Falseif hcausal==Falseand xcausal==True:
temporal = h
h = x
x = temporal
xcausal = False
hcausal = True
intercambia = True
El resultado con el algoritmo es:
xcausal: False
hcausal: True
limites de integral: [ -oo , t ]
x(t)*h(t-tau):
2*tau
e *Heaviside(-tau)*Heaviside(t - tau - 3)
ZSR(t):
/ -6 2*t \
|e *e 1| 1
|-------- - -|*Heaviside(3 - t) + -
\ 2 2/ 2
revisar causalidad de x(t)
En la respuesta se observa que el témino u(3-t)/2 en -∞ se convierte en 1/2 que anula la constante de 1/2 del tercer término. Al revisar la convergencia del primer término se encuentra que tiende a cero, en consecuencia en la señal hacia el lado izquierdo tiende a cero.
gráfica con el algoritmo
Instrucciones en Python
Para los ejercicios se agrupan las instrucciones en una función respuesta_ZSR(x,h) que recibe las dos señales x(t) y h(t), como se desarrolla en la siguiente sección. En la función se analiza con los algoritmos anteriores si son de tipo causal incluida en telg1001.py.
# Integral de convolucion con Sympy# como Respuesta a estado cero ZSR con x(t) y h(t)# https://blog.espol.edu.ec/telg1001/lti-ct-respuesta-a-estado-cero-con-sympy-python/import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t), respuesta impulso h(t)
x = sym.exp(2*t)*u.subs(t,-t)
h = u.subs(t,t-3)
# grafica intervalo [t_a,t_b] plano simétrico
t_b = 7 ; t_a = -t_b
muestras = int(t_b)*10+1
# PROCEDIMIENTO# revisa causalidad de señales
xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)
# intercambia si h(t) no es_causal# con x(t) es_causal por propiedad conmutativa
intercambia = Falseif hcausal==Falseand xcausal==True:
temporal = h
h = x
x = temporal
xcausal = False
hcausal = True
intercambia = True# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
# integral de convolución x(t)*h(t)
xh = x.subs(t,tau)*h.subs(t,t-tau)
xh = sym.expand(xh,t)
ZSR = sym.integrate(xh,(tau,tau_a,tau_b))
ZSR = sym.expand(ZSR,t)
ifnot(ZSR.has(sym.Integral)):
ZSR = fcnm.simplifica_escalon(ZSR)
#ZSR = ZSR.subs(u**2,u) # Heaviside**2=Heaviside
lista_escalon = ZSR.atoms(sym.Heaviside)
ZSR = sym.expand(ZSR,t) # terminos suma
ZSR = sym.collect(ZSR,lista_escalon)
# restaura si hubo intercambio de x(t) y h(t)if intercambia == True:
xcausal = True
hcausal = False# graficar si no tiene Integral o error
cond_graf = not(ZSR.has(sym.Integral))
cond_graf = cond_graf andnot(ZSR.has(sym.oo))
cond_graf = cond_graf andnot(ZSR.has(sym.nan))
# SALIDAprint('xcausal: ',xcausal)
print('hcausal: ',hcausal)
print('limites de integral: [',tau_a,',',tau_b,']')
print('x(t)*h(t-tau):')
if xh.is_Add:
i=0
for term_suma in xh.args:
if i>0:
print('+')
sym.pprint(term_suma)
i=i+1
else:
sym.pprint(xh)
print('\n ZSR(t):')
sym.pprint(ZSR)
if hcausal==False:
print('revisar causalidad de h(t)')
if xcausal==False:
print('revisar causalidad de x(t)')
if intercambia:
print('considere intercambiar h(t)con x(t)')
ifnot(cond_graf):
print('revisar acortar x(t) o h(t) para BIBO')
# GRAFICAif cond_graf:
figura_xh_y = fcnm.graficar_xh_y(x,h,ZSR,t_a,t_b,
muestras,y_nombre='ZSR')
else: # terminos de integrales sin evaluar, infinito o nan
figura_x = fcnm.graficar_ft(x,t_a,t_b,muestras,'x')
figura_h = fcnm.graficar_ft(h,t_a,t_b,muestras,'h')
#plt.show()# grafica animada de convolución# n_archivo = '' # sin crear archivo gif animado
n_archivo = 'convolucion02'# requiere 'imagemagick'if cond_graf:
figura_animada = fcnm.graf_animada_xh_y(x,h,ZSR,t_a,t_b,
muestras, reprod_x = 4,y_nombre='y',
archivo_nombre = n_archivo)
plt.show()
Para el desarrollo analítico, se propone considerar usar las funciones por intervalos separados:
y(t) = \begin{cases} 0 , & t \lt 0 \\ \frac{1}{2} t^2 , & 0 \lt t \lt T \\ Tt-\frac{1}{2}T^2 , & T \lt t \lt 2T \\ - \frac{1}{2}t^2 + Tt + \frac{3}{2} T^2 , & 2T \lt t \lt 3T \\ 0 . & 3T \lt t \end{cases}
Con las indicaciones se propone como tarea realizar el desarrollo analítico. Observando la gráfica animada y considerando los intervalos propuestos.
El resultado con el algoritmo del ejercicio anterior es:
En el resultado del algoritmo se encuentra que es necesario simplificar la expresión resultante al menos para los términos con escalón unitario, tales como u*u o u*u(t-1). Sympy '1.11.1' aún no incorpora esta simplificación, por lo que se realiza una función que revise término a término la expresión.
