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.3 Periodicidad señales analógicas y digitales

Considere la exponencial compleja discreta con frecuencia (ω0 + 2π)

e^{j(\omega _0 + 2 \pi)n} =e^{j 2 \pi n} e^{j \omega _0 n} = e^{j \omega _0 n}

Se puede ver que la exponencial con frecuencia ω0 + 2π es la misma que aquella con la frecuencia ω0.

Existe una situación diferente en el caso contínuo, en el cual las señales e0t son todas distintas para distintos valores de ω0.

En el caso discreto, éstas señales no son distintas, ya que la señal con frecuencia ω0 es idéntica a las señales con frecuencias ω0 ± 2π , ω0 ± 4π, … y las que le siguen. Por lo tanto, al considerar las exponenciales complejas, necesitamos tomar en cuenta solamente un intervalo de frecuencia de longitud 2π dentro del cual se escoge ω0.

ω0N = 2π m

o de otra forma

ω0/2π = m/N

Para el caso de m=1 y N=2

otra prueba con m=7 y N=4

con m=2 y N=1

El código en python para realizar las gráficas es:

# Señales discretas, y(t)  a  y[t]
# propuesta:edelros@espol.edu.ec
import numpy as np
import math
import matplotlib.pyplot as plt

# Definir la funcion para el ejemplo
def analogica(f,t):
    # función matemática CAMBIAR AQUI
    y=np.cos(f*t)
    return(y)

# Programa para graficar la función
# INGRESO
rango = 4 # float(input('rangos en periodos de analogica:'))
m = 1     # float(input('frecuencia de analógica:'))
N = 2    # float(input('frecuencia digital:'))

# PROCEDIMIENTO
# grafica analogica
puntoscontinuos = 500
t0 = -rango*2*np.pi*m/2
tn = rango*2*np.pi*m/2
t  = np.linspace(t0,tn,puntoscontinuos+1)
yanalog = analogica(m,t)

# grafica digital
deltaD = (2*np.pi*m)/(2*N)
muestreo = int((tn-t0)//deltaD +1)
td = np.linspace(t0,tn,muestreo)
ydigital = analogica(m,td)

# SALIDA - GRAFICA
#Escala y para el grafico
margen = 0.1*np.max(yanalog)
ymax = np.max(yanalog)+margen
ymin = np.min(yanalog)-margen

plt.figure(1)
plt.suptitle('Señal Analogica vs digital')

plt.subplot(211)    # grafica de 2x1 arriba
plt.plot(t,yanalog)
plt.axis((t0,tn,ymin,ymax))
plt.ylabel('y(t)')

plt.subplot(212)    # grafica de 2x1 abajo
plt.plot(td,ydigital,'ro')
plt.axis((t0,tn,ymin,ymax))
plt.ylabel('Digital: y[t]')

plt.show()

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 ]

Tablas

Convolución

Convolución Integrales – Tabla de pares

Convolución – Tabla de Propiedades

Convolución de sumas – Tabla

Transformada de Laplace

Transformada de Laplace – Tabla de pares

Transformada de Laplace – Tabla de propiedades

Transformada z

Transformada z – Tabla de pares

Transformada z – Tabla de Propiedades

Transformada de Fourier

Transformada de Fourier – Tabla de pares

Transformada de Fourier – Tabla de Propiedades

 

Tabla de Fórmulas Matemáticas

Transformada de Laplace – Tabla de Propiedades

Referencia: Oppenheim tabla 9.3 p717 pdf 748. Lathi tabla 4.2 p361. Hsu tabla 3-2 p119

Propiedad Señal Transformada de Laplace
x(t)
x1(t)
x2(t)
X(s)
X1(s)
X2(s)
Linealidad ax1(t) + bx2(t) a X1(s) + b X2(s)
Desplazamiento en t x(t-t0) e-st0X(s)
Desplazamiento
en domino de s
es0tx(t) X(s-s0)
Escalamiento
en tiempo
x(at), a>0 \frac{1}{a} X\Big(\frac{s}{a}\Big)
inversión
en tiempo
x(-t) X(-s)
Conjugación x*(t) X*(s)
Convolución
(suponiendo que
x1(t) y x2(t) son cero
para t<0)
x_1 (t) \circledast x_2 (t) X1(s)X2(s)
Diferenciación en
el dominio t
\frac{\delta}{\delta t} x(t) sX(s)-x(0)
\frac{\delta^2}{\delta t^2} x(t) s2X(s)-sx(0)-x(0)
\frac{\delta^3}{\delta t^3} x(t) s3X(s) -s2x(0)-sx(0)-x'(0)
\frac{\delta^n}{\delta t^n} x(t) s^n X(s) -\sum_{k=1}^{n}s^{(n-k)}x^{(k-1)} (0^{-})
Diferenciación en
el dominio s
-tx(t) \frac{\delta}{\delta s} X(s)
Integración en
en el dominio t
\int_{0^-}^t x(\tau) \delta \tau \frac{1}{s} X(s)
\int_{-\infty}^t x(\tau) \delta \tau \frac{1}{s} X(s)+\frac{1}{s} \int_{-\infty}^{0^{-}} x(t) \delta t
Convolución en frecuencia x1(t) x2(t) \frac{1}{2 \pi j} X_1(s) \circledast X_2(s)
valor inicial x(0+) \lim _{s \rightarrow \infty} s X(s)
n>m
valor final x(∞) \lim _{s \rightarrow 0} s X(s)
polos de sX(s) en LHP (izquierda)

Transformada de Laplace – Tabla

Transformada de Laplace – Tabla

Referencia: Lathi Tabla 4.1 p334. Oppenheim Tabla 9.2 p692. Hsu Tabla 3-1 p115

No. x(t) X(s) ROC
1a δ(t) 1 Toda s
1b δ(t-T) e-sT Toda s
2a μ(t) \frac{1}{s} Re{s}>0
2b -μ(-t) \frac{1}{s} Re{s}<0
3 tμ(t) \frac{1}{s^2} Re{s}>0
4a tnμ(t) \frac{n!}{s^{n+1}} Re{s}>0
4b \frac{t^{n-1}}{(n-1)!} \mu (t) \frac{1}{s^n} Re{s}>0
4c -\frac{t^{n-1}}{(n-1)!} \mu (-t) \frac{1}{s^n} Re{s}<0
5 eλtμ(t) \frac{1}{s-\lambda} Re{s}>0
6 teλtμ(t) \frac{1}{(s-\lambda)^2} Re{s}>0
7 tneλtμ(t) \frac{n!}{(s-\lambda)^{n+1}}
8a cos (bt) μ(t) \frac{s}{s^2+b^2} Re{s}>0
8b sin (bt) μ(t) \frac{b}{s^2+b^2} Re{s}>0
9a e-atcos (bt) μ(t) \frac{s+a}{(s+a)^2+b^2} Re{s}>-a
9b e-atsin (bt) μ(t) \frac{b}{(s+a)^2+b^2} Re{s}>-a
10 \mu_n (t) = \frac{\delta ^n}{\delta t^n} \delta (t) sn Toda s
11 \mu_{-n} (t) = \mu (t) \circledast \text{...} \circledast \mu (t)

n veces

\frac{1}{s^n} Re{s}>0
12a re-atcos (bt+θ) μ(t) \frac{(r\cos (\theta)s + (ar \cos (\theta) - br \sin (\theta))}{s^2+2as+(a^2+b^2)}
12b re-atcos (bt+θ) μ(t) \frac{0.5 re^{j \theta}}{s+a-jb} + \frac{0.5 re^{-j \theta}}{s+a+jb}
12c re-atcos (bt+θ) μ(t) \frac{As+B}{s^2+2as+c}
r = \sqrt{\frac{A^2 c +B^2 -2ABa}{c-a^2}} \theta = \tan ^{-1} \Big( \frac{Aa-B}{A\sqrt{c-a^2}}\Big)

b = \sqrt{c-a^2}

12d e^{-at}\Bigg[A \cos (bt) + \frac{B-Aa}{b} \sin (bt) \Bigg] \mu (t) \frac{As+B}{s^2 + 2as+c}
b = \sqrt{c-a^2}

Transformada Laplace – Tabla de Propiedades

Transformada de Laplace – Concepto con Python

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

Referencia: Lathi 4.4-1 p382,  Oppenheim 11.50 p896

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

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

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

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

1. Amplificador Operacional en el dominio s

Referencia: Lathi 4.6-5 p399

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

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

1,1 Amplificador Operacional como multiplicador escalar

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

1,2 Amplificador Operacional como Integrador

Referencia: Oppenheim 11.52 p898

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

1,3 Amplificador Operacional como Sumador

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

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

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


Ejemplo 1. Implementación con Op-Amp

Referencia: Lathi 4.25 p401

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

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

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

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

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

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


Ejemplo 1. Desarrollo analítico

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

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

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

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

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

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

 X(s): 
1

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

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

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

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

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

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

OpAmpEj01 H(s) Polos

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

Las respuestas fueron obtenidas al usar como bloque de entrada,

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

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

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

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

4.4.4 LTI Laplace – Ejercicios con circuitos eléctricos y transformadas de Laplace

Ejemplos de solución a circuitos eléctricos que en el planteamiento se usan ecuaciones diferenciales ordinarias y se desarrollan usando la Transformada de Laplace. En el ejercicio 3 se desarrolla una malla que genera un sistema de ecuaciones EDO y se resuelve con transformadas de Laplace.

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]


Ejercicio 1. Circuito RLC, fuente DC e interruptor con transformada de Laplace

Referencia: Lathi Ej.4.13 p365

En el circuito de la figura, el interruptor se encuentra cerrado mucho tiempo antes de t=0. Cuando se abre en un instante, encuentre la corriente del inductor y(t) para t≥0.

El interruptor se encuentra cerrado por mucho tiempo, la corriente por el inductor es 2 A y el voltaje del capacitor es 10 V, pues el inductor en DC opera como un conductor sin resistencia y el capacitor se encuentra completamente cargado.

Cuando se abre el interruptor, el circuito se equivalente es el mostrado en la derecha, La corriente inicial del inductor es y(0)=2 y el voltaje inicial del capacitor Vc(0)=10.

El voltaje en la entrada es 10 V, empezando en t=0 y puede ser representado como 10 μ(t) para simplificar la representación de la fuente DC y expresar que antes de t=0 se aplica las condiciones iniciales dadas.mientras se den las condiciones iniciales en t=0 solo será necesario saber que la corriente en t≥0 para determinar la respuesta en t≥0 .

La ecuación del circuito en luego de abrir el interruptor es:

\frac{\delta}{\delta t}y(t) + 2 y8t) + 5 \int_{-\infty}^{t}y(\tau) \delta \tau = 10 \mu(t)

Las condiciones iniciales se aplican como:

\frac{\delta}{\delta t}y(t) = sY(s) - y(0^-)= sY(s)-2 \int_{-\infty}^{t} y(\tau) \delta \tau = \frac{1}{s}Y(s) + \frac{1}{s} \int_{-\infty}^{0^-}y(\tau) \delta \tau

y(t) es la corriente del capacitor, por lo que el integral es:

\int_{-\infty}^{0^-}y(\tau) \delta \tau = q_C(0^-) = CV_c (0^-) = \frac{1}{5}10 = 2

por lo que la parte del capacitor es:

\int_{-\infty}^{t}y(\tau) \delta \tau = \frac{1}{s}Y(s) + \frac{2}{s}

La transformada de Laplace de la ecuación integro diferencial del circuito RLC al usar los resultados se reescribe como:

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

El polinomio del denominador tiene raíces complejas

>>> import  sympy as sym
>>> s = sym.Symbol('s')
>>> Qs = s**2+2*s+5
>>> sym.roots(Qs)
{-1 - 2*I: 1, -1 + 2*I: 1}

por lo que es conveniente usar la forma cuadrática de la transformada de Laplace, donde A= 1, B=0, a=1, c=5, siendo,

r=\sqrt{\frac{20}{4}} = \sqrt{5} b=\sqrt{c-a^2} = 2 \theta=\tan^{-1} \Big(\frac{2}{4} \Big) = 0.46364 rad = 26.56^{o}

usando la tabla de transformadas,

y(t) = \sqrt{5} e^{-t} cos (2t+0.46364) \mu (t)

El gráfico de polos de la función de transferencia en el dominio s:

LTIC Laplace Circuito Eléctrico 01 Polos H(s)

que generan una respuesta de salida y(t)

LTI Laplace Circuito Electrico 01 xh_y

si la entrada x(t) es un impulso unitario, la respuesta es la misma que h(t)

La respuesta del algoritmo obtenida es:

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

H(s) parámetros cuadraticos:
2*s/(s**2 + 2*s + 5) : {'A': 2.0, 'B': 0, 'a': 1.0,
 'c': 5.0, 'r': 2.23606797749979, 'b': 2.0,
 'theta': 0.4636476090008061}

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

polosceros:
Q_polos : {-1 - 2*I: 1, -1 + 2*I: 1}
P_ceros : {0: 1}

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

 X(s): 
1

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

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

 ZSR_Qs2 :
 2*s/(s**2 + 2*s + 5) :
  {'A': 2.0, 'B': 0, 'a': 1.0, 'c': 5.0,
 'r': 2.23606797749979, 'b': 2.0,
 'theta': 0.4636476090008061}
yt_ZSR :
/   -t               -t         \             
\- e  *sin(2*t) + 2*e  *cos(2*t)/*Heaviside(t)

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

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

Instrucciones en Python

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

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

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

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

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

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

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

# 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_exp(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()

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]


Ejemplo 2. Circuito RLC con transformada de Laplace

Referencia: Lathi Ejemplo 4.17 p375

Realice el análisis de la corriente i(t) en el circuito mostrado en la figura, suponiendo que todas las condicioes iniciales son cero.

El primer paso es representar el circuito en el dominio de la frecuencia (s), mostrado a la derecha de la figura. Se representan los voltajes y corrientes usando la transformada de Laplace.

El voltaje 10μ(t) se representa como 10/s y la corriente i(t) como I(s). Todos los elementos del circuito se muestran con su respectiva impedancia. El inductor de 1 henrio se representa por s, el capacitor de 1/3 faradio se representa por 2/s. El resistor de 3 ohms es solo 3.

El voltaje en el circuito V(s) en cualquier elemento es I(s) por su impedancia Z(s).

La impedancia del circuito para el lazo es:

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

el voltaje de entrada es V(s)=10/s, por lo que la fórmula básica de corriente es:

I(s) =\frac{V(s)}{Z(s)} = \frac{10/s}{(s^2+3s+2)/s} I(s) = \frac{10}{(s^2+3s+2)}

usando las raíces del denominador,

I(s) =\frac{10}{(s+1)(s+2)}

aplicando fracciones parciales:

I(s) = \frac{10}{s+1} -\frac{10}{s+2}

Aplicando la transformada inversa,

i(t) = 10*(e^{-t}-e^{-2t}) \mu (t)

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]


Ejemplo 3. Circuito RC y RL, fuente DC e interruptor con transformada de Laplace

Referencia: Lathi ejemplo 4.18 p378

El interruptor está en posición cerrada por un largo tiempo antes de que sea abierto en t=0. Encuentre las corrientes y1(t) y y2(t) para t>0

Al revisar el circuito se observa que cuando el interruptor esta cerrado y se han alcanzado las condiciones de estado estable, el voltaje del capacitor Vc=16 V y la corriente del inductor y2= 4 .

Las fuentes de 20V y 4V se pueden también ordenar en forma equivalente

Cuando se abre el interruptor en t=0, las condiciones iniciales son Vc(0)=16 y y2(0)= 4 como se muestra en la versión con transformada de Laplace en la figura anterior. La ecuación diferencial del circuito para suma de voltajes en la malla en el lazo 1 y la tabla de propiedades de la transformada de Laplace para la integración aplicada al capacitor,

-\frac{20}{s} + \frac{1}{s}\Big[Y_1(s)+ \int_{-\infty}^{0^-} y_1(\tau) \delta \tau \Big] + \frac{1}{5} \Big[Y_1(s)-Y_2(s) \Big] = 0 -\frac{20}{s} + \frac{1}{s}\Big[Y_1(s) + 16 \Big] + \frac{1}{5} Y_1(s) - \frac{1}{5}Y_2(s) = 0 -\frac{20}{s} + \frac{1}{s}Y_1(s) +\frac{16}{s} + \frac{1}{5}Y_1(s)-\frac{1}{5}Y_2(s) = 0 \Big[ \frac{1}{s}+ \frac{1}{5}\Big] Y_1(s)-\frac{1}{5}Y_2(s) = \frac{4}{s} \frac{1}{5}\Big[ \frac{5+s}{s}\Big] Y_1(s)-\frac{1}{5}Y_2(s) = \frac{4}{s} \Big[ \frac{5+s}{s}\Big] Y_1(s)-Y_2(s) = \frac{20}{s}
>>> import sympy as sym
>>> s = sym.Symbol('s')
>>> ecuacion1 = -20/s+(1/s)*(Y1+16)+(1/5)*Y1-(1/5)*Y2
>>> ecuacion1.expand()
0.2*Y1 + Y1/s - 0.2*Y2 - 4/s

El voltaje inicial del capacitor se puede representar por un voltaje en serie 16/s y la corriente inicial del inductor de 4 A que representa una fuente de valor VL = Ly2(0) = 1/2(4) = 2. La ecuación diferencial del lazo 2 y usando la transformada de Laplace para la derivada con condiciones iniciales para el inductor,

\frac{1}{5}\Big[Y_2(s)-Y_1(s)\Big] + 1 Y_2(s) + \frac{1}{2} \Big[ sY_2(s) - Y_2(0^-)\big] = 0 \Big[\frac{1}{5}+ 1\Big] Y_2(s) -\frac{1}{5}Y_1(s) + \frac{1}{2} \Big[ sY_2(s) - 4 \Big] = 0 \frac{6}{5} Y_2(s) -\frac{1}{5}Y_1(s) + \frac{1}{2} sY_2(s) - 4\frac{1}{2} = 0 -\frac{1}{5}Y_1(s) +\Big[ \frac{6}{5} + \frac{1}{2} s\Big] Y_2(s) = 2 -Y_1(s) +\Big[ \frac{12+5s}{2}\Big] Y_2(s) = 10

el resultado obtenido com Sympy, como una instrucción adicional a la anterior,

>>> ecuacion2 = (1/5)*(Y2-Y1)+Y2+(1/2)*(s*Y2-4)
>>> ecuacion2.expand()
-0.2*Y1 + 0.5*Y2*s + 1.2*Y2 - 2.0

Con lo que hay que resolver el sistema de ecuaciones:

\begin{cases} \Big[ \frac{s+5}{s}\Big] Y_1(s)-Y_2(s) = \frac{20}{s} \\ -Y_1(s) +\Big[ \frac{12+5s}{2}\Big] Y_2(s) = 10 \end{cases}

multipicando la primera ecuación por (s/(5+s)) y sumando con la segunda

\begin{cases} Y_1(s)-\frac{s}{s+5}Y_2(s) = \frac{20}{s}\frac{s}{s+5} \\ -Y_1(s) +\Big[ \frac{12+5s}{2}\Big] Y_2(s) = 10 \end{cases} -\frac{s}{s+5}Y_2(s) + \frac{12+5s}{2}Y_2(s) = \frac{20}{s+5} + 10 \Big[ -\frac{s}{s+5} + \frac{12+5s}{2} \Big] Y_2(s) = 10\Big[\frac{2+(s+5)}{s+5}\Big] \Big[ \frac{-2s+(s+5)(12+5s)}{2(s+5)}\Big] Y_2(s) = 10\frac{2+s+5}{s+5} (-2s+12s+5s^2 +60+25s) Y_2(s) = 20(s+7) (5s^2 +35s +60) Y_2(s) = 20(s+7) 5(s^2 +7s +12) Y_2(s) = 20(s+7) Y_2(s) = 4\frac{s+7}{s^2 +7s +12}

las raices del denominador son s=-4 y s=-3

Y_2(s) = 4\frac{s+7}{(s+3)(s+4)}

usando fracciones parciales y el método de Heaviside,

k_1 = 4\frac{s+7}{\cancel{(s+3)}(s+4)} \Big|_{s=-3} = 4\frac{(-3)+7}{(-3+4)} = 4\frac{4}{1} = 16 k_2 = 4\frac{s+7}{(s+3)\cancel{(s+4)}}\Big|_{s=-4} = 4\frac{(-4)+7}{(-4+3)}=4\frac{3}{-1} = - 12

reescribiendo en fracciones parciales

Y_2(s) = \frac{16}{s+3} -\frac{12}{s+4}

usando la tabla de transformadas de Laplace:

y_2(t) = (16 e^{-3t}-12e^{-4t}) \mu (t)

de forma semejante se puede resolver para Y1(s) al reemplazar el resultado en la primera ecuación del sistema de ecuaciones,

\Big[ \frac{s+5}{s}\Big] Y_1(s)-Y_2(s) = \frac{20}{s} \Big[ \frac{s+5}{s}\Big] Y_1(s)-\Big[\frac{16}{s+3} -\frac{12}{s+4}\Big] = \frac{20}{s} Y_1(s) = \frac{s}{s+5} \Big[\frac{16}{s+3} -\frac{12}{s+4} + \frac{20}{s}\Big] Y_1(s) = \frac{s}{s+5}\frac{16}{s+3} - \frac{s}{s+5}\frac{12}{s+4} + \frac{s}{s+5}\frac{20}{s} Y_1(s) = 16\frac{s}{(s+5)(s+3)} - 12\frac{s}{(s+5)(s+4)} + \frac{20}{s+5}

realizando las fracciones parciales con método de Heaviside para los dos primeros términos,

>>> import sympy as sym
>>> s = sym.Symbol('s')
>>> ya = 16*s/((s+5)*(s+3))
>>> ya.apart()
40/(s + 5) - 24/(s + 3)
>>> yb = -12*s/((s+5)*(s+4))
>>> yb.apart()
-60/(s + 5) + 48/(s + 4)
>>> 
Y_1(s) = \frac{40}{s+5} -\frac{24}{s+3} - \frac{60}{s+5}+\frac{48}{s+4} + \frac{20}{s+5} Y_1(s) = -\frac{24}{s+3} +\frac{48}{s+4}

aplicando desde la tabla la transformada inversa de Laplace,

y_1(t) = (-24e^{-3t} +48e^{-4t} ) \mu (t)

Instrucciones con Python
la solución del sistema de ecuaciones con transformadas de Laplace se puede realizar con Sympy

Y1 (s):
  24*s + 48  
-------------
 2           
s  + 7*s + 12

Y1 (s) en fracciones parciales:
  48      24 
----- - -----
s + 4   s + 3
Y2 (s):
   4*s + 28  
-------------
 2           
s  + 7*s + 12

Y2 (s) en fracciones parciales:
    12      16 
- ----- + -----
  s + 4   s + 3
>>> 
# Sistemas de ecuaciones con Sympy
import sympy as sym

# INGRESO
s = sym.Symbol('s')
[Y1, Y2] = sym.symbols('Y1 Y2')

ecuacion1 = ((s+5)/s)*Y1 - Y2 -20/s
ecuacion2 = -Y1 + +((12+5*s)/2)*Y2 -10

variables = [Y1,Y2]

# PROCEDIMIENTO
respuesta = sym.solve([ecuacion1,
                       ecuacion2],
                       variables)

# SALIDA
for cadarespuesta in respuesta:
    print(str(cadarespuesta) +'(s):')
    sym.pprint(respuesta[cadarespuesta])
    print('')
    print(str(cadarespuesta)+'(s) en fracciones parciales:') 
    sym.pprint(respuesta[cadarespuesta].apart())

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]