Categoría: Unidad 5 Sistemas LTI – Fourier

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:

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)

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 ]

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}

siendo el periodo T₀ = π y la frecuencia fundamental f₀ = 1/T₀  = 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,

Espectro de frecuencias,

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}

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)

Espectro de frecuencias de f(t)

[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:

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

el espectro de frecuencias,

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

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}

---

Coeficientes a_k y b_k

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 a_k y b_k.

Si requiere revisar el concepto se adjunta el enlace:

Sympy – Fórmulas y funciones simbólicas

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 a_k y b_k, 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 a_k y b_k, 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 t_i.

Para evaluar los puntos t_i 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()

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.

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 ]