3.3.3 LTI CT – Respuesta a impulso unitario h(t) – Diagrama Bloques

Para usar el simulador con diagrama de bloques para respuesta a impulso unitario, es necesario usar la expresión:

h(t) = b_0 \delta (t) +P(D)[y_n(t) \mu(t)]

Donde si el orden de P(D) es menor al orden de Q(D) entonces b0=0.

Para el ejercicio desarrollado en el sistema modelo, P(D)=D, por lo que hay una sola derivada, que se puede obtener de la entrada del primer integrador de y(t).

( D^2 + 3D +2 ) y(t) = Dx(t)

Para obtener h(t), el impulso unitario x(t), es semejante a tener una condicion inicial sin señal de entrada, siguiendo la teoría descrita de sobre condiciones iniciales.

El impulso unitario representa que y'(t) inicia con 1, por ser el término de más alto orden en P(D) y se configura el diagrama de bloques:

Con lo que la respuesta del sistema al impulso unitario que se obtiene es:

El parámetros inicial se configuró en:

Como referencia para observación se muestra la salida y(t):

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]

 

3.3.2 LTI CT – Respuesta a impulso unitario h(t) con Scipy-Python o Runge-Kutta

2. Respuesta a impulso con sistema LTI definido en Scipy.signal

La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación P, Q y las condiciones de inicio. La respuesta a un impulso unitario δ(t) o función de transferencia h(t) se calcula con una versión muestreada xi a partir del tiempo muestreado en el intervalo [a,b].

sistema = senal.lti(P,Q)
# respuesta a impulso
[t_h, hi] = sistema.impulse(T = ti)

Del problema del «Sistema LTIC – Modelo entrada-salida» y la ecuación de segundo orden, se tiene dos opciones para la respuesta:

– scipy.signal.impulse() donde solo se da el vector ti
– aplicar la teoria de los coeficientes donde solo la derivada de mayor orden de P(D) evaluada en cero, tiene el valor de uno. Todos los demás términos tienen el valor de cero. Revisar la parte conceptual de la sección para respuesta a impulso.

Para el ejercicio y como comprobación se calculan las dos formas, con entrada cero y condición inicial [1,0] y con la función de scipy, que son las usadas en la gráfica.

Instrucciones en Python

# Sistema lineal usando scipy.signal
#   Q(D)y(t)=P(D)x(t)
# ejemplo: (D^2+3D+2)y(t)=(D+0)x(t)
import numpy as np
import scipy.signal as senal
import matplotlib.pyplot as plt

# INGRESO
# Señal de entrada x(t)
x = lambda t: t*0

# Sistema LTI
# coeficients Q  P de la ecuación diferencial
Q = [1., 3., 2.]
P = [1., 0.]
# condiciones iniciales [y'(0),y(0)]
cond_inicio = [1,0]

# grafica, parametros
a = 0 ; b = 5 # intervalo observacion [a,b] 
muestras = 101

# PROCEDIMIENTO
# intervalo observado
ti = np.linspace(a, b, muestras)
xi = x(ti)

sistema = senal.lti(P,Q)

# respuesta entrada cero
[t_y, yi, yc] = senal.lsim(sistema, xi, ti, cond_inicio)
# respuesta a impulso
[t_h, hi] = sistema.impulse(T = ti)

# SALIDA - GRAFICA
plt.subplot(211)
plt.suptitle('Respuesta a Impulso Unitario: h(t) = dy/dt')
plt.plot(ti, xi, color='blue', label='x(t)')
plt.plot(t_h, hi, color='red', label='h(t)')
plt.legend()
plt.grid()

plt.subplot(212)
# plt.plot(t_y, yi, color='magenta', label='y(t)')
# plt.plot(t_y, yc[:,0], 'b-.', label='h(t)')
plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')

plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]
..


3. Algoritmo de Runge-Kutta d2y/dx2 de Análisis numérico

Para el caso del algoritmo de Runge-Kutta, se aplica las condiciones iniciales de y(0)= 0 y y'(0)=1,

siguiendo lo descrito en la parte teórica y se toman los valores de z=y’ para obtener el resultado.

Las instrucciones aplicadas al algoritmo de entrada cero se muestran a continuación:

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y

t0 = 0
y0 = 0
z0 = 1

h = 0.1
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)

# GRAFICA
#plt.plot(ti,yi, label='y(t)')
plt.plot(ti,zi, color = 'Red', label='h(t)=dy/dt')

plt.ylabel('dy/dt')
plt.xlabel('t')
plt.title('Respuesta impulso: Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Recuerde que se aplican las modificaciones de acuerdo al planteamiento del problema, por lo que revise los conceptos antes de aplicar el algoritmo.

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]

3.2.3 LTI CT – Respuesta a entrada cero – Diagrama Bloques

Referencia: Lathi 1.8-1 p111. Oppenheim problema 2.61c p164, Ejemplo 9.24 p700

Un diferenciador D representa un sistema inestable, pues para una entrada acotada como el escalón μ(t) resulta en una salida no acotada como el impulso unitario δ(t). Un sistema con diferenciadores expuesto a ruido amplifica el ruido.

Para los diagramas de bloques y simuladores se prefiere usar integradores (1/D) ≈ (1/s).

Con el problema planteado, en ésta sección se cambia la expresión de la ecuación, ubicando el operador D de mayor grado a la izquierda:

D^2 y(t) + 3Dy(t) + 2y(t) = 0 D^2 y(t) = - 3Dy(t) - 2y(t)

La expresión es semejante a lo realizado para el método de Runge-Kutta presentado en la sección con Python.

Se despeja la respuesta para y(t) del lado izquierdo:

\frac{D^2 y(t)}{D^2} = \frac{-3Dy(t) - 2y(t)}{D^2} y(t) = -3 \frac{1}{D} y(t) - 2 \frac{1}{D^2} y(t)

La forma de la ecuación permite crear el diagrama del sistema, usando 1/D como un integrador y 1/D2 representa una segunda integración en serie.

Los coeficientes de los términos se ubican como ganancias retroalimentadas.

Los valores iniciales se configuran en los integradores observando cual salida se usará en el sumador.

Respuesta entrada cero: [Desarrollo Analítico]  [Sympy-Python]  [Scipy-Python ]  [Runge-Kutta d2y/dx2]  [Simulador]

..


3.1 Usando el Simulador XCOS-Scilab

Para XCOS, se añade un visor de señal de salida (Sinks/Scope) y un reloj de muestreo. Para y'(0)=5 siguiendo el diagrama se ubica y'(t)=dy/dt como la salida del primer integrador . El valor se cambia dando doble click al elemento marcado en rojo y escribiendo el valor de -5.

La gráfica se obtiene usando el osciloscopio (scope)

El tiempo de observación se establece con el valor de «Simulatión/Final Integration time» a 5 segundos.

De forma semejante se puede realizar el ejercicio en otros simuladores, como por ejemplo Simulink-Matlab.


3.2 Desarrollo como circuito en simulador

La simulación en forma de circuito usando las herramientas de Simulink/Matlab y SimPowerSystems, para cada componente (Elements) se crea con «Series RLC Branch».

Para extraer las lecturas de las señal se usa un medidor (Measurements) que se ubica en serie con los demás componentes y se visualiza con un osciloscopio (Simulink/Commonly Used Blocks/Scope).

Las conexiones son similares a las que se usarían en el laboratorio. Se puede observar que las respuestas serán las mismas.

Las condiciones iniciales se pueden establecer en los componentes que pueden tener un valor inicial como el capacitor en 5 voltios.

El signo de la corriente de salida ya está incluido en la inversión de polaridad respecto a la polaridad de la fuente.

Respuesta entrada cero: [Desarrollo Analítico]  [Sympy-Python]  [Scipy-Python ]  [Runge-Kutta d2y/dx2]  [Simulador]

3.2.2 LTI CT – Respuesta a entrada cero con Scipy-Python o Runge-Kutta

Referencia: Lathi 1.8-1 p111. Oppenheim problema 2.61c p164, Ejemplo 9.24 p700

Siguiendo con el ejemplo de la referencia planteado en «Sistema LTIC – Modelo entrada-salida» y las instrucciones del desarrollo analítico para «LTIC – Respuesta entrada cero«, se tiene la expresión:

\frac{d^2}{dt^2}y(t) +3\frac{d}{dt}y(t) + 2y(t) = \frac{d}{dt}x(t)

Se prefiere escribir el termino de mayor orden primero y para entrada cero.Se tiene que x(t)=0, con condiciones iniciales y(0)=0, y'(0)=-5.

(D^2 + 3D +2)y(t) = 0

Respuesta entrada cero: [Desarrollo Analítico]  [Sympy-Python]  [Scipy-Python ]  [Runge-Kutta d2y/dx2]  [Simulador]

..


2. Respuesta entrada cero – desarrollo numérico con Scipy-Python

Esta forma de describir el problema simplifica el desarrollo a partir de la descripción de la ecuación del sistema.

Scipy define un sistema LTI usando los coeficientes de la ecuación ordenados como la expresión P(D)/Q(D con la función Scipy.signal.lti().

Se describe una señal de entrada x(t), muestreada en un intervalo [a,b] junto a las condiciones iniciales en un vector [y'(0),y(0)].

El resultado será la simulación del paso de la señal x(t) por el sistema LTI, usando la función scipy.signal.lsim(), con lo que se obtiene la respuesta del sistema y(t) y los componentes yc=[yi(t),y0(t)], donde y0(t) representa la respuesta a entrada cero.

El algoritmo tiene el siguiente resultado:

Instrucciones en Python

# Sistema lineal usando scipy.signal
#   Q(D)y(t)=P(D)x(t)
# ejemplo: (D^2+3D+2)y(t)=(D+0)x(t)

import numpy as np
import scipy.signal as senal
import matplotlib.pyplot as plt

# INGRESO
# tiempo [a,b] intervalo observación
a = 0 ; b = 5 # intervalo tiempo [a,b]
muestras = 41

# Señal de entrada x(t) = 0
x = lambda t: t*0

# Sistema LTI
# coeficientes Q  P de la ecuación diferencial
Q = [1., 3., 2.]
P = [1., 0.]
# condiciones iniciales [y'(0),y(0)]
cond_inicio = [-5,0]

# PROCEDIMIENTO
# intervalo observado
ti = np.linspace(a, b, muestras)
xi = x(ti)

sistema = senal.lti(P,Q)

# respuesta entrada cero
[t_y, yi, yc] = senal.lsim(sistema, xi, ti, cond_inicio)
# respuesta a impulso
[t_h, hi] = sistema.impulse(T = ti)

# SALIDA - GRAFICA
plt.subplot(211)
plt.suptitle('Respuesta a Entrada cero')
plt.plot(ti, xi, color='blue', label='x(t)')
#plt.plot(t_h, hi, color='red', label='h(t)')
plt.legend()
plt.grid()

plt.subplot(212)
# plt.plot(t_y, yi, color='magenta', label='y(t)')
# plt.plot(t_y, yc[:,0], 'g-.', label='y(t)')
plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

Respuesta entrada cero: [Desarrollo Analítico]  [Sympy-Python]  [Scipy-Python ]  [Runge-Kutta d2y/dx2]  [Simulador]
..


3. Respuesta entrada cero – desarrollo numérico con Runge-Kutta de Métodos numéricos en Python

Se reordena el sistema de ecuaciones para plantear una solución con el algoritmo Runge-Kutta de 2do orden para ecuaciones diferenciales con segundas derivadas

(D^2 + 3D +2)y(t) = 0 D^2 y(t) + 3Dy(t) + 2y(t) = 0

Se ubica el término de mayor grado a la izquierda de la ecuación:

D^2 y(t) = - 3Dy(t) - 2y(t)

sutituyendo las expresiones de las derivadas como las funciones f(x) por z y g(x) por z’:

Dy = y’ = z = f(x)
D2y = y" = z’= -3z – 2y = g(x)

Los valores iniciales de t0=0, y0=0, z0=-5 se usan en el algoritmo.

En este caso también se requiere conocer un intervalo de tiempo a observar [0,tn=5] y definir el tamaño de paso o resolución del resultado h=dt=0.1

El algoritmo permite obtener la gráfica y la tabla de datos.

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y

t0 = 0
y0 = 0
z0 = -5

h = 0.1
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('ti, yi, zi')
print(tabla)

# GRAFICA
plt.plot(ti,yi, color = 'orange', label='y0(t)')

plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('Entrada cero con Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Más detalles del algoritmo de Runge-Kutta para segunda derivada en :

6.2.2 EDO Runge-Kutta d2y/dx2

Respuesta entrada cero: [Desarrollo Analítico]  [Sympy-Python]  [Scipy-Python ]  [Runge-Kutta d2y/dx2]  [Simulador]

3.8 Gráfica animada para interpretar el Integral de convolución con matplotlib-Python

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().

x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

Integral de Convolucion 01 animado

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 ------------
def graf_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 in range(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')

    def trama_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)

    def trama_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 video
    if len(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

3.7 Simplifica multiplicación entre impulso δ(t) o escalón unitario μ(t) con Sympy

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.

x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

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 caracter causal.

x(t) = \mu (t) - \mu (t-1) h(t) = t\mu (t) - t\mu (t-2)

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.

>>> sym.DiracDelta(t)*sym.DiracDelta(t-1)
DiracDelta(t)*DiracDelta(t - 1)
>>> sym.simplify(sym.DiracDelta(t)*sym.DiracDelta(t-1))
DiracDelta(t)*DiracDelta(t - 1)
>>> sym.DiracDelta(t)**2
DiracDelta(t)**2
>>> sym.simplify(sym.DiracDelta(t)**2)
DiracDelta(t)**2
>>>

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:

simplifica impulso Sympy 01

otro caso a considerar es:

x = sym.sin(t)
h = d.subs(t,t-2)

xh = x*h

simplifica impulso Sympy 02

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)
# http://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

# PROCEDIMIENTO
def simplifica_impulso(ft):
    ''' simplifica d**2, d(t-1)*d
    '''
    def simplifica_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)
            if len(impulso_en)==0: # anulado por d(t-a)*d(t-b)
               respuesta = 0*t
            elif len(impulso_en)>0: # tiene impulsos
                respuesta = 1
                factor_mul = sym.Mul.make_args(ft)
                for factor_k in factor_mul:
                    if not(factor_k.has(sym.DiracDelta)):
                        if not(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 impulso
                        if 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)

# SALIDA
print('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

>>> sym.Heaviside(t)*sym.Heaviside(t-1)
Heaviside(t)*Heaviside(t - 1)
>>> sym.simplify(sym.Heaviside(t)*sym.Heaviside(t-1))
Heaviside(t)*Heaviside(t - 1)
>>> sym.simplify(sym.Heaviside(-t+1)*sym.Heaviside(t-1))
Heaviside(1 - t)*Heaviside(t - 1)
>>> sym.simplify(sym.Heaviside(t)**2)
Heaviside(t)**2
>>> sym.simplify(sym.DiracDelta(t-1)*sym.Heaviside(t))
DiracDelta(t - 1)*Heaviside(t)
>>>

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).

simplifica escalon Sympy 01

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

simplifica escalon Sympy 02

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)
# http://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

# PROCEDIMIENTO
def simplifica_escalon(ft):
    ''' simplifica multiplicaciones
        Heaviside(t-a)*Heaviside(t-b) en f(t)
    '''
    def simplifica_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 igual
        if 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 diferente
            if 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)

    def simplifica_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)

# SALIDA
print('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()

..

3.6 Busca impulso unitario δ(t) o escalón unitario μ(t) con Sympy-Python

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:

h = 2*sym.pi*d.subs(t,t-2)

con resultados como:

h(t):
2*pi*DiracDelta(t - 2)
busca_impulso(h):  [2]

busca impulso sympy 01

En el caso que se usen términos de varios impulsos, el valor de interés se encuentra usando el mínimo de la lista donde y observando que aún sea ≥0.

h = d.subs(t,t-1) + 2*sym.pi*d.subs(t,t-2)

con resuldado esperado de algoritmo

h(t):
2*pi*DiracDelta(t - 2) + DiracDelta(t - 1)
busca_impulso(h):  [1, 2]

busca impulso Sympy 02

Instrucciones en Python

# Busca impulso unitario en h(t)
# http://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

# PROCEDIMIENTO
def busca_impulso(ft):
    ''' busca en f(t) impulsos sym.DiracDelta
        entrega una lista ordenada de resultados
    '''
    def impulso_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 impulsos
        if ft.has(sym.DiracDelta):
            if not(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)
    
    def impulso_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)
                if len(donde) == 0: # vacio
                    donde.extend(donde_k)
                if len(donde)>0: # sin repetir
                    if not(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)
        if not(donde in respuesta) and len(donde)>0:
            respuesta.extend(donde)
    if len(respuesta)>1: # ordena ascendente
        respuesta.sort()
    respuesta = list(respuesta)
    return(respuesta)

donde_d = busca_impulso(h)

# SALIDA
print('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.

    def escalon_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)
# http://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

# PROCEDIMIENTO
def busca_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)
    '''
    
    def escalon_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)

    def escalon_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 = []
                if not(factor_k.is_Pow): # sin exponente
                    ubicado = escalon_donde(factor_k)
                else:  # con exponentes d**k
                    ubicado = escalon_donde(factor_k.args[0])
                if len(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)
        if not(donde in respuesta) and len(donde)>0:
            respuesta.extend(donde)
    if len(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)

# SALIDA
print('h(t):')
sym.pprint(h)
print('busca_impulso:',fcnm.busca_impulso(h))
if len(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

..

 

3.5 LTI CT Respuesta del Sistema Y(s)=ZIR+ZSR con Sympy-Python

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.

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero
ZSR = h(t) ⊗ x(t)

Referencia: Lathi 1.8-1 p111. Oppenheim problema 2.61c p164 Ejemplo 9.24 p700FIEC05058_RLC

Para el ejemplo, se plantea determinar la corriente de lazo y(t) del circuito mostrado en la imagen.

\frac{dy(t)}{dt} +3 y(t) + 2\int_{-\infty}^t y(\tau)d\tau = x(t)

Como se había indicado en los desarrollos preliminares, para tener todo expresado con un solo operador, se derivan ambos lados de la ecuación:

\frac{d^{2}y(t)}{dt^{2}} + 3\frac{dy(t)}{dt} + 2y(t) = \frac{dx(t)}{dt}

En la entrada de sistema se aplica:

x(t) = 10 e^{-3t} \mu (t)

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.

Resultados con el Algoritmo en Sympy-Python

 ZIR(t):
     -t      -2*t
- 5*e   + 5*e    

 h(t):
/   -t      -2*t\             
\- e   + 2*e    /*Heaviside(t)

 ZSR(t):
     -t                    -2*t                    -3*t             
- 5*e  *Heaviside(t) + 20*e    *Heaviside(t) - 15*e    *Heaviside(t)
xcausal:  True
hcausal:  True
limites de integral: [ 0 , t ]

 y(t) = ZIR(t)+ZSR(t):
     -t                   -t       -2*t                   -2*t       -3*t     
- 5*e  *Heaviside(t) - 5*e   + 20*e    *Heaviside(t) + 5*e     - 15*e    *Heaviside(t)
>>> 

Graficas de y(t) = ZIR(t)+ZSR(t)

LTIC Ejercicio01 Y_total Sympy

LTIC Ejercicio 01 ZIR Sympy

LTIC Ejercicio 01 ht Sympy

LTIC Ejercicio 01 ZSR Sympy

LTIC YTotal ZIR Ej01 animado

Instrucciones en Python

# Respuesta total del sistema
# y(t) = ZIR(t) + ZSR(t)
# http://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 ZSR
if not(sol_ZSR['cond_graf']):
    print('revisar acortar x(t) o h(t) para BIBO')
    
# SALIDA
print('\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.

3.4.3 LTI CT – Respuesta a estado cero ZSR con Sympy-Python

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

Referencia Lathi ejemplo 2.9 p175, Lathi ejemplo 2.6 p166

En la entrada de sistema se aplica:

x(t) = 10 e^{-3t} \mu (t)

El sistema tiene respuesta a impulso:

h(t) = \big( 2e^{-2t} -e^{-t}\big)\mu (t)

El ejemplo es la continuación de lo realizado en respuesta a entrada cero, que tiene la ecuación:

(D^2 + 3D +2)y(t) = Dx(t)

La respuesta del sistema y(t) para un LTIC se determina como:

y(t) = x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

Si la entrada 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

y(t)=\begin{cases}\int_{0^{-}}^{t} x(\tau)h(t-\tau) \delta \tau , & t\ge 0\\ 0, & t<0 \end{cases}

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.

xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)

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 = False
    if hcausal==False and 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:

Resultados con Python

xh :
      -t  -2*tau                                     
- 10*e  *e      *Heaviside(tau)*Heaviside(t - tau) + 

    -2*t  -tau                                  
+ 20*e    *e    *Heaviside(tau)*Heaviside(t - tau)
xcausal : True
hcausal : True
[tau_a,tau_b] : [0, t]
intercambia : False
cond_graf : True
ZSR :
/     -t       -2*t       -3*t\             
\- 5*e   + 20*e     - 15*e    /*Heaviside(t)

respuesta Estado Cero 01 Sympy

Grafica animada de la convolución

LTIC ZSR Ej01 animado

Instrucciones en Python

# Respuesta estado cero ZSR con x(t) y h(t)
# Integral de convolucion con Sympy
# http://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 = 10*sym.exp(-3*t)*u
h = (2*sym.exp(-2*t)-sym.exp(-t))*u

# grafica intervalo [t_a,t_b] plano simétrico
t_b = 5 ; t_a = -t_b
muestras = 101

# PROCEDIMIENTO
def respuesta_ZSR(x,h):
    '''Respuesta a estado cero x(t) y h(t)
    '''
    # 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 = False
    if hcausal==False and 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)
    if not(ZSR.has(sym.Integral)):
        ZSR = fcnm.simplifica_escalon(ZSR)

    lista_escalon = ZSR.atoms(sym.Heaviside)
    ZSR = sym.expand(ZSR,t) # terminos suma
    ZSR = sym.collect(ZSR,lista_escalon)

    if intercambia == True:
        xcausal = True
        hcausal = False

    # graficar si no tiene Integral o error
    cond_graf = ZSR.has(sym.Integral)
    cond_graf = cond_graf or ZSR.has(sym.oo)
    cond_graf = cond_graf or ZSR.has(sym.nan)
    cond_graf = not(cond_graf)
    
    sol_ZSR = {'xh'      : xh,
               'xcausal' : xcausal,
               'hcausal' : hcausal,
               '[tau_a,tau_b]': [tau_a,tau_b],
               'intercambia'  : intercambia,
               'cond_graf'    : cond_graf,
               'ZSR' : ZSR,}
    return(sol_ZSR)

sol_ZSR = respuesta_ZSR(x,h)
ZSR = sol_ZSR['ZSR']

# SALIDA
fcnm.print_resultado_dict(sol_ZSR)
if sol_ZSR['hcausal']==False:
    print('revisar causalidad de h(t)')
if sol_ZSR['xcausal']==False:
    print('revisar causalidad de x(t)')
if sol_ZSR['intercambia']:
    print('considere intercambiar h(t)con x(t)')
if not(sol_ZSR['cond_graf']):
    print('revisar acortar x(t) o h(t) para BIBO')

# GRAFICA
if sol_ZSR['cond_graf']:
    fig_xh_y = fcnm.graficar_xh_y(x,h,ZSR,t_a,t_b,
                            muestras,y_nombre='ZSR')
#plt.show()

# grafica animada de convolución
n_archivo = '' # sin crear archivo gif animado 
# n_archivo = 'LTIC_ZSR_Ej01' # requiere 'imagemagick'
if sol_ZSR['cond_graf']:
    figura_animada = fcnm.graf_animada_xh_y(x,h,ZSR,t_a,t_b,
                      muestras, reprod_x = 4,y_nombre='ZSR',
                      archivo_nombre = n_archivo)
plt.show()

[ ejemplo 1. x(t) y h(t) causal ]  [ ejemplo 2. h(t) causal ]  [ ejemplo 3. x(t) y h(t) causal]
..


Ejemplo 2. Respuesta entrada cero ZSR entre exponencial y escalón unitario

Referencia: Oppenheim Ej 2.6 p98

Sea la y(t) la respuesta a entrada cero entre las siguientes señales:

x(t) = e^{-2t} \mu (t) h(t) = \mu (t)

Para interpretar el integral de convolución en el ejercicio aa desarrollar se adjunta la animación siguiente:

LTIC ZSR Ej02 animado

Las funciones para el algoritmo se describen como:

# entrada x(t), respuesta impulso h(t)
x = sym.exp(-2*t)*u
h = u

Teniendo como resultado:

xh :
 -2*tau                                  
e      *Heaviside(tau)*Heaviside(t - tau)
xcausal : True
hcausal : True
[tau_a,tau_b] : [0, t]
intercambia : False
cond_graf : True
ZSR :
/     -2*t\             
|1   e    |             
|- - -----|*Heaviside(t)
\2     2  /                

Gráfica de algoritmo

LTIC ZSR Ej02 Sympy

[ ejemplo 1. x(t) y h(t) causal ]  [ ejemplo 2. h(t) causal ]  [ ejemplo 3. x(t) y h(t) causal]
..


Ejemplo 3. Ejercicio con Filtros

Tarea:  Lathi 2.6-4 filtros p207

Para pruebas de algoritmo


[ ejemplo 1. x(t) y h(t) causal ]  [ ejemplo 2. h(t) causal ]  [ ejemplo 3. x(t) y h(t) causal]

3.4.2 Integral de Convolución entre x(t) y h(t) causal con Sympy-Python

Referencia: Lathi 2.4-1 p170, Ej 2.8 p173.

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.

x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

Integral de Convolucion 01 animado

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

y(t)=\begin{cases}\int_{0^{-}}^{t} x(\tau)h(t-\tau) \delta \tau , & t\ge 0\\ 0, & t<0 \end{cases}

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.

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectangulo y rampa]

..


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)
if not(y.has(sym.Integral)):
    y = fcnm.simplifica_escalon(y)
#ZSR  = ZSR.subs(u**2,u) # Heaviside**2=Heaviside

# SALIDA
print('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 -----
def graficar_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()

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectángulo y rampa]

..


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

convolucion 02 animado

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 = False
    if hcausal==False and 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

LTI C ZSR Ej03 Sympy

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)
# http://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 = False
if hcausal==False and 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)
if not(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 and not(ZSR.has(sym.oo))
cond_graf = cond_graf and not(ZSR.has(sym.nan))

# SALIDA
print('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)')
if not(cond_graf):
    print('revisar acortar x(t) o h(t) para BIBO')

# GRAFICA
if 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()

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectangulo y rampa]

..


Ejemplo 3. Convolución entre  rectangular y rampa causales

Referencia: Oppenheim Ejemplo 2.7 p99

Realizar la convolución entre las siguientes dos señales:

x(t) = \mu (t) - \mu (t-1) h(t) = t\mu (t) - t\mu (t-2)

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.
convolucion entre rectangulo y rampa truncada
El resultado con el algoritmo del ejercicio anterior es:

xcausal:  True
hcausal:  True
limites de integral: [ 0 , t ]
x(t)*h(t-tau):
t*Heaviside(tau)*Heaviside(t - tau)
+
t*Heaviside(tau - 1)*Heaviside(t - tau - 2)
+
tau*Heaviside(tau)*Heaviside(t - tau - 2)
+
tau*Heaviside(t - tau)*Heaviside(tau - 1)
+
-t*Heaviside(tau)*Heaviside(t - tau - 2)
+
-t*Heaviside(t - tau)*Heaviside(tau - 1)
+
-tau*Heaviside(tau)*Heaviside(t - tau)
+
-tau*Heaviside(tau - 1)*Heaviside(t - tau - 2)

 ZSR(t):
   2                                  2                
  t *Heaviside(t)*Heaviside(t - 1)   t *Heaviside(t)   
- -------------------------------- + --------------- + 
                 2                          2          

   2                                      2                 2   
  t *Heaviside(t - 3)*Heaviside(t - 2)   t *Heaviside(t - 2)    
+ ------------------------------------ - -------------------- + 
                   2                              2            

+ t*Heaviside(t)*Heaviside(t - 1) - t*Heaviside(t - 3)*Heaviside(t - 2) 
                                                                                           
   Heaviside(t)*Heaviside(t - 1)   3*Heaviside(t- 3)*Heaviside(t - 2) 
- ----------------------------- - ----------------------------------- + 
                 2                                  2                

                    2
+ 2*Heaviside(t - 2) 

>>>

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.

La función para simplificar escalón unitario multiplicados se desarrolla en Simplificar multiplicación entre impulso δ(t) o escalón unitario μ(t) con Sympy
De donde se obtiene la instrucción para el algoritmo,

ZSR = fcnm.simplifica_escalon(ZSR)

El resultado luego de aplicar la función sobre ZSR se muestra un poco mas simple:

 ZSR(t):
 2                 2                     2                     2              
t *Heaviside(t)   t *Heaviside(t - 3)   t *Heaviside(t - 2)   t *Heaviside(t - 1)
--------------- + ------------------- - ------------------- - -------------------
       2                   2                     2                     2      

                                                                              
                                            3*Heaviside(t - 3)   Heaviside(t - 1)
- t*Heaviside(t - 3) + t*Heaviside(t - 1) - ------------------ - ----------------
                                                    2                   2        

+ 2*Heaviside(t - 2)

y si reagrupan las expresiones por cada escalón desplazado con sym.collect() se puede obtener:

lista_escalon = ZSR.atoms(sym.Heaviside)
ZSR = sym.expand(ZSR,t) # terminos suma
ZSR = sym.collect(ZSR,lista_escalon)

con el siguiente resultado:

 ZSR(t):
 2                /     2\                    /   2        \                  
t *Heaviside(t)   |    t |                    |  t        1|                  
--------------- + |2 - --|*Heaviside(t - 2) + |- -- + t - -|*Heaviside(t - 1) 
       2          \    2 /                    \  2        2/                  

  / 2        \                 
  |t        3|                 
+ |-- - t - -|*Heaviside(t - 3)
  \2        2/                 

convolucion ejercicio 03

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectangulo y rampa]


But what is a convolution? 3Blue1Brown. 18 nov 2022