6.2.1 LTI DT – Respuesta entrada cero. Ejercicios con Python

Los siguientes ejercicios de LTI Discretos muestran los desarrollos  cuando el denominador de la función de transferencia tiene:

1. raíces reales diferentes (No repetidas)
2. raíces reales repetidas
3. raíces con números complejos


Ejercicio 1. Raíces reales diferentes

Referencia: Lathi Ejemplo 3.13 p271

Para un sistema LTID descrito por la ecuación de diferencias mostrada,

y[n+2] - 0.6 y[n+1] - 0.16 y[n] = 5x[n+2]

encuentre la respuesta a entrada cero y0[n] del sistema dadas las condiciones iniciales y[-1]=0, y[-2] =25/4, ante una entrada

x[n]= 4^{-n} \mu [n]

(El componente de estado cero se calcula en la sección correspondiente)


Para realizar el diagrama, se desplaza la ecuación en dos unidades

y[n] - 0.6 y[n-1] - 0.16 y[n-2] = 5x[n] y[n] = 0.6 y[n-1] + 0.16 y[n-2] + 5x[n]

Para respuesta a entrada cero se usa x[n] = 0

1.1 Desarrollo analítico

y[n+2] - 0.6 y[n+1] - 0.16 y[n] = 5x[n+2]

La ecuación del sistema en notación de operador E es:

E^2 y[n] - 0.6 E y[n] - 0.16 y[n] = 5 E^2 x[n] (E^2 - 0.6 E - 0.16) y[n] = 5 E^2 x[n]

que para entrada cero tiene x[n]=0,

(E^2 - 0.6 E - 0.16) y[n] = 0

tiene el polinomio característico,

\gamma^2 - 0.6 \gamma - 0.16 = 0 (\gamma+0.2) (\gamma - 0.8) = 0

con raíces características γ1=-0.2 y γ2=0.8, obtenidas con  instrucciones con Python:

raices:
[ 0.8 -0.2]
coeficientes:
[0.8 0.2] 

la repuesta a entrada cero tiene la forma

y_0[n] = c_1 (-0.2)^n + c_2 (0.8)^n

se usan las condiciones iniciales para determinar las constantes c1 y c2 aplicadas como:

y_0[-1] = c_1 (-0.2)^{-1} + c_2 (0.8)^{-1} 0 = -5 c_1 + 1.25 c_2 y_0[-2] = c_1 (-0.2)^{-2} + c_2 (0.8)^{-2} \frac{25}{4} = 25 c_1 + 1.5625 c_2

se tiene el sistema de ecuaciones

0 = -5 c_1 + 1.25 c_2 \frac{25}{4} = 25 c_1 + 1.5625 c_2

resolviendo se encuentra que: c1 = 0.2 c2 = 0.8

con lo que la solución tiene la expresión basado en la sumatoria de ci ϒn, para n≥0

y_0 [n] = \frac{1}{5} (-0.2)^n + \frac{4}{5} (0.8)^n

Instrucciones en Python

import numpy as np

coeficientes = [1.,-0.6,-0.16]
raiz = np.roots(coeficientes)
print('raices:')
print(raiz)

A = np.array([[raiz[0]**(-1) , raiz[1]**(-1)],
     [raiz[0]**(-2) , raiz[1]**(-2)]])
B = np.array([0., 25/4])
ci = np.linalg.solve(A,B)
print('coeficientes:')
print(ci)

1.2 Desarrollo con Python

El algoritmo requiere los coeficientes de la ecuación de diferencias Q(E) y P(E) junto a las condiciones iniciales.

QE = [1., -0.6, -0.16]
PE = [5., 0., 0.]
inicial = [25/4, 0.] 

que da como resultado:

respuesta entrada cero: 
raices:  [ 0.8 -0.2]
Matriz: 
[[ 1.5625 25.    ]
 [ 1.25   -5.    ]]
Ci:      [0.8 0.2]
y0:
        n          n
0.2*-0.2  + 0.8*0.8 

Instrucciones en Python

# Sistema LTID. Respuesta entrada cero
# resultado con raices Reales
# ejercicio Lathi ejemplo 3.13 p271
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
# coeficientes E con grado descendente
QE = [1., -0.6, -0.16]
PE = [5., 0., 0.]
# condiciones iniciales ascendente ...,y[-2],y[-1]
inicial = [25/4, 0.] 

# PROCEDIMIENTO
# Respuesta a ENTRADA CERO
# raices, revisa numeros complejos
gamma = np.roots(QE)
revisaImag = np.iscomplex(gamma)
escomplejo = np.sum(revisaImag)

# coeficientes de ecuacion
m_q = len(QE)-1
Ac = np.zeros(shape=(m_q,m_q),dtype=float)

if escomplejo == 0:
    for i in range(0,m_q,1):
        for j in range(0,m_q,1):
            Ac[i,j] = gamma[j]**(-m_q+i)
    ci = np.linalg.solve(Ac,inicial)

# ecuacion y0
n = sym.Symbol('n')
y0 = 0
if escomplejo == 0:
    for i in range(0,m_q,1):
        y0 = y0 + ci[i]*(gamma[i]**n)
    
# SALIDA
print('respuesta entrada cero: ')
print('raices: ', gamma)
if escomplejo == 0:
    print('Matriz: ')
    print(Ac)
    print('Ci:     ', ci)
    print('y0:')
    sym.pprint(y0)
else:
    print(' existen raices con números complejos.')
    print(' usar algoritmo de la sección correspondiente.')

Raíces: [reales No repetidas] [reales repetidas] [números complejos]


2. Raíces repetidas

Referencia: Lathi Ejemplo 3.15 p274

y[n+2] +6 y[n+1] + 9 y[n] = 2 x[n+2] + 6 x[n+1]

con las condiciones iniciales y0[-1] = -1/3 y y0[-2] = -2/9

Para realizar el diagrama de bloques, se desplaza ambos lados de la ecuación en 2 unidades. Luego se despeja y[n]

y[n] +6 y[n-1] + 9 y[n-2] = 2 x[n] + 6 x[n-1] y[n] = -6 y[n-1] - 9 y[n-2] + 2 x[n] + 6 x[n-1]

2.1 Desarrollo analítico

siguiendo los mismos pasos que en los ejercicios anteriores

E^2 y[n] +6 E y[n] + 9 y[n] = 2 E^2 x[n] + 6 E x[n] (E^2 +6 E + 9) y[n] = (2 E^2 + 6 E) x[n]

Para entrada cero se tiene x[n] = 0

(E^2 +6 E + 9) y[n] = (2 E^2 + 6 E) (0) \gamma^2 +6 \gamma + 9 = 0 (\gamma + 3)( \gamma + 3) = 0

de donde se muestra una raiz repetida γ=3, donde los modos característicos hacen que la expresión sea:

y_0[n] = (c_1 + c_2 n) (\gamma)^n

siguiendo el procedimiento realizado para raíces reales se puede determinar los valores de c1=4 y c2=3, con lo que se obtiene la solución:

y_0[n] = (4 + 3n) (-3)^n

los resultados obtenidos con el algoritmo en Python son

respuesta entrada cero: 
raices:  [-3. -3.]
Raices repetidas:  1
Matriz: 
[[ 0.11111111 -0.22222222]
 [-0.33333333  0.33333333]]
Ci:      [4. 3.]
y0:
    n              
-3.0 *(3.0*n + 4.0)
>>>

2.2 Instrucciones en Python, integrando raices reales únicas y repetidas

# Sistema LTID. Respuesta entrada cero
# QE con raices Reales NO repetidas Lathi ejemplo 3.13 pdf271
# QE con raices Reales Repetidas Lathi ejemplo 3.15 p274
# blog.espol.edu.ec/telg1001
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
# coeficientes E con grado descendente
QE = [1., 6., 9.]
PE = [2., 6., 0.]
# condiciones iniciales ascendente ...,y[-2],y[-1]
 inicial = [-2/9, -1/3]

tolera = 1e-6  # casi_cero

# PROCEDIMIENTO
# Respuesta a ENTRADA CERO
# raices, revisa numeros complejos
gamma = np.roots(QE)
revisaImag = np.iscomplex(gamma)
escomplejo = np.sum(revisaImag)

# coeficientes de ecuacion
m_q = len(QE)-1
Ac = np.zeros(shape=(m_q,m_q),dtype=float)

# revisa si parte compleja <tolera o casi_cero 
if escomplejo>0:
    for i in range(0,m_q,1):
        if np.imag(gamma[i]) < tolera:
            gamma[i] = float(np.real(gamma[i]))
    sumaj = np.sum(np.imag(gamma))
    if sumaj <tolera:
        gamma = np.real(gamma)
        escomplejo = 0

# revisa repetidos
unicoscuenta = np.unique(gamma,return_counts=True)
repetidas = np.sum(unicoscuenta[1]-1)

# raices Reales no repetidas
if escomplejo == 0 and repetidas==0:
    for i in range(0,m_q,1):
        for j in range(0,m_q,1):
            Ac[i,j] = gamma[j]**(-m_q+i)
    ci = np.linalg.solve(Ac,inicial)
    
# raices Reales repetidas
if escomplejo == 0 and repetidas > 0:
    for i in range(0,m_q,1):
        for j in range(0,m_q,1):
            Ac[i,j] = ((-m_q+i)**j)*gamma[j]**(-m_q+i)
    ci = np.linalg.solve(Ac,inicial)

# ecuacion y0 entrada cero
n = sym.Symbol('n')
y0 = 0*n

if escomplejo == 0 and repetidas==0:
    for i in range(0,m_q,1):
        y0 = y0 + ci[i]*(gamma[i]**n)
        
if escomplejo == 0 and repetidas > 0:
    for i in range(0,m_q,1):
        y0 = y0 + ci[i]*(n**i)*(gamma[i]**n)
    y0 = y0.simplify()
    
# SALIDA
print('respuesta entrada cero: ')
print('raices: ', gamma)
if escomplejo == 0:
    if repetidas>0:
        print('Raices repetidas: ', repetidas)
    print('Matriz: ')
    print(Ac)
    print('Ci:     ', ci)
    print('y0:')
    sym.pprint(y0)
else:
    print(' existen raices con números complejos.')
    print(' usar algoritmo de la sección correspondiente.')

Raíces: [reales No repetidas] [reales repetidas] [números complejos]


Ejercicio 3. Raices con números complejos

Referencia: Lathi Ejemplo 3.16 p275

Para el caso de tener raíces con números complejos se muestra el siguiente ejemplo:

y[n+2] - 1.56 y[n+1] + 0.81 y[n] = x[n+1] + 3x[n]

con condiciones iniciales y0[-1]=2 y y0[-2]=1

Para realizar el diagrama de bloques, se desplaza ambos lados de la ecuación en 2 unidades. Luego se despeja y[n]

y[n] - 1.56 y[n-1] + 0.81 y[n-2] = x[n-1] + 3x[n-2] y[n] = 1.56 y[n-1] - 0.81 y[n-2] + x[n-1] + 3x[n-2]

3.1 Desarrollo analítico

y[n+2] - 1.56 y[n+1] + 0.81 y[n] = x[n+1] + 3x[n]

La expresión usando operador E:

E^2 y[n] - 1.56 E y[n] + 0.81 y[n] = E x[n] + 3x[n] (E^2 - 1.56 E + 0.81) y[n] = (E+ 3)x[n]

para entrada cero se hace x[n] = 0

(E^2 - 1.56 E + 0.81) y[n] = (E+ 3)(0)

por lo que, el polinomio característico es,}

(\gamma^2 - 1.56 \gamma + 0.81) = 0

con raíces de números complejos

(\gamma - (0.78 + j 0.45))(\gamma - (0.78-j 0.45)) = 0

la raíz también se puede escribir como:

\gamma = 0.9 e^{\pm j\frac{\pi}{6}}

obtenida mediante:

raices:
[0.78+0.44899889j 0.78-0.44899889j]
magnitud:  [0.9 0.9]
angulo:  [ 0.52231482 -0.52231482]
>>>  

Se puede escribir la solución como:

y_0[n] = c(0.9)n e^{j\frac{\pi}{6} } + c^{*} (0.9)n e^{-j\frac{\pi}{6}}

Luego aplicando oas condiciones iniciales  se encuentra que

c = 2.34 e^{-j0.17} c^{*} = 2.34 e^{j0.17}

Otra manera de desarrollo permite encontrar la solución usando la forma real de la solución, conociendo que la solución esta dada por la expresión:

y_0[n] = c(0.9)^n \cos \Big( \frac{\pi}{6}n + \theta \Big)

con lo que se puede crear el sistema de ecuaciones aplicando las condiciones iniciales, ademas de considerar que:

\cos(a \pm b) = \cos(a) \cos(b) \mp \sin(a) \sin(b)

Para y0[-2] = 1

1 = c(0.9)^{-2} \cos \Big( \frac{\pi}{6}(-2) + \theta \Big) 1 = c\frac{1}{(0.9)^2} \Big[ \cos \big(-2\frac{\pi}{6} \big) \cos ( \theta ) - \sin \big(-2\frac{\pi}{6}\big) \sin (\theta) \Big] 1 = c\frac{1}{0.81} \Big[ \frac{1}{2} \cos ( \theta ) - \frac{-\sqrt{3}}{2} \sin (\theta) \Big] 1 = c \Big[ \frac{1}{1.62} \cos ( \theta ) + \frac{\sqrt{3}}{1.62} \sin (\theta) \Big]

Para y0[-1] = 2

2 = c(0.9)^{-1} \cos \Big( \frac{\pi}{6}(-1) + \theta \Big) 2 = c\frac{1}{0.9} \Big[ \cos \big(-\frac{\pi}{6} \big) \cos ( \theta ) - \sin \big(-\frac{\pi}{6}\big) \sin (\theta) \Big] 2 = c\frac{1}{0.9} \Big[ \frac{\sqrt{3}}{2} \cos ( \theta ) - \frac{-1}{2} \sin (\theta) \Big] 2 = c\Big[ \frac{\sqrt{3}}{1.8} \cos ( \theta ) + \frac{1}{1.8} \sin (\theta) \Big]

con lo que se puede determinar los valores para [c cos(θ)] y [c sin(θ)]

usando una matriz A y vector B a partir de las ecuaciones

c \cos ( \theta ) = 2.31 c \sin ( \theta ) = -0.4049

al dividir  [c sin(θ)]/[c cos(θ)] = tan (θ), con lo que se obtiene el ángulo en radianes.

\frac{ c \sin ( \theta )} {c \cos ( \theta )} = \frac{-0.4049}{2.31} = -0.17528 \theta = \tan ^{-1} (-0.17528) = -0.1735

al sustituir el valor en c cos(θ) = 2.31 se tiene,

c = \frac{2.31}{ \cos (-0.1735)} = 2.3452

con lo que la solución se puede escribir mediante:

y_0[n] = 2.3452 (0.9)^n \cos \Big( \frac{\pi}{6}n - 0.1735 \Big)

para n≥0

Instrucciones de cálculo

import numpy as np

coeficientes = [1.,-1.56,0.81]
raiz = np.roots(coeficientes)
print('raices:')
print(raiz)

g_magnitud = np.absolute(raiz)
g_angulo = np.angle(raiz)

print('magnitud: ',g_magnitud)
print('angulo: ',g_angulo)

3.2 Algoritmo en Python, integrando raices reales y complejas

usando el algoritmo se obtienen los siguientes resultados:

Respuesta entrada cero: 
raices :  [0.78+0.44899889j 0.78-0.44899889j]
raices complejas:  2
magnitud: [0.9 0.9]
theta: [ 0.52231482 -0.52231482]
Matriz: 
[[0.62002743 1.06757851]
 [0.96296296 0.55431961]]
Cj:      [ 2.31       -0.40490078] 
Ci:  2.345217397781524
y0: 
                    n                                             
2.34521739778152*0.9 *cos(0.522314821806049*n + 0.173519005551295)
>>> 

Para el caso de raíces complejas, se añade al algoritmo anterior las instrucciones en Python:

# Sistema LTID. Respuesta entrada cero
# QE con raices Reales NO repetidas Lathi ejemplo 3.13 pdf271
# QE con raices Reales Repetidas Lathi ejemplo 3.15 p274
# QE con raices Complejas Lathi ejemplo 3.16 p275
# blog.espol.edu.ec/telg1001
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
# coeficientes E con grado descendente
QE = [1.0, -1.56, 0.81]
PE = [0., 1., 3.]
# condiciones iniciales ascendente ...,y[-2],y[-1]
inicial = [1.,2.]

tolera = 1e-6  # casi_cero
muestras = 10  # para grafica

# PROCEDIMIENTO
# Respuesta a ENTRADA CERO
# raices, revisa numeros complejos
gamma = np.roots(QE)
revisaImag = np.iscomplex(gamma)
escomplejo = np.sum(revisaImag)

# coeficientes de ecuacion
m_q = len(QE)-1
Ac = np.zeros(shape=(m_q,m_q),dtype=float)

# revisa si parte compleja <tolera o casi_cero 
if escomplejo>0:
    for i in range(0,m_q,1):
        valorimag = np.imag(gamma[i])
        if np.abs(valorimag)<tolera:
            gamma[i] = float(np.real(gamma[i]))
    sumaj = np.sum(np.abs(np.imag(gamma)))
    if sumaj <tolera:
        print(sumaj)
        gamma = np.real(gamma)
        escomplejo = 0

# revisa repetidos
unicoscuenta = np.unique(gamma,return_counts=True)
repetidas = np.sum(unicoscuenta[1]-1)

# Determina coeficientes ci de Y[n]
# raices Reales no repetidas
if escomplejo == 0 and repetidas==0:
    for i in range(0,m_q,1):
        for j in range(0,m_q,1):
            Ac[i,j] = gamma[j]**(-m_q+i)
    ci = np.linalg.solve(Ac,inicial)
    
# raices Reales repetidas
if escomplejo == 0 and repetidas > 0:
    for i in range(0,m_q,1):
        for j in range(0,m_q,1):
            Ac[i,j] = ((-m_q+i)**j)*gamma[j]**(-m_q+i)
    ci = np.linalg.solve(Ac,inicial)

# raices Complejas
if escomplejo > 0:
    g_magnitud = np.absolute(gamma)
    g_angulo = np.angle(gamma)

    for i in range(0,m_q,1):
        k = -(m_q-i)
        a = np.cos(np.abs(g_angulo[i])*(k))*(g_magnitud[i]**(k))
        b = -np.sin(np.abs(g_angulo[i])*(k))*(g_magnitud[i]**(k))
        Ac[i] = [a,b]
    Ac = np.array(Ac)
    cj = np.linalg.solve(Ac,inicial)

    theta = np.arctan(cj[1]/cj[0])
    ci = cj[0]/np.cos(theta)

# ecuacion y0 entrada cero
n = sym.Symbol('n')
y0 = 0*n

if escomplejo == 0 and repetidas==0:
    for i in range(0,m_q,1):
        y0 = y0 + ci[i]*(gamma[i]**n)
        
if escomplejo == 0 and repetidas > 0:
    for i in range(0,m_q,1):
        y0 = y0 + ci[i]*(n**i)*(gamma[i]**n)
    y0 = y0.simplify()

if escomplejo > 0:
    y0 = ci*(g_magnitud[0]**n)*sym.cos(np.abs(g_angulo[i])*n - theta)
    
# SALIDA
print('respuesta entrada cero: ')
print('raices: ', gamma)
if escomplejo == 0:
    if repetidas>0:
        print('Raices repetidas: ', repetidas)
    print('Matriz: ')
    print(Ac)
    print('Ci:     ', ci)
    print('y0:')
    sym.pprint(y0)
if escomplejo > 0:
    print('raices complejas: ', escomplejo)
    print('magnitud:',g_magnitud)
    print('theta:',g_angulo)
    print('Matriz: ')
    print(Ac)
    print('Cj: ', cj)
    print('Ci: ',ci)
    print('y0: ')
    sym.pprint(y0)

# grafica datos
ki = np.arange(-m_q,muestras,1)
y0i = np.zeros(muestras+m_q)
# valores iniciales
y0i[0:m_q] = inicial[0:m_q]
# evaluación de y0[n]
y0n = sym.lambdify(n,y0)
y0i[m_q:] = y0n(ki[m_q:])

# grafica y[n] ante x[n]=0
plt.stem(ki,y0i)
plt.xlabel('ki')
plt.ylabel('y0[n]')
plt.title('y0[n]='+str(y0))
plt.grid()
plt.show()

Raíces: [reales No repetidas] [reales repetidas] [números complejos]

6.1.3 Sistema LTI DT – Solución iterativa con Python

El ejercicio planteado en Sistema Discreto – Solución iterativa, se desarrolla con un algoritmo en Python.


Ejercicio 1

Referencia: Lathi ejemplo 3.11 p266

Un sistema LTID tiene la siguiente expresión:

y[n] - 0.5 y[n-1] = x[n]

con condiciones iniciales y[-1]=16 y entrada causal x[n] = n2 ,para n≥0.


El algoritmo tiene el siguiente resultado:

LTI DT Solución iterativaEj01

muestras:  10
[[  n,  xi,  yi]]
[[ -1.     0.    16.  ]
 [  0.     0.     8.  ]
 [  1.     1.     5.  ]
 [  2.     4.     6.5 ]
 [  3.     9.    12.25]
 [  4.    16.    22.12]
 [  5.    25.    36.06]
 [  6.    36.    54.03]
 [  7.    49.    76.02]
 [  8.    64.   102.01]
 [  9.    81.   132.  ]]
>>> 

Algoritmo en Python

Las condiciones iniciales se ingresan como un vector con valores de posición ascendente. También es necesario indicar el número de muestras a observar. esto permite construir el vector de salida yi.

Las fórmulas de entrada y salida se muestran en la sección «calcula los siguientes valores»

# Sistema LTID. Solución iterativa
# Ejercicio Lathi ejemplo 3.11 p266
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
muestras = 10
# condiciones iniciales ascendente ...,y[-2],y[-1]
y0  = [16.]

# PROCEDIMIENTO
m0 = len(y0)
m  = muestras + m0
# vectores
n  = np.arange(-m0,muestras,1)
xi = np.zeros(m, dtype=float)
yi = np.zeros(m, dtype=float)

# añade condiciones iniciales
xi[0:m0] = 0
yi[0:m0] = y0

# calcula los siguientes valores
for k in range(m0,m,1):
    xi[k] = n[k]**2
    yi[k] = 0.5*yi[k-1]+xi[k]

# tabla de valores
tabla = np.concatenate([[n],[xi],[yi]], axis=0)
tabla = np.transpose(tabla)
# SALIDA
np.set_printoptions(precision=2)
print('muestras: ', muestras)
print('[[  n,  xi,  yi]]')
print(tabla)

# Gráfica
plt.subplot(211)
plt.stem(n,xi,label='x[n]=n**2')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.axvline(0,color='gray')
plt.legend()
plt.grid()

plt.subplot(212)
plt.stem(n,yi,label='y[n]=0.5y[n-1]+x[n]',
         markerfmt='go',
         linefmt='green')
plt.xlabel('n')
plt.ylabel('y[n]')
plt.axvline(0,color='gray')
plt.legend()
plt.grid()
plt.show()

Ejercicio 2

Se aplica el mismo algoritmo a otro ejercicio para comparar resultados con lo realizado en forma analítica.

Referencia: Lathi ejemplo 3.12 p267

Resolver de forma iterativa,

y[n+2] - y[n+1] +0.24y[n] = x[n+2] - 2x[n+1]

con las condiciones iniciales de y[-1] = 2, y[-2] =1 y entrada causal x[n]=n que inicia en n=0

Para usar la ecuación en el algoritmo, se desplaza dos posiciones en ambos lados de la ecuacion para despejar y[n]

y[n] - y[n-1] +0.24y[n-2] = x[n0] - 2x[n-1] y[n] = y[n-1] -0.24y[n-2] + x[n0] - 2x[n-1]

LTI DT Sol iterativa Ej02

muestras:  10
[[  n,  xi,  yi]]
[[ -2.     0.     1.  ]
 [ -1.     0.     2.  ]
 [  0.     0.     1.76]
 [  1.     1.     2.28]
 [  2.     2.     1.86]
 [  3.     3.     0.31]
 [  4.     4.    -2.14]
 [  5.     5.    -5.21]
 [  6.     6.    -8.7 ]
 [  7.     7.   -12.45]
 [  8.     8.   -16.36]
 [  9.     9.   -20.37]]
>>> 

Algoritmo en Python

# Sistema LTID. Solución iterativa
# Ejercicio Lathi ejemplo 3.12 p267
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
muestras = 10
# condiciones iniciales ascendente ...,y[-2],y[-1]
y0  = [1,2]

# PROCEDIMIENTO
m0 = len(y0)
m  = muestras + m0
# vectores
n  = np.arange(-m0,muestras,1)
xi = np.zeros(m, dtype=float)
yi = np.zeros(m, dtype=float)

# añade condiciones iniciales
xi[0:m0] = 0
yi[0:m0] = y0

# calcula los siguientes valores
for k in range(m0,m,1):
    xi[k] = n[k]
    yi[k] = yi[k-1]-0.24*yi[k-2]+xi[k]-2*xi[k-1]

# tabla de valores
tabla = np.concatenate([[n],[xi],[yi]], axis=0)
tabla = np.transpose(tabla)
# SALIDA
np.set_printoptions(precision=2)
print('muestras: ', muestras)
print('[[  n,  xi,  yi]]')
print(tabla)

# Gráfica
plt.subplot(211)
plt.stem(n,xi,label='x[n]=n')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.axvline(0,color='gray')
plt.legend()
plt.grid()

plt.subplot(212)
plt.stem(n,yi,label='y[n]=yi[k-1]-0.24*yi[k-2]+xi[k]-2*xi[k-1]',
         markerfmt='go',
         linefmt='green')
plt.xlabel('n')
plt.ylabel('y[n]')
plt.axvline(0,color='gray')
plt.legend()
plt.grid()
plt.show()

5.3 LTI-CT Fourier – H(ω) Funciones de transferencia

Referencia: Oppenheim 3.2 p182, 4.4 p314, Lathi p714 Hsu 5.5.A p223

Las transformadas de Fourier presentan una facilidad con la propiedad de la convolución expresada como:

y(t)= h(t) \circledast x(t) \leftrightarrow Y(\omega) = H(\omega) X(\omega)

Lo que para señales y sistemas representa un componente importante en el análisis de sistemas LTI, al poder representar los sistemas por bloques en serie y paralelo semejante a lo realizado con la transformada de Laplace:

LIT CT Bloques Serie

La convergencia de la transformada de Fourier se garantiza bajo ciertas condiciones. El análisis de Fourier para el estudio de sistemas LTI se restringe a los que sus respuestas al impulso tienen transformada de Fourier. Los detalles se mostrarán en cada ejercicio presentado.


Ejemplo 1. H(ω) desde una ecuación diferencial

Referencia: Oppenheim Ej 4.25 p331

Considere un sistema LTI estable que se caracteriza por la ecuación diferencial,

\frac{\delta^2}{\delta t^2}y(t) + 4\frac{\delta}{\delta t} y(t) + 3 y(t) = \frac{\delta}{\delta t}x(t)+2x(t)

usando la propiedad de la diferenciación en el dominio ω, se tiene que:

(j\omega)^2 Y(\omega) + 4(j\omega) Y(\omega) + 3 Y(\omega) = (j\omega)X(\omega)+2X(\omega) \Big((j\omega)^2 + 4(j\omega) + 3\Big) Y(\omega) = \Big((j\omega)+2\Big) X(\omega)

la función de transferencia se obtiene como,

H(\omega) = \frac{Y(\omega)}{X(\omega)} H(\omega) = \frac{(j\omega)+2}{(j\omega)^2 + 4(j\omega) + 3}

Para facilitar encontrar h(t) se usan fracciones parciales respecto a jω, de forma semejante a lo realizado para el dominio s con instrucciones de Sympy-Python, Las instrucciones las puede recordar observando el ejercicio 2, donde se adjunta el algoritmo.

Hw:
(jw+2)/(jw**2+4*jw+3)
Hw en fracciones parciales:
    1            1     
---------- + ----------
2*(jw + 3)   2*(jw + 1)

con la tabla de transformadas de Fourier se obtiene que,

h(t) = \frac{1}{2} e^{-t} \mu(t) + \frac{1}{2} e^{-3t} \mu(t)

Ejemplo 2. respuesta total Y(ω)=X(ω)H(ω)

Referencia: Oppenheim Ej 4.27 p332

Considere el sistema del ejemplo anterior y suponga una entrada x(t) dada por

x(t) = e^{-t} \mu(t) X(\omega) = \frac{1}{j\omega +1}

la salida del sistema usando los componentes en el dominio de la frecuencia

Y(\omega) = H(\omega)X(\omega) = \Big[ \frac{j\omega+2}{(j\omega+1)(j\omega + 3)} \Big]\Big[\frac{1}{j\omega +1}\Big] = \frac{j\omega+2}{(j\omega+1)^2 (j\omega + 3)}

separando en fracciones parciales

Hw:
          jw + 2         
-------------------------
         /  2           \
(jw + 1)*\jw  + 4*jw + 3/

Hw fracciones parciales jw:
      1            1             1     
- ---------- + ---------- + -----------
  4*(jw + 3)   4*(jw + 1)             2
                            2*(jw + 1) 
>>> 

con la tabla de transformadas de Fourier se obtiene que,

y(t) = \Big[ \frac{1}{4} e^{-t} +\frac{1}{2}te^{-t} - \frac{1}{4} e^{-3t} \Big] \mu (t)

Instrucciones en Python

# H(w) Funciones de transferencia
# blog.espol.edu.ec/telg1001/
import sympy as sym

# INGRESO
w = sym.Symbol('w', real=True)
j = sym.I
jw = sym.Symbol('jw',real=True)

Xw = 1/(jw+1)
Hw = (jw+2)/(jw**2+4*jw+3)

# PROCEDIMIENTO
Yw = Hw*Xw
Ywp = sym.apart(Yw,jw)

print('Hw:')
sym.pprint(Yw)
print('\nHw fracciones parciales jw:')
sym.pprint(Ywp)

Ejemplo 3. y(t) con h(t) y x(t) en dominio de frecuencia

Referencia: Oppenheim Ej 4.19 p320,  Hsu 5.4.I p220

Considere la respuesta al impulso de un sistema LTI como:

h(t) = e^{-at} \mu (t) \text{ , } a>0

y una señal de entrada:

x(t) = e^{-bt} \mu (t) \text{ , } b>0

En lugar de calcular la convolución se resolverá el problema en el dominio de la frecuencia.

Para observar una gráfica, se supondrán valores de a=3 y b=2

H(\omega) = \frac{1}{a+ j \omega} X(\omega) = \frac{1}{b+ j \omega} Y(\omega) = H(\omega)X(\omega) = \Big[ \frac{1}{a+ j \omega} \Big] \Big[ \frac{1}{b+ j \omega} \Big]

Para hacer la transformada inversa de Fourier, se usa fracciones parciales, semejante a lo realizado en la unidad 4

Y(\omega) = \frac{k_1}{a+ j \omega} + \frac{k_2}{b+ j \omega} k_1 = \frac{1}{\cancel{(a+ j \omega)} (b+ j \omega)} \Big|_{j\omega=-a} = \frac{1}{b-a} k_2 = \frac{1}{(a+ j \omega) \cancel{(b+ j \omega)}} \Big|_{j\omega=-b} = \frac{1}{a-b}

con lo que k1= -k2

Y(\omega) = \frac{1}{b-a} \Big[\frac{1}{a+ j \omega} - \frac{1}{b+ j \omega} \Big]

la transformada inversa se obtiene para cada término de la suma como:

y(t) = \frac{1}{b-a} \Big [ e^{-at} \mu(t) - e^{-bt} \mu(t) \Big]

Resultado usando el algoritmo del ejercicio 3

 Y(w) = H(w)*X(w)
         1         
-------------------
(a + I*w)*(b + I*w)

 Y(w) en fracciones parciales:
        1                   1        
----------------- - -----------------
(a - b)*(b + I*w)   (a - b)*(a + I*w)

Ejemplo 4. h(t) como respuesta a un impulso desplazado

Referencia: Oppenheim Ej 4.16 p317, Lathi 7.14 708. Hsu 5.4.B p219

Considere un sistema LTI CT (continuo en el tiempo) con respuesta a impulso dado por:

h(t) = \delta (t-t_0)

La respuesta a este sistema esta dada por la respuesta desplazada del impulso.

H(\omega) = 1 e^{-j \omega t_0} = e^{-j \omega t_0}

El sistema ante una entrada x(t) con transformada de Fourier X(ω) tendrá la salida:

Y(\omega) =H(\omega)X(\omega) = e^{-j \omega t_0} X( \omega)

resultado con algoritmo Python

 h(t):
DiracDelta(-a + t)

 H(w):
 -I*a*w
e      
>>> 

Instrucciones en Python

# H(w) Funciones de transferencia
# blog.espol.edu.ec/telg1001/
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True,)
w = sym.Symbol('w', real=True)
j = sym.I
a = sym.Symbol('a', real=True,positive=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)

#ht = d
ht = d.subs(t,t-a)
#ht = sym.exp(-a*sym.Abs(t))
#ht = u.subs(t,t+a)-u.subs(t,t-a)
#ht = sym.exp(-a*t)*u

# PROCEDIMIENTO
ht = sym.expand(ht)  # expresion de sumas
Hw = sym.fourier_transform(ht,t,w/(2*sym.pi))

# SALIDA
print('\n h(t):')
sym.pprint(ht)
print('\n H(w):')
sym.pprint(Hw)

Ejemplo 5. h(t) de un exponencial decreciente

Referencia: Oppenheim Ejercicio 4.1. p290, Lathi ejemplo 7.1 p685, Hsu Ejemplo 5.2 p218

x(t) =e^{-at} \mu (t) \text{ ; } a \gt 0


realizado con el algoritmo anterior, se puede comprobar el resultado con la tabla de transformadas de Fourier

 h(t):
 -t              
 ---             
  2              
e   *Heaviside(t)

 H(w):
    2    
---------
2*I*w + 1
>>>

Ejemplo 6. h(t) como un diferenciador

Referencia: Oppenheim Ej 4.16 p317, Lathi 7.14 716. Hsu 5.4.G p220

El bloque de respuesta a impulso es un diferenciador, por lo que las funciones de entrada y salida se relacionan por

y(t) = \frac{\delta}{\delta t}x(t)

usando la diferenciación de las propiedades de la transformada de Fourier

Y(\omega) = j \omega X( \omega )

en consecuencia:

H(\omega) = j \omega

Ejemplo 7. h(t) como un integrador

Referencia: Oppenheim Ej 4.17 p318,  Hsu 5.4.I p220

El bloque de respuesta a impulso es una integración, por lo que las funciones de entrada y salida se relacionan por

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

usando la propiedad de Integración:

H(\omega) = \frac{1}{j \omega} +\pi \delta(\omega) Y(\omega) = \Big[\frac{1}{j \omega} +\pi \delta(\omega) \Big]X( \omega) Y(\omega) = \frac{1}{j \omega}X(\omega) +\pi \delta(\omega) X(0)

5.2 Transformada de Fourier – Señales aperiódicas continuas con Sympy-Python

Referencia: Oppenheim 4.1. p285, Lathi 7.2 p680

El planteamiento de Fourier es que una señal aperiódica (no periódica) puede observarse como una señal periódica con periodo infinito.

[ exponencial decreciente ] [ exponencial decreciente con |t| ] [ rectangular ] [ Pulso unitario ]


Ejercicio 1. exponencial decreciente con μ(t)

Referencia: Oppenheim Ejercicio 4.1 p290, Lathi ejemplo 7.1 p685, Hsu Ejemplo 5.2 p218

Considere la señal contínua exponencial decreciente, desarrolle la transformada de Fourier

x(t) =e^{-at} \mu (t) \text{ ; } a \gt 0

el integral a desarrollar,

X(\omega) =\int_{-\infty}^{\infty} e^{-at} \mu (t) e^{-j \omega t} \delta t X(\omega) =\int_{0}^{\infty} e^{-(a+j\omega)t} \delta t =-\frac{1}{a+j\omega} e^{-(a+j\omega)t} \Big|_{0}^{\infty} X(\omega) =-\frac{1}{a+j\omega} e^{-(a+j\omega)(\infty)} +\frac{1}{a+j\omega} e^{-(a+j\omega)(0)} X(\omega) = \frac{1}{a+j\omega}

para a>0

Se obtiene una parte real y otra imaginaria al evaluar F(w) en un intervalo que incluya -a y a.

También es posible observar la magnitud y fase de F(w) en una gráfica.

Se observa por ejemplo que el valor de la magnitud |F(0)| = 1/a. También que la fase se encuentra acotada en π/2 y que fase(F(-a)) = -π/4.

Usando el algoritmo con Python se obtiene el siguiente resultado:

 expresion f(t):
 -a*t             
e    *Heaviside(t)

 expresion F(w):
   1   
-------
a + I*w

 |F(w)|:
     1      
------------
   _________
  /  2    2 
\/  a  + w  

 F(w)_fase:
     /w\
-atan|-|
     \a/
>>>

Instrucciones en Python

# Transformada de Fourier de señales aperiodicas
# blog.espol.edu.ec/telg1001/
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True,)
w = sym.Symbol('w', real=True)
a = sym.Symbol('a', real=True,positive=True)
T = sym.Symbol('T', real=True,positive=True)
j = sym.I
u = sym.Heaviside(t)
d = sym.DiracDelta(t)

ft = sym.exp(-a*t)*u
#ft = sym.exp(-a*sym.Abs(t))
#ft = sym.Heaviside(t+T)-sym.Heaviside(t-T)
#ft = d

# PROCEDIMIENTO
unilateral = False
ftw = ft*sym.exp(-j*w*t) # f(t,w) para integrar
ftw = sym.expand(ftw)  # expresion de sumas
ftw = sym.powsimp(ftw) # simplifica exponentes

lim_a = 0 # unilateral
if not(unilateral):
    lim_a = -sym.oo
    
# integral de Fourier
Fw_F = sym.integrate(ftw,(t,lim_a,sym.oo))
if Fw_F.is_Piecewise:
    Fw = Fw_F.args[0] # primera ecuacion e intervalo
    Fw = Fw[0] # solo expresion
else:
    Fw = Fw_F
Fw = Fw.simplify()

# Magnitud y Fase
Fw_magn = sym.Abs(Fw)
Fw_fase = sym.atan(sym.im(Fw)/sym.re(Fw))

# SALIDA
print('\n expresion f(t):')
sym.pprint(ft)
print('\n expresion F(w):')
sym.pprint(Fw)
print('\n |F(w)|:')
sym.pprint(Fw_magn)
print('\n F(w)_fase:')
sym.pprint(Fw_fase)

Instrucción simplificada con Sympy

La librería Sympy incorpora una función para determinar en forma simbólica la transformada de Fourier, se usa la instrucción para simplificar el algoritmo anterior.

# Transformada de Fourier de señales aperiodicas
# blog.espol.edu.ec/telg1001/
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True,)
w = sym.Symbol('w', real=True)
a = sym.Symbol('a', real=True,positive=True)
T = sym.Symbol('T', real=True,positive=True)
j = sym.I
u = sym.Heaviside(t)
d = sym.DiracDelta(t)

ft = sym.exp(-a*t)*u
#ft = sym.exp(-a*sym.Abs(t))
#ft = u.subs(t,t+T)-u.subs(t,t-T)
#ft = d

# PROCEDIMIENTO
ft = sym.expand(ft)  # expresion de sumas
Fw = sym.fourier_transform(ft,t,w/(2*sym.pi))

# Magnitud y Fase
Fw_magn = sym.Abs(Fw)
Fw_fase = sym.atan(sym.im(Fw)/sym.re(Fw))

# SALIDA
print('\n expresion f(t):')
sym.pprint(ft)
print('\n expresion F(w):')
sym.pprint(Fw)
print('\n |F(w)|:')
sym.pprint(Fw_magn)
print('\n F(w)_fase:')
sym.pprint(Fw_fase)

Instrucciones para añadir la gráfica

Para las gráficas, se dan valores a las constantes a, T, se define el intervalo de la gráfica y el número de muestras a usar.

Luego se sustituye en las ecuaciones resultantes los valores de ‘a‘ con ‘a_k‘ obteniendo valores para las gráficas. La gráfica para f(t) se realiza de forma semejante en la unidad 1. Para la gráfica F(w) se usan una gráfica con parte Real e Imaginaria, otra gráfica para magnitud y fase.

# GRAFICA ----------------
import numpy as np
import matplotlib.pyplot as plt
import telg1001 as fcnm
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
a_k = 2 # valor de 'a' constante
# Grafica, intervalo tiempo [t_a,t_b]
T_k = 1    # valor T o periodo
t_a = -2*T_k ; t_b = 2*T_k
muestras = 52     # 51 resolucion grafica

ft = ft.subs({a:a_k,T:T_k}) # a tiene valor a_k
Fw = Fw.subs({a:a_k,T:T_k})
Fw_magn = Fw_magn.subs({a:a_k,T:T_k})
Fw_fase = Fw_fase.subs({a:a_k,T:T_k})

figura_ft = fcnm.graficar_ft(ft,t_a,t_b,muestras)

# F(w) real e imaginaria
w_a = -a_k*4 ; w_b = a_k*4
wi  = np.linspace(w_a,w_b,muestras)

# convierte a sympy una constante
Fw = sym.sympify(Fw,w) 
if Fw.has(w): # no es constante
    F_w = sym.lambdify(w,Fw,
                       modules=equivalentes)
else:
    F_w = lambda w: Fw + 0*w
Fwi = F_w(wi) # evalua wi

# F(w) magnitud y fase
# convierte a sympy una constante
Fw_magn = sym.sympify(Fw_magn,w)
if Fw_magn.has(w): # no es constante
    F_w_magn = sym.lambdify(w,Fw_magn,
                            modules=equivalentes)
else:
    F_w_magn = lambda w: Fw_magn + 0*w

# convierte a sympy una constante
Fw_fase = sym.sympify(Fw_fase,w) 
if Fw_fase.has(w): # no es constante
    F_w_fase = sym.lambdify(w,Fw_fase,
                            modules=equivalentes)
else:
    F_w_fase = lambda w: Fw_fase + 0*w

Fwi_magn = F_w_magn(wi) # evalua wi
Fwi_fase = F_w_fase(wi) 

# F(w) real e imaginaria
fig_Fw, graf_Fwi = plt.subplots(2,1)

graf_Fwi[0].plot(wi,np.real(Fwi),label='Re(F(w))',
                 color='orange')
graf_Fwi[0].legend()
graf_Fwi[0].set_ylabel('Re (F(w)) ')
graf_Fwi[0].grid()

graf_Fwi[1].plot(wi,np.imag(Fwi),label='Imag(F(w))',
                 color='brown')
graf_Fwi[1].legend()
graf_Fwi[1].set_xlabel('w')
graf_Fwi[1].set_ylabel('imag(F(w))')
graf_Fwi[1].grid()

plt.suptitle(r'F(w) = $'+str(sym.latex(Fw))+'$')
# plt.show()

# F(w) magnitud y fase
fig_Fwmg, graf_Fw = plt.subplots(2,1)

graf_Fw[0].plot(wi,Fwi_magn,label='F(w)_magn',
                color='orange')
Fwa0 = F_w_magn(0)
Fwa1 = F_w_magn(-a_k)
Fwa2 = F_w_magn(a_k)
if Fw_magn.has(w): # no es constante
    graf_Fw[0].stem([-a_k,a_k],[Fwa1,Fwa2],
                basefmt=' ',linefmt ='--')
    etiqueta1 = '('+str(a_k)+','+str(round(Fwa2,4))+')'
    graf_Fw[0].annotate(etiqueta1, xy=(a_k,Fwa2))
etiqueta0 = '('+str(0)+','+str(round(Fwa0,4))+')'
graf_Fw[0].scatter(0,Fwa0)
graf_Fw[0].annotate(etiqueta0, xy=(0,Fwa0))
graf_Fw[0].legend()
graf_Fw[0].set_ylabel('F(w) magnitud ')
graf_Fw[0].grid()

graf_Fw[1].plot(wi,Fwi_fase,label='F(w)_fase',
                 color='brown')
graf_Fw[1].axhline(np.pi/2,linestyle ='--')
graf_Fw[1].axhline(-np.pi/2,linestyle ='--')
if Fw_magn.has(w): # no es constante
    Fwa1f = F_w_fase(-a_k)
    Fwa2f = F_w_fase(a_k)
    graf_Fw[1].stem([-a_k,a_k],[Fwa1f,Fwa2f],
                basefmt=' ',linefmt ='--')
    etiqueta3 = '('+str(a_k)+','+str(round(Fwa2f,4))+')'
    graf_Fw[1].annotate(etiqueta3, xy=(a_k,Fwa2f))
graf_Fw[1].legend()
graf_Fw[1].set_xlabel('w')
graf_Fw[1].set_ylabel('F(w) fase')
graf_Fw[1].grid()

plt.suptitle(r'F(w) = $'+str(sym.latex(Fw))+'$')

plt.show()

[ exponencial decreciente ] [ exponencial decreciente con |t| ] [ rectangular ]  [ Pulso unitario ]


Ejercicio 2. exponencial decreciente con |t|, función par

Referencia: Oppenheim Ejercicio 4.2 p291, Lathi ejemplo 7. p685, Hsu 5.21 p248

Considere la señal contínua exponencial decreciente, desarrolle la transformada de Fourier

x(t) =e^{-a|t|} \text{ ; } a \gt 0

X(\omega) = \int_{-\infty}^{\infty}e^{-a|t|} e^{-j \omega t} \delta t X(\omega) = \int_{-\infty}^{0}e^{at} e^{-j \omega t} \delta t + \int_{0}^{\infty}e^{-at} e^{-j \omega t} \delta t = \int_{-\infty}^{0}e^{at-j \omega t} \delta t + \int_{0}^{\infty}e^{-at-j \omega t} \delta t = \frac{1}{at-j \omega} e^{(a-j \omega) t}\Big|_{-\infty}^{0} - \frac{1}{a+j\omega}e^{-(at+j \omega) t} \Big| _{0}^{\infty} = \Big[ \frac{1}{at-j \omega} e^{(a-j \omega) (0)} - \frac{1}{at-j \omega t} e^{(a-j \omega) (-\infty)} \Big] + - \Big[ \frac{1}{a+j\omega}e^{-(at+j \omega) (\infty)} - \frac{1}{a+j\omega}e^{-(at+j \omega)(0)}\Big] = \frac{1}{a-j\omega} +\frac{1}{a+j\omega} = \frac{a+j\omega+ a-j\omega}{(a-j\omega)(a+j\omega)} X(\omega) = \frac{2a}{(a^2+\omega^2)}

Para desarrollar el ejercicio con el algoritmo, la entrada de señal se expresaría en el algoritmo como:

ft = sym.exp(-a*sym.Abs(t)) # Ej 4.2 Oppenheim

el resultado a comparar con el desarrollo analítico es:

 expresion f(t):
 -a*|t|
e      

 expresion F(w):
  2*a  
-------
 2    2
a  + w 

 |F(w)|:
  2*a  
-------
 2    2
a  + w 

 F(w)_fase:
0

Grafica de F(w) magnitud y fase

El algoritmo es el mismo que el ejercicio 1, modificando el bloque de ingreso para el problema

[ exponencial decreciente ] [ exponencial decreciente con |t| ] [ rectangular ] [ Pulso unitario ]


Ejercicio 3. Rectangular centrada en origen

Referencia: Oppenheim Ejercicio 4.4 p293, Lathi ejemplo 7.2 p689, Hsu 5.19 p247

Considere la señal pulso rectangular o pulso compuerta (gate), desarrolle la transformada de Fourier

x(t) =\begin{cases}1 && |t|<T_1, \\ 0 && |t|>T_1\end{cases} X(\omega) = \int_{-T_1}^{T_1} e^{-j\omega t} \delta t = -\frac{1}{j \omega} e^{-j\omega t}\Big|_{-T_1}^{T_1} = -\frac{1}{j \omega} \Big[ e^{-j\omega (T_1)} - e^{-j\omega (-T_1)}\Big] = -\frac{1}{j \omega} \Big[ e^{-T_1 j\omega } - e^{jT_1\omega } \Big] = \frac{1}{j \omega} e^{T_1 j\omega} - \frac{1}{j \omega} e^{-T_1 j\omega }

en este punto es conveniente usar la forma trigonometrica de un exponencial con exponente complejo

= \frac{1}{j \omega} (\cos (T_1\omega)+j \sin(T_1\omega)) - \frac{1}{j \omega} (cos(T_1 \omega) -jsin(T_1\omega)) = \frac{1}{j \omega}\cos (T_1\omega)+j\frac{1}{j \omega} \sin(T_1\omega) - \frac{1}{j \omega} cos(T_1 \omega) +j\frac{1}{j \omega} sin(T_1\omega)) X(\omega) = 2\frac{\sin(T_1\omega)}{\omega}

Para desarrollar el ejercicio con el algoritmo, la señal se expresaría como:

ft = sym.Heaviside(t+T)-sym.Heaviside(t-T) # Ej 4.4 Oppenheim

obteniendo el siguiente resultado

 expresion f(t):
-Heaviside(-T + t) + Heaviside(T + t)

 expresion F(w):
2*sin(T*w)
----------
    w     

 |F(w)|:
  |sin(T*w)|
2*|--------|
  |   w    |

 F(w)_fase:
0

con la siguiente gráfica f(T)

gráfica de F(w) parte real e imaginaria

El algoritmo es el mismo que el ejercicio 1, modificando el bloque de ingreso para el problema

[ exponencial decreciente ] [ exponencial decreciente con |t| ] [ rectangular ] [ Pulso unitario ]


Ejercicio 4. Pulso unitario

Referencia: Oppenheim Ejercicio 4.3 p292, Lathi ejemplo 7.3 p693, Hsu 5.1 p218

Considere la señal pulso unitario, desarrolle la transformada de Fourier

x(t) = \delta (t)

X(j\omega) = \int_{-\infty}^{\infty} \delta (t) e^{-j \omega t} \delta t = 1

Es decir la transformada de Fourier tiene componentes de todas las frecuencias.

El algoritmo entrega el siguiente resultado:

 expresion f(t):
DiracDelta(t)

 expresion F(w):
1

 |F(w)|:
1

 F(w)_fase:
0

El pulso unitario se define en Sympy como:

ft = sym.DiracDelta(t) # Ej 4.3 Oppenheim

En la parte gráfica, el pulsos unitarios se grafican con plt.stem(0,1), no requiriendo el resto de las graficas para observar el resultado.

El algoritmo es el mismo que el ejercicio 1, modificando el bloque de ingreso para el problema

Tarea: Realizar otros ejercicios con exponenciales para comprobar la operación del algoritmo.

[ exponencial decreciente ] [ exponencial decreciente con |t| ] [ rectangular ] [ Pulso unitario ]

5.1.1 Series de Fourier – Señales periódicas ejemplos con Python

Ejemplos para interpretar resultados de series de Fourier de señales periódicas contínuas a partir del algoritmo de la sección anterior.

[exponencial periódica] [triangular periódica] [rectangular periódica]


Ejemplo 1. Señal exponencial periódica con n términos de Fourier

Referencia: Lathi Ej 6.1 p599

Encuentre serie de Fourier para la señal x(t) y muestre el espectro de amplitud y fase de x(t)

x(t) = e^{-\frac{t}{2}} \text{ , }{0 \leq t \leq \pi}

Serie Fourier Ej01 f(t) Periodica

siendo el periodo T0 = π y la frecuencia fundamental f0 = 1/T0  = 1/π Hz,

\omega_0 = \frac{a \pi}{T_0} = 2 \text{ rad/s}

1.1 Desarrollo analítico

x(t) = a_0 +\sum_{n=1}^{\infty} a_n \cos (2 n t) +b_n \sin (2 n t) a_0 = \frac{1}{\pi} \int_{0}^{\pi} e^{-t/2} \delta t = 0.504 a_n = \frac{2}{\pi} \int_{0}^{\pi} e^{-t/2} \cos(2 n t) \delta t = 0.504 \Big( \frac{2}{1+16 n^2}\Big) b_n = \frac{2}{\pi} \int_{0}^{\pi} e^{-t/2} \sin (2 n t) \delta t = 0.504\Big( \frac{8n}{1+16 n^2}\Big) x(t) = 0.504\Big[ 1 +\sum_{n=1}^{\infty} \Big( \frac{2}{1+16 n^2} \Big)(\cos (2 n t) + 4n \sin (2 n t))\Big]

Serie de Fourier en la forma exponencial,

C_0 = a_0 = 0.504 C_n = \sqrt{a_n^2 + b_n^2} = 0.504\sqrt{ \Big(\frac{2}{1+16 n^2}\Big) ^2 + \Big( \frac{8n}{1+16 n^2} \Big) ^2 } C_n = 0.504\Big( \frac{8}{\sqrt{1+16 n^2}} \Big) \theta_n = \tan ^{-1} \Big( \frac{-b_n}{a_n}\Big) = \tan^{-1} (-4n) = -\tan^{-1}(4n)

con lo que se puede expresar x(t) como,

 serie de Fourier fs(t), n términos no cero: 
+ 0.504279523791290
+ 0.0593270027989753*cos(2*t)
+ 0.015516293039732*cos(4*t)
+ 0.00695557963850055*cos(6*t)
+ 0.00392435427074934*cos(8*t)
+ 0.00251510984434559*cos(10*t)
+ 0.00174793595768211*cos(12*t)
+ 0.0012847885956466*cos(14*t)
+ 0.00098396004642203*cos(16*t)
+ 0.000777609134604919*cos(18*t)
+ 0.000629955682437589*cos(20*t)
+ 0.237308011195901*sin(2*t)
+ 0.124130344317856*sin(4*t)
+ 0.0834669556620066*sin(6*t)
+ 0.0627896683319894*sin(8*t)
+ 0.0503021968869117*sin(10*t)
+ 0.0419504629843708*sin(12*t)
+ 0.0359740806781048*sin(14*t)
+ 0.0314867214855049*sin(16*t)
+ 0.0279939288457771*sin(18*t)
+ 0.0251982272975036*sin(20*t)

 coeficientes: 
['k_i', 'ak', 'bk', 'ck', 'cfs']
0 ak: 0.50427952379129 bk: 0
   ck: 0.50427952379129 cfs: 0
1 ak: 0.0593270027989753 bk: 0.2373080111959012
   ck: 0.24461149899148976 cfs: -1.3258176636680326
2 ak: 0.015516293039732001 bk: 0.12413034431785602
   ck: 0.1250963537844502 cfs: -1.4464413322481353
3 ak: 0.006955579638500553 bk: 0.08346695566200663
   ck: 0.08375627006732632 cfs: -1.4876550949064553
4 ak: 0.003924354270749339 bk: 0.06278966833198943
   ck: 0.06291218487450254 cfs: -1.5083775167989393
5 ak: 0.0025151098443455863 bk: 0.05030219688691173
   ck: 0.05036503538347567 cfs: -1.5208379310729538
6 ak: 0.0017479359576821145 bk: 0.04195046298437075
   ck: 0.04198686252526162 cfs: -1.5291537476963082
7 ak: 0.001284788595646599 bk: 0.03597408067810477
   ck: 0.035997016020363336 cfs: -1.5350972141155728
8 ak: 0.0009839600464220293 bk: 0.031486721485504944
   ck: 0.031502092109552245 cfs: -1.5395564933646284
9 ak: 0.0007776091346049191 bk: 0.02799392884577709
   ck: 0.02800472689009217 cfs: -1.5430256902014756
10 ak: 0.000629955682437589 bk: 0.025198227297503564
   ck: 0.025206100513536184 cfs: -1.5458015331759765
>>>

Instrucciones en Python:

# Serie de Fourier, espectro con n coeficientes
# Lathi ej 6.1 p599
import numpy as np
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True,)

T0 = sym.pi ; t0 = 0 # periodo ; t_inicio
ft = sym.exp(-t/2)
n = 11  # número de coeficientes

t_a = t0   # intervalo de t =[t_a,t_b]
t_b = t0 + T0

# PROCEDIMIENTO
serieF = sym.fourier_series(ft,(t,t0,t0+T0))
serieFk = serieF.truncate(n)

def fourier_series_coef(serieF,n,T0, casicero=1e-10):
    ''' coeficientes de serie de Fourier
        ak,bk,ck_mag, ck_fase
    '''
    w0 = 2*sym.pi/T0
    ak = [float(serieF.a0.evalf())]
    bk = [0]; k_i=[0]
    ak_coef = serieF.an
    bk_coef = serieF.bn
    ck_mag  = [ak[0]] ; ck_fase = [0]
    for i in range(1,n,1):
        ak_valor = ak_coef.coeff(i).subs(t,0)
        ak_valor = float(ak_valor.evalf())
        ak.append(ak_valor)
        bk_term = bk_coef.coeff(i).evalf()
        bk_valor = 0
        term_mul = sym.Mul.make_args(bk_term)
        for term_k in term_mul:
            if not(term_k.has(sym.sin)):
                bk_valor = float(term_k)
            else: # sin(2*w0*t)
                ki = term_k.args[0].subs(t,1)/w0
        bk.append(bk_valor)
        k_i.append(i)
        
        # magnitud y fase
        ak_signo = 1 ; bk_signo = 1
        if abs(ak_valor)>casicero:
            ak_signo = np.sign(ak_valor)
        if abs(bk_valor)>casicero:
            bk_signo = np.sign(bk_valor)
        signo_ck = ak_signo*bk_signo
        ck_mvalor = signo_ck*np.sqrt(ak_valor**2 + bk_valor**2)
        ck_mag.append(ck_mvalor)
        pendiente = np.nan
        if (abs(ak_valor)>=casicero):
            pendiente = -bk_valor/ak_valor
        ck_fvalor = np.arctan(pendiente)
        ck_fase.append(ck_fvalor)
    coef_fourier = {'k_i': k_i,'ak': ak,'bk': bk,
                    'ck_mag' : ck_mag,'ck_fase': ck_fase}
    return (coef_fourier)

coef_fourier = fourier_series_coef(serieF,n,T0)

# SALIDA
#print('\n serie de Fourier fs(t): ')
#sym.pprint(serieF)

print('\n serie de Fourier fs(t), n términos no cero: ')
term_suma = sym.Add.make_args(serieFk)
for term_k in term_suma:
    operador = '+'
    factor_mul = sym.Mul.make_args(term_k)
    for factor_k in factor_mul:
        if not(factor_k.has(t)):
            if sym.sign(factor_k)<0:
                operador = ''
    print(operador,term_k.evalf())

print('\n coeficientes: ')
m = len(coef_fourier['k_i'])
print(list(coef_fourier.keys()))
for i in range(0,m,1):
    print(str(coef_fourier['k_i'][i]),
          'ak:',str(coef_fourier['ak'][i]),
          'bk:',str(coef_fourier['bk'][i]) )
    print('  ','ck_mag:',str(coef_fourier['ck_mag'][i]),
          'ck_fase:',str(coef_fourier['ck_fase'][i]) )

gráfica de la serie de Fourier comparada con f(t) en un periodo,

Serie Fourier Ej01 SerieFt

Espectro de frecuencias,

Serie Fourier Ej01 Espectro k
Graficas en Python

Complementarias para realizar las gráficas, en un periodo, en varios periodos con periodo central destacado y espectro de frecuencias

# GRAFICA ----------------
import matplotlib.pyplot as plt
import telg1001 as fcnm
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
# Evaluación para gráfica
muestras  = 101
t_a = float(t_a) ; t_b = float(t_b)
figura_ft = fcnm.graficar_ft(ft,t_a,t_b,muestras)

ti = np.linspace(t_a,t_b,muestras)
# convierte a sympy una constante
serieF_k = sym.sympify(serieFk,t) 
if serieFk.has(t): # no es constante
    serieF_t = sym.lambdify(t,serieF_k,modules=equivalentes)
else:
    serieF_t = lambda t: serieF_k + 0*t
SerieFi = serieF_t(ti) # evalua ti

etiqueta = 'coeficientes = '+ str(n)
figura_ft.axes[0].plot(ti,SerieFi,'orange',
                    label = etiqueta)
figura_ft.axes[0].legend()
ft_etq = ''
if not(ft.has(sym.Piecewise)):
    ft_etq = '$ = '+str(sym.latex(ft)) +'$'
etiq_1 = 'f(t) '+ft_etq+' ; $T_0='+str(sym.latex(T0))+'$'
figura_ft.axes[0].set_title(r'Serie de Fourier '+etiq_1)
#plt.show()

def graficar_ft_periodoT0(ft,t0,T0,muestras=51,
                          n_periodos=4,f_nombre='f'):
    ''' grafica f(t) en intervalo[t0,t0+T0] para n_periodos
    '''
    # convierte a sympy una constante
    ft = sym.sympify(ft,t) 
    if ft.has(t): # no es constante
        f_t = sym.lambdify(t,ft,modules=equivalentes)
    else:
        f_t = lambda t: ft + 0*t
    # intervalo de n_periodos
    ti = np.linspace(float(t0-T0*n_periodos/2),
                     float(t0+T0*n_periodos/2),
                     n_periodos*muestras)
    fk = np.zeros(n_periodos*muestras)
    # ajuste de intervalo por periodos
    ti_T0 = (ti-t0)%float(T0)+t0
    fi = f_t(ti_T0)
    
    # intervalo de UN periodo
    ti0 = np.linspace(float(t0),float(t0+T0),muestras)
    fi0 = f_t(ti0)
    
    fig_fT0, graf_fT0 = plt.subplots()
    graf_fT0.plot(ti,fi,label=f_nombre+'(t)',
                 color= 'blue',linestyle='dashed')
    graf_fT0.plot(ti0,fi0,label=f_nombre+'(t)',color= 'blue')
    graf_fT0.axvline(t0, color ='red')
    graf_fT0.axvline(t0+T0, color ='red')
    graf_fT0.set_xlabel('t')
    graf_fT0.set_ylabel(f_nombre+'(t)')
    graf_fT0.legend()
    graf_fT0.grid()
    ft_etq = '' ; ft_ltx =''
    if not(ft.has(sym.Piecewise)):
        ft_etq = '$ = '+str(sym.latex(ft)) +'$'
    etiq_1 = r''+f_nombre+'(t) '+ft_etq+' ; $T_0='+str(sym.latex(T0))+'$'
    graf_fT0.set_title(etiq_1)
    return(fig_fT0)

def graficar_w_espectro(coef_fourier,T0,f_nombre='f'):
    ''' coef_fourier es diccionario con entradas
        ['k_i','ck_mag','ck_fase'] indice, Ck_magnitud, Ck_fase
    '''
    # espectro de frecuencia
    k_i = coef_fourier['k_i']
    ck  = coef_fourier['ck_mag']
    cfs = coef_fourier['ck_fase']

    # grafica de espectro de frecuencia
    fig_espectro_w, graf_spctr = plt.subplots(2,1)
    graf_spctr[0].stem(k_i,ck,label='Ck_magnitud')
    graf_spctr[0].set_ylabel('|Ck|')
    graf_spctr[0].legend()
    graf_spctr[0].grid()
    ft_etq = '' ; ft_ltx =''
    if not(ft.has(sym.Piecewise)):
        ft_etq = '$ = '+str(sym.latex(ft)) +'$'
    etiq_2 = ft_etq+' ; $T_0='+str(sym.latex(T0))+'$'
    etiq_1 = r'Espectro frecuencia '+f_nombre+'(t) '
    graf_spctr[0].set_title(etiq_1+etiq_2)
    graf_spctr[1].stem(k_i,cfs,label='Ck_fase')
    graf_spctr[1].legend()
    graf_spctr[1].set_ylabel('Ck_fase')
    graf_spctr[1].set_xlabel('k')
    graf_spctr[1].grid()
    return(fig_espectro_w)

fig_fT0 = graficar_ft_periodoT0(ft,t0,T0,muestras,f_nombre='x')
fig_espectro_w = graficar_w_espectro(coef_fourier,T0,f_nombre='x')
plt.show()

[exponencial periódica] [triangular periódica] [rectangular periódica]


Ejemplo 2. Señal triangular periódica

Referencia: Lathi Ej 6.2 p602

Encuentre serie de Fourier para la señal x(t) y muestre el espectro de amplitud y fase de x(t)

x(t) = \begin{cases} 4t & t< -\frac{1}{4} \\ -4t+4 & t\leq \frac{3}{4} \end{cases}

Serie Fourier Ej02 f(t) Periodica

en el algoritmo el bloque de ingreso se escribe el periodo, intervalo y la función como:

T0 = 2 ; t0 = -T0/4 # periodo ; t_inicio
A  = 2
ft = sym.Piecewise((0, t<(-T0/4)),
                   ((4*A/T0)*t, t<(T0/4)),
                   (-(4*A/T0)*t+2*A, t<=(3*T0/4)),
                   (0, True))
n = 9 # número de coeficientes

Serie de Fourier en forma trigonométrica:

 serie de Fourier fs(t), n términos no cero: 
+ 1.6211389382774*sin(pi*t)
+ 0.0200140609663877*sin(9*pi*t)
+ 0.00560947729507752*sin(17*pi*t)
+ 0.0648455575310962*sin(5*pi*t)
+ 0.009592538096316*sin(13*pi*t)
 -0.0133978424651025*sin(11*pi*t)
 -0.180126548697489*sin(3*pi*t)
 -0.00720506194789958*sin(15*pi*t)
 -0.0330844681281103*sin(7*pi*t)

 coeficientes: 
['k_i', 'ak', 'bk', 'ck', 'cfs']
0 ak: 0.0 bk: 0
   ck: 0.0 cfs: 0
1 ak: 0.0 bk: 1.6211389382774044
   ck: 1.6211389382774044 cfs: nan
2 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan
3 ak: 0.0 bk: -0.18012654869748937
   ck: 0.18012654869748937 cfs: nan
4 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan
5 ak: 0.0 bk: 0.06484555753109618
   ck: 0.06484555753109618 cfs: nan
6 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan
7 ak: 0.0 bk: -0.03308446812811029
   ck: 0.03308446812811029 cfs: nan
8 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan

Comparando la serie de Fourier con f(t)

Serie Fourier Ej02 Serie f(t)

Espectro de frecuencias de f(t)

Serie Fourier Ej02 Espectro

[exponencial periódica] [triangular periódica] [rectangular periódica]


Ejemplo 3. Señal rectangular periódica

Referencia: Lathi Ej 6.4 p605, Oppenheim Ej 4.5 p294

Encuentre serie de Fourier para la señal x(t) y muestre el espectro de amplitud y fase de x(t)

x(t) = \begin{cases} 0 & t<-\frac{\pi}{2} \\ 1 & -\frac{\pi}{2} \leq t \lt \frac{\pi}{2} \\ 0 & t>\frac{\pi}{2}\end{cases}
T0 = 2*sym.pi ; t0 = -T0/2 # periodo ; t_inicio
ft = sym.Piecewise((0, t<-(T0/4)),
                   (1, t<(T0/4)),
                   (0, t<(T0/2)),
                   (0, True))
n = 9 # numero de coeficientes

la gráfica de f(t) dentro de un periodo:

Serie Fourier Ej03 f(t) Periodica

la serie de Fourier en forma trigonométrica

 serie de Fourier fs(t), n términos no cero: 
+ 0.500000000000000
+ 0.636619772367581*cos(0.318309886183791*pi*t)
+ 0.0489707517205832*cos(4.13802852038928*pi*t)
+ 0.0707355302630646*cos(2.86478897565412*pi*t)
+ 0.127323954473516*cos(1.59154943091895*pi*t)
 -0.0909456817667973*cos(2.22816920328654*pi*t)
 -0.0578745247606892*cos(3.5014087480217*pi*t)
 -0.0424413181578388*cos(4.77464829275686*pi*t)
 -0.212206590789194*cos(0.954929658551372*pi*t)

 coeficientes: 
['k_i', 'ak', 'bk', 'ck', 'cfs']
0 ak: 0.5 bk: 0
   ck: 0.5 cfs: 0
1 ak: 0.6366197723675814 bk: 0.0
   ck: 0.6366197723675814 cfs: -0.0
2 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan
3 ak: -0.21220659078919377 bk: 0.0
   ck: -0.21220659078919377 cfs: 0.0
4 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan
5 ak: 0.12732395447351627 bk: 0.0
   ck: 0.12732395447351627 cfs: -0.0
6 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan
7 ak: -0.09094568176679733 bk: 0.0
   ck: -0.09094568176679733 cfs: 0.0
8 ak: 0.0 bk: 0.0
   ck: 0.0 cfs: nan
>>>

comparando la serie de Fourier con f(t):

Serie Fourier Ej03 Serie f(t)

el espectro de frecuencias,

Serie Fourier Ej03 Espectro

[exponencial periódica] [triangular periódica] [rectangular periódica]

5.1 Series de Fourier de señales periódicas con Sympy-Python

Referencia: Lathi 6.1 p593, Oppenheim 3.3 p186 Chapra 5ed 19.2 p546

La serie de Fourier aproxima una señal o función contínua mediante una serie infinita de sinusoides.

f(t) = a_0 + \sum_{k=1}^{\infty}[a_k \cos( k \omega_0 t) + b_k \sin(k \omega_0t)]

donde los coeficientes de la ecuación se calculan como:

a_k = \frac{2}{T} \int_0^T f(t) \cos(k \omega_0 t) \delta t b_k = \frac{2}{T} \int_0^T f(t) \sin(k \omega_0 t) \delta t

algoritmo: [ integral con Sympy ] [ instrucción Sympy ]


Ejemplo. Serie de Fourier de rectángulo(t) con Integral de Sympy

Referencia: Lathi 6.5 p617, Chapra 5ed Ejemplo 19.2 p548

Utilice la serie de Fourier continua para aproximar la función de onda cuadrada o rectangular. Para el cálculo del ejemplo se usará hasta 4 términos de la serie.

f(t) = \begin{cases} -1 && -T/2<t<-T/4 \\ 1 && -T/4<t<T/4 \\ -1 && T/4<t<T/2 \end{cases}

Serie Fourier Cuadrado 01


Coeficientes ak y bk

Para determinar las expresiones de los coeficientes, en Python se usa la librería Sympy que nos facilita el desarrollo las fórmulas simbólicas ak y bk.

Si requiere revisar el concepto se adjunta el enlace:

Fórmulas y funciones simbólicas con Python – Sympy

El desarrollo a papel y lápiz se deja como tarea.

Usando Python se obtiene los siguientes resultados:

 expresión ak:
/  /     /pi*k\            \                                                  
|2*|2*sin|----| - sin(pi*k)|                                                  
|  \     \ 2  /            /                                                  
<--------------------------- for And(k > -oo, k < oo, k != 0
|            pi*k                                                             
|                                                                             
\             0                                                     otherwise 


 expresión bk:
0

usando formato latex para la expresión, generado por Sympy se obtiene:

\begin{cases} \frac{2 \cdot \left(2 \sin{\left(\frac{\pi k}{2} \right)} - \sin{\left(\pi k \right)}\right)}{\pi k} & \text{for}\: \left(k > -\infty \vee k > 0\right) \wedge \left(k > -\infty \vee k < \infty\right) \wedge \left(k > 0 \vee k < 0\right) \wedge \left(k < 0 \vee k < \infty\right) \\0 & \text{otherwise} \end{cases} b_k = 0

Instrucciones en Python

# Serie de Fourier, con n coeficientes
# Ref.Chapra 5ed Ejemplo 19.2 p548/pdf572
import sympy as sym

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

T0  = 2*sym.pi ; t0 = -T0/2 # periodo ; t_inicio
ft = sym.Piecewise((-1,t <-T0/2),
                   (-1,t <-T0/4),
                   ( 1,t < T0/4),
                   (-1,t < T0/2),
                   (-1,True),)
n = 4 # número de coeficientes

# PROCEDIMIENTO
w0 = 2*sym.pi/T0
k  = sym.Symbol('k')

# Términos ak para coseno()
enintegral  = ft*sym.cos(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,t0,t0 + T0))
Fak = (2/T0)*yaintegrado
Fak = sym.simplify(Fak)

# Términos bk para seno()
enintegral = ft*sym.sin(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,t0,t0 + T0))
Fbk = (2/T0)*yaintegrado
Fbk = sym.simplify(Fbk)

# SALIDA
print(' expresión ak:')
sym.pprint(Fak)
print('\n ak formato latex')
print(sym.latex(Fak))

print('\n expresión bk:')
sym.pprint(Fbk)
print('\n bk formato latex')
print(sym.latex(Fbk))


Evaluación de coeficientes

Con las expresiones obtenidas en el bloque anterior, se evalúan los n coeficientes ak y bk, substituyendo en las expresiones los valores en cada índice i y se obtiene como resultado:

 coeficientes ak,bk: 
k_i : [0, 1, 2, 3]
ak : [0, 4/pi, 0, -4/(3*pi)]
bk : [0, 0, 0, 0]]

Instrucciones Python

Las instrucciones son adicionales al bloque anterior. La evaluación se mantiene usando las expresiones simbólicas usando la instrucción .subs()

# evalua los coeficientes
a0  = Fak.subs(k,0)/2
b0  = Fbk.subs(k,0)
ak = [a0] ; bk = [b0] ; k_i = [0]
i = 1
while not(i>=n):
    ak_valor = Fak.subs(k,i)
    bk_valor = Fbk.subs(k,i)
    ak.append(ak_valor)
    bk.append(bk_valor)
    k_i.append(i)
    i = i+1
print('\n coeficientes ak,bk: ')
for uncoef in coef_fourier:
    print(uncoef,':',coef_fourier[uncoef])

Expresión de la Serie de Fourier

Encontrados los coeficientes ak y bk, se los usa en la expresión principal

f(t) = a_0 + \sum_{k=1}^{\infty}[a_k \cos( k \omega_0 t) + b_k \sin(k \omega_0t)]

obteniendo la siguiente expresión para la serie de Fourier  como fs(t)

 serie de Fourier fs(t): 
4*cos(t)   4*cos(3*t)
-------- - ----------
   pi         3*pi  

Instrucciones en Python

Las instrucciones se añaden a continuación de los bloques anteriores,

# serie de Fourier
serieF = ak[0] + 0*t 
i = 1
while not(i>=n):
    serieF = serieF + ak[i]*sym.cos(i*w0*t)
    serieF = serieF + bk[i]*sym.sin(i*w0*t)
    i = i+1

print('\n serie de Fourier f(t): ')
sym.pprint(serieF)

Gráficas de f(t) y Serie de Fourier

Para comparar la función f(t) con la aproximación en series de Fourier, se usa el intervalo [a,b] con una cantidad de muestras usando la instrucción np.linspace() y guardando el resultado en ti.

Para evaluar los puntos ti en cada expresión, por simplicidad se convierte la expresión de su forma simbólica a numérica lambda, usando sym.lambdify()

Instrucciones en Python

Las instrucciones para graficar las series de Fourier van a continuación de las anteriores,

# Grafica ---------------------
import numpy as np
import matplotlib.pyplot as plt

# Evaluación para gráfica
f_t = sym.lambdify(t,ft)
serie_ft = sym.lambdify(t,serieF)
ti = np.linspace(float(t_a),float(t_b),muestras)
fi = f_t(ti)
serie_fi = serie_ft(ti) 

# Grafica de Serie de Fourier
plt.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
plt.plot(ti,serie_fi,label = etiqueta)
plt.xlabel('ti')
plt.legend()
plt.grid()
plt.title('Serie de Fourier f(t); T0='+str(T0))
plt.show()

Serie Fourier Cuadrado 01

Tarea: Realizar el ejercicio, aumentando el número de términos a 8

Espectro de Fourier de amplitud y fase

Referencia: Lathi 6.1 p599, Chapra 19.3 p551

Los espectros ofrecen información que no aparece en el dominio del tiempo. La gráfica de frecuencias ofrece una representación rápida de la estructura de armónicas, que son como las huellas dactilares que ayudan a caracterizar y entender la forma de una señal complicada.

La gráfica considera todas las magnitudes como positivas.  Posteriormente se puede observar  como en algunos textos que se incorpora el signo en la grafica de magnitud. En la siguiente sección se trata más éste detalle.

Fourier Cuadrado 01 freq

Para mostrar este espectro de frecuencias se incrementó el número de términos de la serie a n=11, para observar mejor la forma de la gráfica.

# espectro de frecuencias k
ak_i = np.array(ak,dtype=float)
bk_i = np.array(bk,dtype=float)

ck_mag = np.sqrt(ak_i**2 + bk_i**2)
ck_fase = np.arctan(-bk_i/ak_i)

revisa0 = (abs(ak_i)>=casicero)
pendiente = -bk_i[revisa0]/ak_i[revisa0]
ck_fase[revisa0] = np.arctan(pendiente)

coef_fourier['ck_mag']  = ck_mag
coef_fourier['ck_fase'] = ck_fase

print('\n coeficientes ak,bk,Ck_mag,Ck_fase: ')
for uncoef in coef_fourier:
    print(uncoef,':',coef_fourier[uncoef])

# grafica de espectro de frecuencia
plt.subplot(211)
plt.stem(k_i,ck_mag,label='|Ck|')
plt.ylabel('ck_mag')
plt.title('Espectro de frecuencia ; T0='+str(T0))
plt.legend()
plt.grid()
plt.subplot(212)
plt.stem(k_i,ck_fase,label='Ck_fase')
plt.legend()
plt.ylabel('ck_fase')
plt.xlabel('k_i')
plt.grid()
plt.show()

El algoritmo final como una integración de las partes presentadas se usa en la página siguiente con algunos ejemplos tradicionales de la transformada de Fourier.


Instrucciones Python  para Integral de la serie de Fourier con n coeficientes

Se integran todas las partes anteriores en un algoritmo

# Serie de Fourier, con n coeficientes
# Ref.Chapra 5ed Ejemplo 19.2 p548/pdf572
import sympy as sym
import numpy as np

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

T0  = 2*sym.pi ; t0 = -T0/2 # periodo ; t_inicio
ft = sym.Piecewise((-1,t <-T0/2),
                   (-1,t <-T0/4),
                   ( 1,t < T0/4),
                   (-1,t < T0/2),
                   (-1,True),)
n = 4 # número de coeficientes

t_a = t0   # intervalo de t =[t_a,t_b]
t_b = t0 + T0
muestras = 101 # 51 resolucion grafica
casicero = 1e-10

# PROCEDIMIENTO
w0 = 2*sym.pi/T0
k  = sym.Symbol('k')

# Términos ak para coseno()
enintegral  = ft*sym.cos(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,t0,t0 + T0))
Fak = (2/T0)*yaintegrado
Fak = sym.simplify(Fak)

# Términos bk para seno()
enintegral = ft*sym.sin(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,t0,t0 + T0))
Fbk = (2/T0)*yaintegrado
Fbk = sym.simplify(Fbk)

# evalua los coeficientes
a0  = Fak.subs(k,0)/2
b0  = Fbk.subs(k,0)
ak = [a0] ; bk = [b0] ; k_i = [0]
i = 1
while not(i>=n):
    ak_valor = Fak.subs(k,i)
    bk_valor = Fbk.subs(k,i)
    ak.append(ak_valor)
    bk.append(bk_valor)
    k_i.append(i)
    i = i+1
coef_fourier = {'k_i': k_i,
                'ak' : ak, 'bk': bk}

# serie de Fourier
serieF = ak[0] + 0*t 
i = 1
while not(i>=n):
    serieF = serieF + ak[i]*sym.cos(i*w0*t)
    serieF = serieF + bk[i]*sym.sin(i*w0*t)
    i = i+1

# espectro de frecuencias k
ak_i = np.array(ak,dtype=float)
bk_i = np.array(bk,dtype=float)

ck_mag = np.sqrt(ak_i**2 + bk_i**2)
ck_fase = np.arctan(-bk_i/ak_i)

revisa0 = (abs(ak_i)>=casicero)
pendiente = -bk_i[revisa0]/ak_i[revisa0]
ck_fase[revisa0] = np.arctan(pendiente)

coef_fourier['ck_mag']  = ck_mag
coef_fourier['ck_fase'] = ck_fase

# SALIDA
print(' expresión ak:')
sym.pprint(Fak)
print('\n ak formato latex')
print(sym.latex(Fak))

print('\n expresión bk:')
sym.pprint(Fbk)
print('\n bk formato latex')
print(sym.latex(Fbk))

print('\n serie de Fourier f(t): ')
sym.pprint(serieF)

print('\n coeficientes ak,bk,Ck_mag,Ck_fase: ')
for uncoef in coef_fourier:
    print(uncoef,':',coef_fourier[uncoef])

# Grafica ---------------------
import matplotlib.pyplot as plt

# Evaluación para gráfica
f_t = sym.lambdify(t,ft)
serie_ft = sym.lambdify(t,serieF)
ti = np.linspace(float(t_a),float(t_b),muestras)
fi = f_t(ti)
serie_fi = serie_ft(ti) 

# Grafica serie de Fourier
fig_serieF, graf_sF = plt.subplots()
graf_sF.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
graf_sF.plot(ti,serie_fi,label = etiqueta)
graf_sF.set_xlabel('ti')
graf_sF.legend()
graf_sF.grid()
graf_sF.set_title('Serie de Fourier f(t); T0='+str(T0))
# plt.show()

# grafica de espectro de frecuencia
fig_espectro, graf_sptr = plt.subplots(2,1)
graf_sptr[0].stem(k_i,ck_mag,label='|Ck|')
graf_sptr[0].set_ylabel('ck_mag')
graf_sptr[0].set_title('Espectro de frecuencia ; T0='+str(T0))
graf_sptr[0].legend()
graf_sptr[0].grid()

graf_sptr[1].stem(k_i,ck_fase,label='Ck_fase')
graf_sptr[1].legend()
graf_sptr[1].set_ylabel('ck_fase')
graf_sptr[1].set_xlabel('k_i')
graf_sptr[1].grid()
plt.show()

algoritmo: [ integral con Sympy ] [ instrucción Sympy ]



Instrucción Sympy para Serie de Fourier

Usando el ejercicio 1, se muestra el resultado obtenido usando la instrucción incorporada en Sympy.  La función truncate(n) aplicada a la serie, permite obtener los n términos no cero. Existen operaciones adicionales para desplazamiento y escala en x que evitan tener que realizar las operaciones nuevamente y optimizan el tiempo del algoritmo

 serie de Fourier f(t): 
4*cos(t)   4*cos(3*t)   4*cos(5*t)      
-------- - ---------- + ---------- + ...
   pi         3*pi         5*pi         

 serie de Fourier f(t), n términos: 
             /pi\                  /3*pi\
4*cos(t)*sinc|--|   4*cos(3*t)*sinc|----|
             \4 /                  \ 4  /
----------------- - ---------------------
        pi                   3*pi        

 serie de Fourier f(t), k términos no cero: 
4*cos(t)   4*cos(3*t)   4*cos(5*t)   4*cos(7*t)
-------- - ---------- + ---------- - ----------
   pi         3*pi         5*pi         7*pi   
>>> 

Instrucciones en Python

# Serie de Fourier, con n coeficientes
# Ref.Chapra 5ed Ejemplo 19.2 p548/pdf572
import sympy as sym

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

T0  = 2*sym.pi ; t0 = 0 # periodo ; t_inicio
ft = sym.Piecewise((-1,t <-T0/2),
                   (-1,t <-T0/4),
                   ( 1,t < T0/4),
                   (-1,t < T0/2),
                   (-1,True),)
n = 4 # número de coeficientes

t_a = t0   # intervalo de t =[t_a,t_b]
t_b = t0 + T0

# PROCEDIMIENTO
serieF = sym.fourier_series(ft,(t,t0,t0+T0))
serieFn = sym.expand(serieF.sigma_approximation(n))
serieFk = serieF.truncate(n)

# SALIDA
print('\n serie de Fourier f(t): ')
sym.pprint(serieF)

print('\n serie de Fourier f(t), n términos: ')
sym.pprint(serieFn)

print('\n serie de Fourier f(t), k términos no cero: ')
sym.pprint(serieFk)

Adicionalmente incorpora otras operaciones para el dominio de la frecuencia como desplazamiento y escalamiento en x

>>> sym.pprint(serieF)
4*cos(t)   4*cos(3*t)   4*cos(5*t)      
-------- - ---------- + ---------- + ...
   pi         3*pi         5*pi         

>>> sym.pprint(serieF.scalex(2))
4*cos(2*t)   4*cos(6*t)   4*cos(10*t)      
---------- - ---------- + ----------- + ...
    pi          3*pi          5*pi         

>>> sym.pprint(serieF.shiftx(2))
4*cos(t + 2)   4*cos(3*t + 6)   4*cos(5*t + 10)      
------------ - -------------- + --------------- + ...
     pi             3*pi              5*pi           
>>> 

algoritmo: [ integral con Sympy ] [ instrucción Sympy ]

4.4.2 LTI Laplace – Ejercicio para Y(s)=H(s)X(s) con constantes símbolo con Sympy-Python

Los ejercicios incluyen constantes tipo símbolo, que por simplicidad de la ecuación en la resolución se usan solo las operaciones básicas para la transformada Inversa de Laplace.

Ejemplo 1.  Y(s) dado H(s) y X(s) con una constante α

Referencia: Oppenheim ejemplo 9.37 p718 pdf749, semejante al ejercicio sistema-modelo LTI Lathi ejemplo 1.16 p111. Oppenheim problema 2.61c p164 pdf195

Suponga un sistema LTI causal que se describe mediante la ecuación diferencial:

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

que en el dominio ‘s’ se escribe como

s^2 Y(s) + 3 s Y(s) + 2 Y(s) = X(s) (s^2 + 3 s + 2) Y(s) = X(s) H(s) = \frac{X(s)}{Y(s)} = \frac{1}{s^2 + 3 s + 2}

Que es semejante a lo realizado en el ejemplo 2 de H(s) Diagramas de bloques, dominio ‘s’ y ecuaciones diferenciales, de donde se tiene que:

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

Suponga que la señal de entrada es x(t) = α μ(t) , usando la tabla de transformadas de Laplace,  la entrada del sistema es,

X(s) = \frac{\alpha}{s} = \frac{\alpha}{s+0}

Siempre que las señales de entrada son cero para t<0, la convolución es la multiplicación de los términos. (no aplica para valores t<0)

Y(s) = H(s) X(s) = \Big[\frac{1}{(s+1)(s+2)} \Big]\Big[\frac{\alpha}{s+0}\Big]

usando fracciones parciales con Python para encontrar el equivalente, la expresión se puede separar en,

Y(s) = \frac{\alpha}{2}\frac{1}{s+0} - \alpha\frac{1}{s+1} +\frac{\alpha}{2}\frac{1}{s+2}

Luego, aplicando la transformada inversa de Laplace de 1/(s+a) que es una exponencial e-at μ(t), se tiene la respuesta en el dominio del tiempo.

y(t) = \frac{\alpha}{2}e^{0} u(t) - \alpha e^{-t}u(t) + \frac{\alpha}{2} e^{- 2t} u(t) y(t) = \alpha \Big [ \frac{1}{2} - e^{-t} + \frac{1}{2} e^{- 2t} \Big] u(t)

Instrucciones con Sympy-Python

Se crea la expresión H(s) y X(s) en forma simbólica usando las variables, s y a.

Dado que la expresión es tipo polinomio, y que no existen exponenciales, el desarrollo del algoritmo puede realizarse con operaciones básicas realizadas en las  secciones anteriores. Añadiendo la operación para encontrar la salida Y(s)

# Señal de salida Y(s)
Ys = Hs*Xs
Ys_parcial = Ys.apart(s)

Además de realizar la transformada inversa de Laplace en cada término de Y(s).  Con lo que se obtiene el resultado:

Ys:
       a        
----------------
  / 2          \
s*\s  + 3*s + 2/

 Ys = Hs*Xs :
    a         a      a 
--------- - ----- + ---
2*(s + 2)   s + 1   2*s

 y(t) :
                                         -2*t             
a*Heaviside(t)      -t                a*e    *Heaviside(t)
-------------- - a*e  *Heaviside(t) + --------------------
      2                                        2          
>>> sym.pprint(yt.simplify())
  / 2*t      t    \  -2*t             
a*\e    - 2*e  + 1/*e    *Heaviside(t)
--------------------------------------
                  2                   
>>>   

Instrucciones con Sympy-Python

# Y(s) = H(s)*X(s)
# Oppenheim ejemplo 9.37 p718/pdf746
import sympy as sym

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)
a = sym.Symbol('a', real=True)

Hs = 1/(s**2+3*s+2)
Xs = a/s

# PROCEDIMIENTO
# Señal de salida Y(s)
Ys = Hs*Xs
Ysp = Ys.apart(s)

yt = sym.inverse_laplace_transform(Ysp,s,t)

# SALIDA
print('Ys:')
sym.pprint(Ys)
print('\n Ys = Hs*Xs :')
sym.pprint(Ysp)
print('\n y(t) :')
sym.pprint(yt)

Ejemplo 2. Y(s) con condiciones iniciales tipo constante simbolo

Referencia: Oppenheim ejemplo 9.38 p719

Considere el sistema caracterizado por la ecuación diferencial con condiciones iniciales.

\frac{\delta^2}{\delta t^2}y(t) + 3 \frac{\delta}{\delta t}y(t) + 2y(t) = x(t) y(0^{-}) = \beta \text{, } y'(0^{-}) = \gamma

sea x(t) = αμ(t)

X(s) = \frac{\alpha}{s} = \frac{\alpha}{s+0}

Aplicando la transformada de Laplace en ambos lados de la ecuación con condiciones iniciales y aplicando la propiedad de diferenciación:

\frac{\delta^2}{\delta t^2}y(t) \Rightarrow s^2Y(s) - y(0^{-})s - y'(0^{-}) = s^2Y(s) - \beta s - \gamma \frac{\delta}{\delta t}y(t) \Rightarrow sY(s) - y(0^{-}) = sY(s) - \beta

resultados que se insertan en la ecuación original de derivadas de t:

[s^{2}Y(s) - \beta s - \gamma] + 3[sY(s) - \beta] + 2Y(s) = \frac{\alpha}{s}

que se reordenan para despejar Y(s),

[s^{2}+3s +2] Y(s) - \beta s - (\gamma + 3 \beta) = \frac{\alpha}{s} [s^{2}+3s +2] Y(s) = \frac{\alpha}{s}+ \beta s + (\gamma + 3 \beta) (s+1)(s+2)Y(s) = \frac{\alpha}{s}+ \beta s + (\gamma + 3 \beta) Y(s) = \frac{\alpha}{s(s+1)(s+2)}+ \frac{\beta s}{(s+1)(s+2)} + \frac{(\gamma + 3 \beta)}{(s+1)(s+2)} Y(s) = \frac{\beta (s+3)}{(s+1)(s+2)} + \frac{\gamma}{(s+1)(s+2)} + \frac{\alpha}{s(s+1)(s+2)}

donde Y(s) es la transformada unilateral de Laplace de y(t)

el último término, representa la respuesta del sistema LTI causal y la condición de reposo inicial, conocido también como la respuesta en estado cero.

Una interpretación análoga se aplica a los primeros dos términos del miembro derecho de la ecuación, que representa la respuesta del sistema cuando la entrada es cero (α=0) , conocida también como respuesta a entrada cero.

Observe que la respuesta a entrada cero es una función lineal de los valores de las condiciones iniciales. Demostrando que la respuesta es la superposición del estado cero y respuesta a entrada cero.

si α=2, β = 3 y γ=-5, al realizar la expansión en fracciones parciales encontramos que:

Y(s) = \frac{1}{s} - \frac{1}{s+1} + \frac{3}{s+2}

con la aplicación de las tablas de transformadas se obtiene:

y(t) = \Big[ 1 - e^{-t} + 3 e^{-2t} \Big] \mu (t) \text{, para } t\gt 0

Instrucciones con Sympy-Python

En este caso, se ingresan las propiedades de diferenciación para Laplace, aunque pueden crearse a partir de una fórmula mas general en la siguiente sección.

La ecuación se define como una expresión de (lado_izquierdo)=(lado_derecho) contruida a partir la ecuación del problema incluyendo las condiciones iniciales diferentes de cero.

Luego se busca la expresión solución para Ys, por simplicidad se la cambia a fracciones parciales, asi es más sencillo usar la tabla de transformadas de Laplace.

Para una sustituir los valores de a,b y g, se usa un diccionario {} con las parejas de {variable: valor}.

ecuacion:
   2                                 a
Y*s  + 3*Y*s + 2*Y - b*s - 3*b - g = -
                                     s

 Y(s): 
       2              
a + b*s  + 3*b*s + g*s
----------------------
     / 2          \   
   s*\s  + 3*s + 2/   

 Y(s) en Fracciones parciales
 a    a - 2*b - 2*g   a - 2*b - g
--- + ------------- - -----------
2*s     2*(s + 2)        s + 1   

sustituye valores de a, b y g
 Ys:
  3       1     1
----- - ----- + -
s + 2   s + 1   s

 y(t) con constantes simbolo 
/   2*t                                     t\  -2*t             
\a*e    + a - 2*b - 2*g + 2*(-a + 2*b + g)*e /*e    *Heaviside(t)
-----------------------------------------------------------------
                                2                                

sustituye valores de a, b y g

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

 y(t) expande expresion
                -t                   -2*t             
Heaviside(t) - e  *Heaviside(t) + 3*e    *Heaviside(t)
>>> 

Compare los resultados con lo expuesto en la parte teórica.

# Y(s) con condiciones iniciales como constantes
# Oppenheim ejemplo 9.38 p719
import numpy as np
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True)
s = sym.Symbol('s')
a,b,g = sym.symbols('a b g')
Y = sym.Symbol('Y')

# propiedad de diferenciacion
d2Y = (s**2)*Y - b*s - g
dY  = s*Y - b

Xs = a/s    # señal de entrada

# PROCEDIMIENTO
# ecuacion del sistema
Lizq = d2Y + 3*dY + 2*Y # Lado izquierdo
Lder = Xs
ecuacion = sym.Eq(Lizq, Lder)

# despeja la ecuacion para Y(s)
Ys = sym.solve(ecuacion,Y)[0]

Yspc = sym.apart(Ys,s) # incluye constante simbolo

Ysp = Yspc.subs({a:2,b:3,g:-5})

# y(t) con inversa transformada Laplace
ytc = sym.inverse_laplace_transform(Yspc,s,t)
yt = sym.inverse_laplace_transform(Ysp,s,t)

# SALIDA
print('ecuacion:')
sym.pprint(ecuacion)
print('\n Y(s): ')
sym.pprint(Ys)
print('\n Y(s) en Fracciones parciales')
sym.pprint(Yspc)
print('\nsustituye valores de a, b y g')
print(' Ys:')
sym.pprint(Ysp)
print('\n y(t) con constantes simbolo ')
sym.pprint(ytc)
print('\nsustituye valores de a, b y g')
print('\n y(t): ')
sym.pprint(yt)
print('\n y(t) expande expresion')
sym.pprint(sym.expand(yt))

4.4.1 LTI Laplace – Ejercicio resuelto para Y(s)=H(s)X(s) con Sympy-Python

Se presenta un ejercicio que integra los temas tratados en la Unidad 4 para el análisis del Sistema con Transformada de Laplace. El algoritmo se muestra como una referencia para adaptar en los ejercicios integradores usando las funciones resumidas en telg1001.py .

Referencia: Lathi Ejemplo 4.12 p361

Resuelva la ecuación diferencial de segundo orden

(D^2 + 5D + 6) y(t) = (D+1) x(t)

para las condiciones iniciales y(0-)=2 y y'(0-)=1, cuando se presenta la entrada

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

Desarrollo

La ecuación diferencial del sistema, interpretando los operadores D por derivadas se muestra a continuación.

\frac{\delta^2}{\delta t^2}y(t) + 5\frac{\delta}{\delta t}y(t) + 6y(t) = \frac{\delta^2}{\delta t^2}x(t) + x(t)

La ecuación se representa en el dominio s de las transformadas de Laplace para usarlas en el algoritmo:

s^2 Y(s) + 5s Y(s) + 6 Y(s) = s X(s) + X(s) (s^2 + 5s + 6)Y(s) = (s+1) X(s)

con lo que se escribe la función de transferencia H(s) que se ingresa para el análisis del algoritmo.

H(s) = \frac{X(s)}{Y(s)} = \frac{s+1}{s^2 + 5s + 6}

Para el caso de x(t) se usa la transformada de Laplace con el algoritmo para obtener X(s), mostrando que se puede realizar directamente para términos simples. Si la expresión es mas elaborada se usaría las funciones realizadas para f(t) con Sympy-Python

Con el algoritmo se obtienen los siguientes resultados, con lo que se puede analizar la estabilidad del sistema usando la gráfica de polos. No se observan polos en el lado derecho del plano RHP, por lo que el sistema es asintóticamente estable, BIBO-estable.

LTIC Laplace Algoritmo Integrador 01 polos

la gráfica de H(s) junto a los polos:

LTIC Laplace Algoritmo Integrador 01 polos Hs

la gráfica de h(t) como respuesta al impulso,

LTIC Laplace Algoritmo Integrador 01 ht

la gráfica de la señal de entrada x(t) y la respuesta y(t) es:

LTIC Laplace Algoritmo Integrador 01 xh_y

Las expresiones para Y(s), H(s), análisis de estabilidad que entrega el algoritmo son:

 H(s) = P(s)/Q(s):
  2       1  
----- - -----
s + 3   s + 2
 H(s) en factores:
     s + 1     
---------------
(s + 2)*(s + 3)

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

polosceros:
Q_polos : {-2: 1, -3: 1}
P_ceros : {-1: 1}

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

 X(s): 
  1  
-----
s + 4

Respuesta entrada cero ZIR H(s) y condiciones iniciales
term_cero : 2*s + 11
ZIR :
    5       7  
- ----- + -----
  s + 3   s + 2
yt_ZIR :
/   -2*t      -3*t\             
\7*e     - 5*e    /*Heaviside(t)

 ZSR respuesta estado cero:
ZSR :
      3         2         1    
- --------- + ----- - ---------
  2*(s + 4)   s + 3   2*(s + 2)
yt_ZSR :
/   -2*t                -4*t\             
|  e          -3*t   3*e    |             
|- ----- + 2*e     - -------|*Heaviside(t)
\    2                  2   /             

 Y(s)_total = ZIR + ZSR:
      3         3         13   
- --------- - ----- + ---------
  2*(s + 4)   s + 3   2*(s + 2)

 y(t)_total = ZIR + ZSR:
/    -2*t                -4*t\             
|13*e          -3*t   3*e    |             
|-------- - 3*e     - -------|*Heaviside(t)
\   2                    2   /             
>>>

Con los datos de H(s) y H(s) con fracciones parciales, se presentan dos formas de realizar el diagrama de bloques. El primer diagrama corresponde a H(s), mientras qu el segundo corresponde a la interpretación de las fracciones parciales.

Instrucciones Python

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

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

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

# H(s) respuesta impulso
Ps = s+1
Qs = s**2 + 5*s + 6
Hs = Ps/Qs
#Hs = 1+0*s cuando es constante

# X(s) Señal de entrada
xt = sym.exp(-4*t)*u

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

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

# PROCEDIMIENTO
Hs = fcnm.apart_s(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_s_parametros(Hs_fc)

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

estable = fcnm.estabilidad_asintotica_s(Q_polos)

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

# PROCEDIMIENTO Respuesta ZIR, ZSR
Xs = fcnm.laplace_transform_suma(xt)

# ZIR_s respuesta entrada cero de s
sol_ZIR = fcnm.respuesta_ZIR_s(Hs,cond_inicio)
ZIR = sol_ZIR['ZIR']
yt_ZIR = sol_ZIR['yt_ZIR']

# ZSR respuesta estado cero, Y(s) a entrada X(s)
sol_ZSR = fcnm.respuesta_ZSR_s(Hs,Xs)
ZSR = sol_ZSR['ZSR']
yt_ZSR = sol_ZSR['yt_ZSR']

# Respuesta total Y(s) y y(t)
Ys = ZIR + ZSR
Ys = fcnm.apart_s(Ys)
yt = yt_ZIR + yt_ZSR
lista_escalon = yt.atoms(sym.Heaviside)
yt = sym.collect(yt,lista_escalon)

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

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

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

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

print('\n X(s): ')
sym.pprint(Xs)
print('\nRespuesta entrada cero ZIR H(s) y condiciones iniciales')

if not(sol_ZIR == sym.nan): # existe resultado
    fcnm.print_resultado_dict(sol_ZIR)
else:
    print(' insuficientes condiciones iniciales')
    print(' revisar los valores de cond_inicio[]')

print('\n ZSR respuesta estado cero:')
fcnm.print_resultado_dict(sol_ZSR)

print('\n Y(s)_total = ZIR + ZSR:')
sym.pprint(Ys)
print('\n y(t)_total = ZIR + ZSR:')
sym.pprint(yt)

# Graficas polos, H(s), con polos h(t) --------
figura_s  = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,muestras,f_nombre='H')
figura_ht = fcnm.graficar_ft(ht,t_a,t_b,muestras,f_nombre='h')
# GRAFICAS y(t),x(t),h(t) ---------------------
figura_ft = fcnm.graficar_xh_y(xt,ht,yt,t_a,t_b,muestras)
plt.show()

4.4 LTI Laplace – Y(s)=ZIR+ZSR Respuesta entrada cero y estado cero con Sympy-Python

Referencia: Lathi Ejemplo 4.12 p361

Resuelva la ecuación diferencial de segundo orden usando respuesta a entrada cero y respuesta a estado cero.

(D^2 + 5D + 6) y(t) = (D+1) x(t)

para las condiciones iniciales y(0-)=2 y y'(0-)=1, ante la entrada

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

El ejercicio muestra nuevamente el concepto sobre la respuesta total del sistema tiene dos componentes:

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]
..


1. Respuesta a entrada cero y condiciones iniciales

Si el sistema bajo observación se no tiene señal de entrada X(s) = 0, la parte de estado cero no aporta componente de salida.

Respuesta
total
= respuesta a
entrada cero ZIR
+ respuesta a
estado cero ZSR

Sin embargo, si existen valores iniciales en los componentes del sistema (circuito con capacitores), se producirá una respuesta a entrada cero que básicamente responde a las condiciones iniciales al sistema.

Referencia: Lathi Ejemplo 4.12 p361

Un sistema LTI tiene la ecuación diferencial de segundo orden expresada con operadores D:

(D^2 + 5D + 6) y(t) = (D+1) x(t)

El sistema tiene las condiciones iniciales y(0-)=2 y y'(0-)=1.

1.1 Desarrollo Analítico para Respuesta a entrada cero y condiciones iniciales

La ecuación se escribe usando la transformada de Laplace:

(s^2 + 5s + 6) Y(s) = (s+1) X(s)

El diagrama básico de la expresión LTI, sin simplificar, se muestra en la figura. La sección de color naranja representa Q(s)  y la sección en verde representa P(s)

Para el análisis de entrada cero X(s)=0, la ecuación se simplifica a

(s^2 + 5s + 6) Y(s) = 0

Las condiciones iniciales se aplican usando la tabla de propiedades de Laplace, donde se muestra que los valores iniciales se añaden a los componentes de primera y segunda derivada:

\frac{\delta}{\delta t} y(t) \Longleftrightarrow sY(s) - y(0^{-}) = sY(s) -2 \frac{\delta^2}{\delta t^2} y(t) \Longleftrightarrow s^2 Y(s) - sy(0^{-}) - y'(0^{-}) = s^2 Y(s) -2s - 1

Se recuerda que para la entrada x(t) = 0, no se crean términos en el lado derecho de la expresión

[s^2 Y(s) -2s - 1] + 5[sY(s) -2] + 6Y(s) = 0

Agrupando los terminos Y(s), se tienen los componentes que no dependen de Y(s) que son los resultantes de las condiciones iniciales.

[s^2 + 5s + 6]Y(s) - (2s+11) = 0 [s^2 + 5s + 6]Y(s) = (2s+11) Y(s)= \frac{2s+11}{s^2 + 5s + 6}

El resultado Y(s) corresponde al término de entrada cero con condiciones iniciales aplicadas.

Aplicando un algoritmo en Python para añadir los componentes de condiciones iniciales, el resultado es:

Respuesta entrada cero ZIR H(s) y condiciones iniciales

 term_cero :
2*s + 11

 ZIR :
    5       7  
- ----- + -----
  s + 3   s + 2

 yt_ZIR :
/   -2*t      -3*t\             
\7*e     - 5*e    /*Heaviside(t)
>>>

1.2 Desarrollo con Sympy-Python

Se inicia con los resultados de algoritmos anteriores, donde se ingresan las expresiones de Q(s) y P(s) correspondientes a H(s). Las condiciones iniciales  se ingresan en un arreglo en orden de derivada descendente, por ejemplo para  el ejercicio [y'(0), y(0)] = [1., 2.]

Par generar las expresiones de condicion inicial, se genera el término de condición cero o term_cero, compuesto de los términos de orden tomados del vector inicial.

# Y_ZIR(s) Respuesta entrada cero con Transformada de Laplace
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# Ejemplo Lathi 4.12 p361
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

# H(s) respuesta impulso
Ps = s + 1
Qs = s**2 + 5*s + 6
Hs = Ps/Qs
#Hs = 6*sym.exp(-2*s)/s**2

# X(s) Señal de entrada
Xs = 0*s

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [1, 2]

# PROCEDIMIENTO
def respuesta_ZIR_s(Hs,cond_inicio):
    '''respuesta a entrada cero ZIR con H(s) y condiciones de inicio
       Si ZIR es sym.nan existen insuficientes condiciones
    '''
    polosceros = fcnm.busca_polosceros(Hs)
    Q_polos = polosceros['Q_polos']
    
    # coeficientes y grado de Q
    Q = 1+0*s
    for raiz in Q_polos:
        veces = Q_polos[raiz]
        Q = Q*((s-raiz)**veces)
    Q = sym.poly(Q,s)
    
    Q_coef  = Q.all_coeffs() # coeficientes Q
    Q_grado = Q.degree()  # grado polinomio Q
    
    term_cero = 0*s 
    if len(cond_inicio) == Q_grado:
        for j in range(0,Q_grado,1):
            term_orden = 0
            for i in range(j,Q_grado,1):
                term_orden = term_orden + cond_inicio[i]*(s**(i-j))
            term_cero = term_cero + term_orden*Q_coef[j]
        ZIR = term_cero/Q
        ZIR = fcnm.apart_s(ZIR)
        ZIR_Qs2 = fcnm.Q_cuad_s_parametros(ZIR)
    else:
        ZIR = sym.nan # insuficientes condiciones iniciales

    if not(ZIR==sym.nan):
        # entrada_cero en t
        yt_ZIR = sym.inverse_laplace_transform(ZIR,s,t)
        # simplifica log(exp()) ej: e**(-2t)/(s**2)
        if yt_ZIR.has(sym.log):
            yt_ZIR = sym.simplify(yt_ZIR,inverse=True)
        lista_escalon = yt_ZIR.atoms(sym.Heaviside)
        yt_ZIR = sym.expand(yt_ZIR,t) # terminos suma
        yt_ZIR = sym.collect(yt_ZIR,lista_escalon)
        
        sol_ZIR ={'term_cero' : term_cero}
        if len(ZIR_Qs2)>0: # añade si Q es cuadratico
            sol_ZIR['ZIR_Qs2'] = ZIR_Qs2
        sol_ZIR = sol_ZIR | {'ZIR'   : ZIR,
                             'yt_ZIR': yt_ZIR}
    else:
        sol_ZIR = sym.nan
    
    return(sol_ZIR)

# respuesta a entrada_cero ZIR en s
sol_ZIR = respuesta_ZIR_s(Hs,cond_inicio)

# SALIDA - polinomio
print('\nRespuesta entrada cero ZIR H(s) y condiciones iniciales')
if not(sol_ZIR == sym.nan): # existe resultado
    for k in sol_ZIR:
        print('\n',k,':')
        sym.pprint(sol_ZIR[k])
else:
    print(' insuficientes condiciones iniciales')
    print(' revisar los valores de cond_inicio[]')

El siguiente caso sería realizar el mismo sistema LTI CT, con respuesta Y(s) ante una extrada X(s), conocido también como estado cero.

Para mantenerlo simple, se consideraría condiciones iniciales en cero.

Finalmente se sumarían la entrada cero y estado cero para obtener la respuesta total.

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]

..


2. Respuesta a estado cero con x(t) de Sistema con ecuación diferencial de 2do orden

El ejercicio continua  el desarrollo del concepto sobre la respuesta total del sistema tiene dos componentes:

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

Referencia: Lathi Ejemplo 4.12 p361

Un sistema LTI tiene la ecuación diferencial de segundo orden expresada con operadores D:

(D^2 + 5D + 6) y(t) = (D+1) x(t)

(viene de la sección anterior…) si se somete a la entrada,

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

2.1. Desarrollo Analítico para estado cero

Ahora, todos los componentes tienen condición inicial cero y una señal de entrada x(t). El diagrama de bloques del sistema se obtuvo en la sección anterior de «entrada cero»:

En este caso se considera que ningún componente del sistema tiene valor inicial [0., 0.], interpretado por ejemplo como que todos los capacitores estan descargados, la salida corresponde solo a la parte conformada por x(t) y la respuesta al impulso h(t).

para la entrada x(t),

x(t) = e^{-4t} \mu(t) \Longleftrightarrow X(s) = \frac{1}{s+4}

Usando transformadas de Laplace, la ecuación del sistema se convierte en:

[s^2 Y(s)] + 5[sY(s)] + 6Y(s) = (s+1)X(s)

que al aplicar la parte de X(s),

[s^2 Y(s)] + 5[sY(s)] + 6Y(s) = (s+1)\frac{1}{s+4} (s^2 + 5s + 6) Y(s) = \frac{s+1}{s+4} Y(s) = \frac{s+1}{(s+4)(s^2 + 5s + 6)}

que es equivalente a y(t) = h(t)⊗x(t) como convolución en dominio de tiempo o Y(s)=H(s)X(s) en el dominio de s.

Para simplificar la respuesta de la inversa de transformada de Laplace, se separa en fracciones parciales. Al usar el algoritmo en Python se obtiene:

Y(s) = -\frac{1}{2}\frac{1}{s+2} +2\frac{1}{s+3} - \frac{3}{2} \frac{1}{s+4}

Tomando la referencia de la tabla de transformadas de Laplace, fila 5, se tiene:

y(t) = -\frac{1}{2} e^{-2t} \mu(t) + 2 e^{-3t} \mu(t) -\frac{3}{2} e^{-4t} \mu (t)

2.2. Instrucciones con Sympy-Python

Para este caso se añade la transformada de Laplace X(s) de x(t). La expresión μ(t) se representa como Heaviside en Sympy.

Las instrucciones son semejantes a las usadas en algoritmos anteriores, tal como Hs = Ps/Qs, considerando que respuesta a estado cero es ZSR = Hs*Xs

Se refleja el resultado obtenido en el desarrollo analítico,

 ZSR respuesta estado cero:

 ZSR :
      3         2         1    
- --------- + ----- - ---------
  2*(s + 4)   s + 3   2*(s + 2)

 yt_ZSR :
/   -2*t                -4*t\             
|  e          -3*t   3*e    |             
|- ----- + 2*e     - -------|*Heaviside(t)
\    2                  2   /             
>>> 

las instrucciones usadas son:

# Y_ZSR(s) Respuesta estado cero con Transformada Laplace
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# Ejemplo Lathi 4.12 p361
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

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

# X(s) Señal de entrada
Xs = 1/(s+4)

# PROCEDIMIENTO
def respuesta_ZSR_s(Hs,Xs):
    '''respuesta a estado cero ZSR con H(s) y X(s)
    '''
    ZSR = Hs*Xs
    ZSR = fcnm.apart_s(ZSR)
    ZSR_Qs2 = fcnm.Q_cuad_s_parametros(ZSR)

    # ZSR respuesta estado cero(t)
    yt_ZSR = sym.inverse_laplace_transform(ZSR,s,t)

    # simplifica log(exp()) ej: e**(-2t)/(s**2)
    if yt_ZSR.has(sym.log):
        yt_ZSR = sym.simplify(yt_ZSR,inverse=True)
    lista_escalon = yt_ZSR.atoms(sym.Heaviside)
    yt_ZSR = sym.expand(yt_ZSR,t) # terminos suma
    yt_ZSR = sym.collect(yt_ZSR,lista_escalon)

    sol_ZSR = {'ZSR' : ZSR}
    if len(ZSR_Qs2)>0:
        sol_ZSR['ZSR_Qs2'] = ZSR_Qs2
    sol_ZSR = sol_ZSR | {'yt_ZSR' : yt_ZSR}
    return(sol_ZSR)

sol_ZSR = respuesta_ZSR_s(Hs,Xs)

# SALIDA - polinomio
print('\n ZSR respuesta estado cero:')
for k in sol_ZSR:
    print('\n',k,':')
    sym.pprint(sol_ZSR[k])

El siguiente bloque consiste en hacer el ejercicio completo del sistema LTIC, considerando condiciones iniciales y una entrada x(t), que da la respuesta total del sistema.

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]

..


3. Y(s) respuesta total con condiciones iniciales y entrada x(t)

Finalmente se obtiene el resultado sumando los dos componentes desarrollados en  los numerales anteriores

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero

3.1 Desarrollo Analítico

Se incluye para mostrar que la respuesta total Y(s) = entrada_cero+estado_cero

Para resolver con transformadas de Laplace, la ecuación se expresa como :

\frac{\delta ^2}{\delta t^2} y(t) + 5 \frac{\delta}{\delta t} y(t) + 6 y(t) = \frac{\delta}{\delta t}x(t) +x(t)

haciendo y(t) ⇔Y(s) y usando la tabla de propiedades

\frac{\delta}{\delta t} y(t) \Longleftrightarrow sY(s) - y(0^{-}) = sY(s) -2 \frac{\delta^2}{\delta t^2} y(t) \Longleftrightarrow s^2 Y(s) - sy(0^{-}) - y'(0^{-}) = s^2 Y(s) -2s - 1

para la entrada x(t),

x(t) = e^{-4t} \mu(t) \Longleftrightarrow X(s) = \frac{1}{s+4} \frac{\delta}{\delta t}x(t) \Longleftrightarrow sX(s) - x(0^{-}) = s\frac{1}{s+4} -0 = \frac{s}{s+4}

con lo que la ecuación diferencial en transformada de Laplace se escribe,

[s^2 Y(s) -2s - 1] + 5[sY(s) -2] + 6Y(s) = \frac{s}{s+4}+\frac{1}{s+4}

Agrupando los terminos Y(s) y los términos de la derecha,

[s^2 + 5s + 6]Y(s) - (2s+11) = \frac{s+1}{s+4} [s^2 + 5s + 6]Y(s) = (2s+11) + \frac{s+1}{s+4} Y(s) = \frac{2s+11}{s^2 + 5s + 6} + \frac{s+1}{(s+4)(s^2 + 5s + 6)}

despejando Y(s) se puede observar que la expresión es la suma de,

Y(s) = \text{entrada cero} + \text{estado cero} Y(s) = \frac{(s+4)(2s+11)+s+1}{(s+4)(s^2 + 5s + 6)} Y(s)= \frac{2s^2+20s+45}{(s^2 + 5s + 6)(s+4)} Y(s)= \frac{2s^2+20s+45}{(s+2)(s+3)(s+4)}

expandiendo en fracciones parciales:

Y(s)= \Big(\frac{13}{2}\Big)\frac{1}{s+2} - 3\frac{1}{s+3} - \Big(\frac{3}{2}\Big)\frac{1}{s+4}

Finalmente, usando la transformada inversa de Laplace obtener la expresión en el dominio de t,

Y(t)= \Big[ \frac{13}{2}e^{-2t} - 3 e^{-3t} - \frac{3}{2}e^{-4t} \Big] \mu(t)

El ejercicio muestra la facilidad de resolver el sistema usando transformadas de Laplace.

Respuesta entrada cero ZIR H(s) y condiciones iniciales

 term_cero :
2*s + 11

 ZIR :
    5       7  
- ----- + -----
  s + 3   s + 2

 yt_ZIR :
/   -2*t      -3*t\             
\7*e     - 5*e    /*Heaviside(t)

 ZSR respuesta estado cero:

 ZSR :
      3         2         1    
- --------- + ----- - ---------
  2*(s + 4)   s + 3   2*(s + 2)

 yt_ZSR :
/   -2*t                -4*t\             
|  e          -3*t   3*e    |             
|- ----- + 2*e     - -------|*Heaviside(t)
\    2                  2   /             

 Y0(s) total = ZIR + ZSR:
      3         3         13   
- --------- - ----- + ---------
  2*(s + 4)   s + 3   2*(s + 2)

 y(t) total = ZIR + ZSR:
/    -2*t                -4*t\             
|13*e          -3*t   3*e    |             
|-------- - 3*e     - -------|*Heaviside(t)
\   2                    2   /             
>>> 

3.2 Instrucciones Sympy-Python

Las instrucciones para respuesta a entrada cero ZIR y respuesta estado cero ZSR se integran en un algoritmo.

# Y(s) Respuesta total con entrada cero y estado cero
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# Ejemplo Lathi 4.12 p361
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

# H(s) respuesta impulso
Ps = s + 1
Qs = s**2 + 5*s + 6
Hs = Ps/Qs
#Hs = 6*sym.exp(-2*s)/s**2

# X(s) Señal de entrada
Xs = 1/(s+4)

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [1, 2]

# PROCEDIMIENTO
# ZIR_s respuesta entrada cero de s
sol_ZIR = fcnm.respuesta_ZIR_s(Hs,cond_inicio)
ZIR = sol_ZIR['ZIR']
yt_ZIR = sol_ZIR['yt_ZIR']

# ZSR respuesta estado cero, Y(s) a entrada X(s)
sol_ZSR = fcnm.respuesta_ZSR_s(Hs,Xs)
ZSR = sol_ZSR['ZSR']
yt_ZSR = sol_ZSR['yt_ZSR']

# Respuesta total Y(s) y y(t)
Ys = ZIR + ZSR
Ys = fcnm.apart_s(Ys)
yt = yt_ZIR + yt_ZSR
lista_escalon = yt.atoms(sym.Heaviside)
yt = sym.collect(yt,lista_escalon)

# SALIDA - polinomio
print('\nRespuesta entrada cero ZIR H(s) y condiciones iniciales')
if not(sol_ZIR == sym.nan): # existe resultado
    for k in sol_ZIR:
        print('\n',k,':')
        sym.pprint(sol_ZIR[k])
else:
    print(' insuficientes condiciones iniciales')
    print(' revisar los valores de cond_inicio[]')

print('\n ZSR respuesta estado cero:')
for k in sol_ZSR:
    print('\n',k,':')
    sym.pprint(sol_ZSR[k])

print('\n Y0(s) total = ZIR + ZSR:')
sym.pprint(Ys)
print('\n y(t) total = ZIR + ZSR:')
sym.pprint(yt)

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]

4.3.3 H(s) – polos y Estabilidad del sistema con Sympy-Python

Referencia: Lathi 4.3-3 p371, Lathi 2.4-2 p198

En un sistema acotado (bounded) estable, sometido a una señal de entrada en la que se conoce que su salida siempre es la misma, se lo considera un sistema BIBO-estable. (Bounded input, Boundes output)

La función de transferencia H(s) de un sistema es una descripción externa con la forma:

H(s) = \frac{P(s)}{Q(s)}

El polinomio Q(s) representa una descripción interna, en la que si todos los polos de H(s) se encuentran en la parte izquierda del plano (left half plane, LHP), todos los términos h(t) son exponenciales decrecientes y h(t) es absolutamente integrable. En consecuencia el sistema es BIBO-estable.

Se asume que H(s) es una función en la que M ≤ N, siendo M el grado de P(s) y N e grado de Q(s), caso contrario el sistema es es BIBO-inestable.

Laplace estabilidad de un polo real animación

1. Un sistema LTI CT es asintóticamente estable si y solo si todos los polos de su función de transferencia H(s) se encuentran en el lado izquierdo del plano, pudiendo ser simples o repetidos.

2. Un sistema LTI CT es inestable si y solo si se dan una o ambas de las condiciones:

2.1 Existe al menos un polo de H(s) en el lado derecho del plano RHP.
2.2 Existen polos repetidos de H(s) en el eje imaginario.

3. Un sistema LTI CT es marginalmente estable si y solo si, no hay polos de H(s) en el lado derecho del plano RHP y algunos polos NO repetidos sobre en el eje imaginario (parte Real = 0).

Laplace Establilidad Dos Polos Imaginarios

 


polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]
..


Ejemplo 1. H(s) con polos no repetidos, RHP=0

Referencia: Lathi 2.14a p201

Revise la estabilidad BIBO de los siguientes sistemas LTI¨CT descritos por las ecuaciones:

(D+1)(D^2 + 4D +8) y(t) = (D-3) x(t)

Ejemplo 1 Desarrollo analítico

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s+1)(s^2 + 4s +8) Y(s) = (s-3) X(s) \frac{Y(s)}{X(s)} = \frac{s-3}{(s+1)(s^2 + 4s +8)}

El grado del polinomio Q es N=3 y es mayor que el grado del polinomio P es M=1.

La expresión de entrada se escribe en el algoritmo como:

Hs = ((s-3)/(s+1))*(1/(s**2+4*s+8))

Los polos se obtienen como las raíces del denominador Q, necesario para plantear las fracciones parciales y simplificar la expresión como la suma de componentes más simples.

polos: {-1: 1, -2 - 2*I: 1, -2 + 2*I: 1}

Las fracciones parciales de H(s) usando el algoritmo de la sección anterior se obtiene

H(s) = \frac{1}{5}\frac{4}{s+1}+\frac{1}{5} \frac{4s+17}{s^2 + 4s +8}

Ejemplo 1 Desarrollo con Sympy-Python

Con el algoritmo se tiene el siguiente resultado:

busca_PolosCeros:
 Q_polos : {-1: 1, -2 - 2*I: 1, -2 + 2*I: 1}
 P_ceros : {3: 1}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  0  unicos: 0

 asintoticamente:  estable

 h(t): 
     -t                   -2*t                            -2*t                
  4*e  *Heaviside(t)   9*e    *sin(2*t)*Heaviside(t)   4*e    *cos(2*t)*Heaviside(t)
- ------------------ + ----------------------------- + -----------------------------
          5                          10                              5        
>>>

En la gráfica se observa que todos los polos se encuentran en el lado izquierdo del plano LHP. Por lo que se considera que existen asintotas en la función correspondiente a tiempo, y el sistema es «asintóticamente estable«. También se considera BIBO estable por no tener los polos en el lado derecho RHP=0

Hs estabilidad polos Ej01

la gráfica de h(t) muestra lo indicado por el análisis de polos

Hs estabilidad ht Ej01

Instrucciones en Python Ejemplo 1

El algoritmo es una continuación del manejo de Polos y Ceros realizado en la sección anterior, donde se añade la revisión de las raíces de Q cuando son de tipo imaginarias.

Por simplicidad para el análisis se usan solo las partes esenciales, asi el enfoque es sobre la posición de los polos y su relación con h(t).

En el proceso, se agrupan las partes reales en Q_real y las imaginaria en Q_imag para realizar el conteo de lo repetido.

# Polos y ceros de funcion de transferencia H(s)
# Estabilidad del sistema, H(s) respuesta al impulso
# Ps es numerador, Qs es denominador
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
import telg1001 as fcnm  

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t',real=True)

Hs = ((s-3)/(s+1))*(1/(s**2+4*s+8))

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

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

# h(t) desde H_fp(s) con inversa Laplace 
ht = sym.inverse_laplace_transform(Hs,s,t)

# SALIDA

# Analiza estabilidad asintotica
casi_cero = 1e-12
# Separa parte real e imaginaria de raices
cuenta_real = 0; cuenta_imag = 0
unicos = 0 ; repetidos = 0 ; enRHP = 0
for raiz in Q_polos:
    [r_real,r_imag] = raiz.as_real_imag()
    if abs(r_real)>casi_cero and abs(r_imag)<casi_cero :
        cuenta_real = cuenta_real+1
    # para estabilidad asintotica
    conteo = Q_polos[raiz]
    if conteo==1 and r_real==0 and abs(r_imag)>0:
        unicos = unicos + 1
    if conteo>1  and r_real==0 and abs(r_imag)>0:
        repetidos = repetidos + 1
    if r_real>0:
        enRHP = enRHP + 1
cuenta_imag = len(Q_polos)-cuenta_real

# Revisa lado derecho del plano RHP
asintota = ""
if enRHP==0:
    asintota = 'estable'
if enRHP>0 or repetidos>0:
    asintota = 'inestable'
if enRHP==0 and unicos>0:
    asintota = 'marginalmente estable'

# SALIDA
print('busca_PolosCeros:')
for k in polosceros.keys():
    eq_lista = ['H','H_fp']
    if k in eq_lista:
        print('',k,':')
        sym.pprint(polosceros[k])
    else:
        print('\n',k,':',polosceros[k])

print(' polos reales: ',cuenta_real)
print(' polos complejos: ',cuenta_imag)
print(' sobre lado derecho RHP:',enRHP)
print(' sobre Eje Imaginario, repetidos: ',
      repetidos,' unicos:', unicos)
print('\n asintoticamente: ', asintota)

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

# GRAFICA -----------
muestras_H = 161
fig_Hs = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,
                          muestras=muestras_H,
                          f_nombre='H',
                          solopolos=True)
figura_h = fcnm.graficar_ft(ht,t_a,t_b,f_nombre='h')
plt.show()

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]
..


Ejemplo 2 H(s) con polos no repetidos, uno en RHP

Referencia: Lathi 2.14b p201

Revise la estabilidad BIBO de los siguientes sistemas LTIC descritos por las ecuaciones:

(D-1)(D^2 + 4D +8) y(t) = (D+2) x(t)

Ejemplo 2 Desarrollo analítico

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s-1)(s^2 + 4s +8) Y(s) = (s+2) X(s) \frac{Y(s)}{X(s)} = \frac{s+2}{(s-1)(s^2 + 4s +8)}

El grado del polinomio Q es N=3 y es mayor que el grado del polinomio P es M=1.

Los polos se obtienen como las raíces del denominador Q, necesario para plantear las fracciones parciales y simplificar la expresión en componentes más simples.

La expresión para el algoritmo se escribe como:

Hs = ((s+2)/(s-1))*(1/(s**2+4*s+8))

y el resultado de polos obtenido es:

polos: {1: 1, -2 - 2*I: 1, -2 + 2*I: 1}

Las fracciones parciales de H(s) usando el algoritmo de la sección anterior se obtiene

H(s) = \frac{3}{13}\frac{1}{s-1}+\frac{1}{13} \frac{3s+2}{s^2 + 4s +8}

Ejemplo 2 Desarrollo con Sympy-Python

Con el algoritmo se tiene el siguiente resultado:

busca_PolosCeros:
 H :
        s + 2         
----------------------
        / 2          \
(s - 1)*\s  + 4*s + 8/
 H_fp :
       3*s + 2            3     
- ----------------- + ----------
     / 2          \   13*(s - 1)
  13*\s  + 4*s + 8/             

 P_ceros : {-2: 1}

 Q_polos : {1: 1, -2 - 2*I: 1, -2 + 2*I: 1}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 1
 sobre Eje Imaginario, repetidos:  0  unicos: 0

 asintoticamente:  inestable

 h(t): 
   t                   -2*t                            -2*t                   
3*e *Heaviside(t)   2*e    *sin(2*t)*Heaviside(t)   3*e    *cos(2*t)*Heaviside(t)
----------------- + ----------------------------- - -----------------------------
        13                        13                              13          
>>>

En la gráfica se observa que un polo se encuentra en el lado derecho del plano RHP. Un polo en RHP genera un componente creciente en el tiempo, en consecuencia, el sistema es «asintóticamente inestable». También el sistema es BIBO inestable.
Hs estabilidad polos Ej02

la gráfica de h(t) complementa lo interpretado con la posición de los polos en el plano s

Hs estabilidad ht Ej02

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]

..


Ejemplo 3 H(s) con polos complejos sobre eje imaginario únicos

Referencia: Lathi 2.14c p201

Revise la estabilidad BIBO de los siguientes sistemas LTIC descritos por las ecuaciones:

(D+2)(D^2 +4)y(t) = (D^2+D+1)x(t)

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s+2)(s^2 + 4) Y(s) = (s^2+s+1) X(s)
Ps = s**2+s+1
Qs = (s+2)*(s**2+4)
Hs = Ps/Qs

El grado del polinomio Q es N=3 y el grado del polinomo P  es M=2, los polos del denominador muestran raíces complejas, por lo que se tendrá componentes cuadráticos para la transformada inversa de Laplace.

 Q_polos : {-2: 1, -2*I: 1, 2*I: 1}

Usando fracciones parciales con el algoritmo se tiene que,

H(s) = \frac{3}{8}\frac{1}{s+2}+\frac{1}{8} \frac{5s-2}{s^2 + 4}

resultados:

busca_PolosCeros:
 Q_polos : {-2: 1, -2*I: 1, 2*I: 1}
 P_ceros : {-1/2 - sqrt(3)*I/2: 1, -1/2 + sqrt(3)*I/2: 1}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  0  unicos: 2

 asintoticamente:  marginalmente estable

 h(t): 
                                                       -2*t             
  sin(2*t)*Heaviside(t)   5*cos(2*t)*Heaviside(t)   3*e    *Heaviside(t)
- --------------------- + ----------------------- + --------------------
            8                        8                       8          
>>>

al tener polos sobre el eje imaginario el sistema es asintóticamente marginal estable, pero BIBO-inestable.

Hs estabilidad polos Ej03

con gráfica de h(t)

Hs estabilidad ht Ej03

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]


Ejemplo 4. H(s) con polos repetidos sobre eje imaginario

Referencia: Lathi 2.14d p201

Revise la estabilidad BIBO de los siguientes sistemas LTIC descritos por las ecuaciones:

(D+1)(D^2 + 4)^2 y(t) = (D^2 + 2D + 8) x(t)

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s+1)(s^2 + 4)^2 Y(s) = (s^2 + 2s + 8) X(s)

Las raíces del denominador Q muestran que tienen valores complejos, mostrando que se usaría un componente cuadratico para la transformada inversa de Laplace.

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

al presentar dos raíces repetidas sobre el eje imaginario, se considera asintóticamente inestable y también BIBO inestable.

Revise el resultado de la inversa de Laplace usando las tablas.

los resultados con el algoritmo son:

busca_PolosCeros:
 H :
    2            
   s  + 2*s + 8  
-----------------
                2
        / 2    \ 
(s + 1)*\s  + 4/ 
 H_fp :
   2*(s - 6)     7*(s - 1)        7     
- ----------- - ----------- + ----------
            2      / 2    \   25*(s + 1)
    / 2    \    25*\s  + 4/             
  5*\s  + 4/                            

 P_ceros : {-1 - sqrt(7)*I: 1, -1 + sqrt(7)*I: 1}

 Q_polos : {-1: 1, -2*I: 2, 2*I: 2}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  2  unicos: 0

 asintoticamente:  inestable        
>>>

gráfica de polos

Hs estabilidad polos Ej04

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]


Tarea

Lathi ejercicios 2.15 p202


Ejemplo 5. Con suma de términos y exponenciales de retraso en tiempo

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

la expresión para el algoritmo se escribe como:

Hs = (s+2)/(s-1) + sym.exp(-2*s)/(s-2)

con lo que se obtiene el resultado, usando intervalo de tiempo [-1,3] para las gráficas.

exp(-2*s) : {'Q_polos': {2: 1}, 'P_ceros': {}, 'Hs_k': 1/(s - 2)}
1 : {'Q_polos': {1: 1}, 'P_ceros': {-2: 1}, 'Hs_k': (s + 2)/(s - 1)}
Q_polos : {2: 1, 1: 1}
P_ceros : {-2: 1}
n_polos_real : 2
n_polos_imag : 0
enRHP : 2
unicos : 0
repetidos : 0
asintota : inestable

 h(t): 
 -4  2*t                       t                             
e  *e   *Heaviside(t - 2) + 3*e *Heaviside(t) + DiracDelta(t)
>>>

con gráfica de polos

Hs estabilidad polos Ej05

con gráfica h(t)

Hs estabilidad ht Ej05

Instrucciones con Python

# Estabilidad del sistema, H(s) respuesta al impulso
# Polos y ceros de funcion de transferencia H(s)
# Ps es numerador, Qs es denominador
import sympy as sym
import telg1001 as fcnm  

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t',real=True)

Hs = (s+2)/(s-1) + sym.exp(-2*s)/(s-2)

# Ps = s**2+2*s+8
# Qs = (s+1)*((s**2+4)**2)
# Hs = Ps/Qs

#Ps = s**2+s+1
#Qs = (s+2)*(s**2+4)
#Hs = Ps/Qs

#Hs = ((s+2)/(s-1))*(1/(s**2+4*s+8))
#Hs = ((s-3)/(s+1))*(1/(s**2+4*s+8))

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

# PROCEDIMIENTO
Hs_fp = fcnm.apart_s(Hs)
polosceros = fcnm.busca_polosceros(Hs_fp)
P_ceros = polosceros['P_ceros']
Q_polos = polosceros['Q_polos']

def estabilidad_asintotica_s(Q_polos, casi_cero=1e-8):
    ''' Analiza estabilidad asintotica con Q_raiz
        Separa parte real e imaginaria de raices
        casicero es la tolerancia para considerar cero
    '''
    # Separa parte real e imaginaria de raices
    cuenta_real = 0; cuenta_imag = 0
    unicos = 0 ; repetidos = 0 ; enRHP = 0
    for raiz in Q_polos:
        [r_real,r_imag] = raiz.as_real_imag()
        if abs(r_real)>casi_cero and abs(r_imag)<casi_cero :
            cuenta_real = cuenta_real+1
        # para estabilidad asintotica
        conteo = Q_polos[raiz]
        if conteo==1 and r_real==0 and abs(r_imag)>0:
            unicos = unicos + 1
        if conteo>1  and r_real==0 and abs(r_imag)>0:
            repetidos = repetidos + 1
        if r_real>0:
            enRHP = enRHP + 1
    cuenta_imag = len(Q_polos)-cuenta_real

    # Revisa lado derecho del plano RHP
    asintota = ""
    if enRHP==0:
        asintota = 'estable'
    if enRHP>0 or repetidos>0:
        asintota = 'inestable'
    if enRHP==0 and unicos>0:
        asintota = 'marginalmente estable'
        
    estable = {'n_polos_real': cuenta_real,
               'n_polos_imag': cuenta_imag,
               'enRHP'     : enRHP,
               'unicos'    : unicos,
               'repetidos' : repetidos,
               'asintota'  : asintota,}
    return(estable)

estable = estabilidad_asintotica_s(Q_polos)

# h(t) desde H_fp(s) con inversa Laplace
ht = 0*t
term_suma = sym.Add.make_args(Hs_fp.expand())
for term_k in term_suma:
    ht_k = sym.inverse_laplace_transform(term_k,s,t)
    ht = ht +ht_k
ht = sym.expand(ht,t) # terminos suma

# SALIDA
fcnm.print_resultado_dict(polosceros)
fcnm.print_resultado_dict(estable)

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

# GRAFICA -----------
import matplotlib.pyplot as plt
muestras_H = 161
fig_Hs = fcnm.graficar_Fs(Hs_fp,Q_polos,P_ceros,
                          muestras=muestras_H,
                          f_nombre='H',solopolos=True)
figura_h = fcnm.graficar_ft(ht,t_a,t_b,f_nombre='h')
plt.show()

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]


H(s) Implicaciones de estabilidad y Osciladores

Referencia: Lathi 2.5-3 p203

Todos los sistemas de procesamiento de señales deben ser asintóticamente estables.

Los sistemas inestables no son requeridos debido a que con condiciones iniciales, intencionales o no, llevan a respuestas fuera de rango. Estas respuestas inestables destruyen el sistema o lo llevan a condiciones de saturación que cambian la naturaleza del sistema debido al tipo de crecimiento exponencial.

Sistemas marginalmente estables, aunque BIBO inestables, tienen aplicación como osciladores, que es un sistema que genera su propia señal sin la aplicación de entradas externas.

La salida de un oscilador es una respuesta de entrada cero. Si la respuesta es una sinusiode con frecuencia ω0  el sistema es marginalmente estable en ±jω0. Para diseñar un oscilador de frecuencia  se debe usar un sistema descrito por la ecuación diferencial:

(D^2 + \omega _0 ^2 ) y(t) = x(t)

Sin embargo, osciladores en la práctica se construyen con sistemas no lineales.