Transformadas de Fourier – Tabla de Propiedades

Referencia: Schaum Hsu Tabla 5-1 p222. Lathi Tabla 7.2 p717, Oppenheim Tabla 4.1p328 pdf356

Transformada Fourier – Tabla

Series de Fourier de señales periódicas con Python

Operación x(t) X(ω)
x(t)
x1(t)
x2(t)
X(ω)
X1(ω)
X2(ω)
Multiplicación por escalar k x(t) k X(ω)
Aditiva x1(t) +x2(t) X1(ω) + X2(ω)
Linealidad a1 x1(t) + a2 x2(t) a1 X1(ω) + a2 X2(ω)
Conjugada x*(t) X*(-ω)
Inversión en tiempo x(-t) X(-ω)
Dualidad X(t) 2π x(-ω)
Escalamiento
(a Real)
x(at) \frac{1}{|a|} X \Big(\frac{\omega}{a} \Big)
Desplazamiento en tiempo (t Real) x(t-t0) X(ω) e-jωt0
Deplazamiento en frecuencia (ω0 Real) x(t)e0t X(ω-ω0)
Convolución en tiempo x_1(t) \circledast x_2 (t) X1(ω)X2(ω)
Convolución en frecuencia x1(t)x2(t) \frac{1}{2 \pi}X_1(\omega) \circledast X_2 (\omega)
Derivada en tiempo \frac{\delta^n}{\delta t^n}x(t) (jω)n X(ω)
Derivada en frecuencia (-jt)x(t) \frac{\delta}{\delta \omega}X(\omega)
Integral en tiempo \int_{- \infty}^{t}x(u) \delta u \frac{X(\omega)}{j \omega} + \pi X(0) \delta (\omega)
Señal Real x(t) = x_e(t) + x_o(t) X(\omega) = A(\omega)+jB(\omega)
componente par
componente impar
x_e(t)
x_o(t)
Re[X(\omega)] = A(\omega)
j Im[X(\omega)] = j B(\omega)

Relación de Parseval

\int_{-\infty}^{\infty} x_1 (\lambda) X_2 (\lambda) \delta \lambda = \int_{-\infty}^{\infty} X_1 (\lambda) x_2 (\lambda) \delta \lambda \int_{-\infty}^{\infty} x_1 (t) x_2 (t) \delta t = \frac{1}{2 \pi}\int_{-\infty}^{\infty} X_1 (\omega) X_2 (-\omega) \delta \omega \int_{-\infty}^{\infty} |x (t)|^2 \delta t = \frac{1}{2 \pi} \int_{-\infty}^{\infty} |X (\omega) |^2 \delta \omega

Transformada Fourier – Tabla

Series de Fourier de señales periódicas con Python

Transformadas de Fourier – Tabla

Referencia: Lathi Tabla 7.1 p699. Shaum Hsu Tabla 5-2 p223. Oppenheim Tabla 4.2 p329 pdf357

Transformada Fourier – Propiedades

Series de Fourier de señales periódicas con Python

No x(t) X(ω)
1 e-at μ (t) \frac{1}{a + j \omega} a>0
2 eat μ (-t) \frac{1}{a - j \omega} a>0
3 e-a|t| \frac{2a}{a^2+ \omega ^2} a>0
4 t e-at μ (t) \frac{1}{(a+j \omega)^2} a>0
5 tn  e-at μ (t) \frac{n!}{(a+j \omega)^{n+1}} a>0
6a δ(t) 1
6b δ(t-t0) e-jωt0
7 1 2\pi \delta (\omega)
8 e0t 2\pi \delta (\omega-\omega_0)
9 cos (ω0 t) \pi [\delta (\omega - \omega_0) +\delta (\omega + \omega_0)]
10 sin (ω0 t) \pi [\delta (\omega + \omega_0) -\delta (\omega - \omega_0)]
11a μ(t) \pi \delta (\omega ) +\frac{1}{j \omega }
11b μ(-t) \pi \delta (\omega ) - \frac{1}{j \omega }
12 sgn (t) \frac{2}{j \omega}
13 cos (ω0 t) μ(t) \frac{\pi}{2} [\delta (\omega - \omega_0) +\delta (\omega + \omega_0)] + \frac{j \omega}{\omega_0^2 - \omega ^2}
14 sin (ω0 t) μ(t) \frac{\pi}{2j} [\delta (\omega - \omega_0) - \delta (\omega + \omega_0)] + \frac{\omega_0}{\omega_0^2 - \omega ^2}
15 e-at sin (ω0 t) μ(t) \frac{\omega_0}{(a+j\omega)^2 + \omega_0^2} a>0
16 e-at cos (ω0 t) μ(t) \frac{a + j\omega}{(a+j\omega)^2 + \omega_0^2} a>0
17 rect \Big(\frac{1}{\tau}\Big) \tau sinc \Big( \frac{\omega \tau}{2} \Big)
18 \frac{W}{\pi} sinc (Wt) rect \Big(\frac{\omega}{2W} \Big)
19 \Delta \Big( \frac{t}{\tau} \Big) \frac{\tau}{2}sinc ^2 \Big( \frac{\omega \tau}{4} \Big)
20 \frac{W}{2\pi} sinc ^2 \Big(\frac{Wt}{2} \Big) \Delta \Big(\frac{\omega}{2W} \Big)
21 \sum_{n=- \infty}^{\infty} \delta(t-nT) \omega_0 \sum_{n=-\infty}^{\infty} \delta(\omega-n \omega_0) \omega_0 = \frac{2 \pi}{T}
22 \frac{e^{-t^2}}{2 \sigma ^2} \sigma \sqrt{2 \pi} e^{-\sigma^2 \omega^2 /2}
23 \frac{1}{a^2 + t^2} e ^{-a|\omega|}
24 e ^{-at^2} \sqrt{\frac{\pi}{a}}e ^{-\omega^2 /4a} a>0
25 p_a(t) = \begin{cases} 1 & |t|<a \\ 0 & |t|>a \end{cases} 2a\frac{\sin (\omega a)}{(\omega a)}

Transformada Fourier – Propiedades

Series de Fourier de señales periódicas con Python

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

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(j\omega) = H(j \omega) X(j\omega)

Lo que para señales y sistemas representa un componente importante en el análisis de sistemas LTI, al poder representar los sitemas 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(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(j \omega) = 1 e^{-j \omega t_0} = e^{-j \omega t_0}

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

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

Ejemplo 2. 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 propiedad de diferenciación

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

en consecuencia:

H(j \omega) = j \omega

Ejemplo 3. 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(j \omega) = \frac{1}{j \omega} +\pi \delta(\omega) Y(j \omega) = \Big[\frac{1}{j \omega} +\pi \delta(\omega) \Big]X(j \omega) Y(j \omega) = \frac{1}{j \omega}X(j \omega) +\pi \delta(\omega) X(0)

Ejemplo 4. 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(j \omega) = \frac{1}{a+ j \omega} X(j \omega) = \frac{1}{b+ j \omega} Y(j \omega) = H(j \omega)X(j \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(j \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(j \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]

En el algoritmo se procede de forma semejante a usar H(s). La parte gráfica se trata como un elemento separado del problema.

Instrucciones en Python

# H(w) fracciones parciales
# Ejemplo Oppenjeim 4.17 p318
# Qw Y(w) = Pw X(w)
# H(w)=Pw/Qw, numerador Pw, denominador Qw
# http://blog.espol.edu.ec/telg1001/

import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

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

# H(w) respuesta impulso
Pw = 1+0*w
Qw = a + sym.I*w

# X(w) Señal de entrada
Xw = 1/(b+sym.I*w)

# PROCEDIMIENTO
Pw = Pw.as_poly(w)
Qw = Qw.as_poly(w)

# Respuesta impulso
Hw = Pw/Qw
Hwp = sym.apart(Hw,w)

# Y(w) a entrada X(w)
Yw = Hw*Xw
Ywp = Yw.apart(w)

# SALIDA - polinomio
print('H(w):')
sym.pprint(Hw)

print('\nH(w) fracciones parciales:')
sym.pprint(Hwp)

print('\n Y(w)')
sym.pprint(Yw)

print('\n Y(w) total en fracciones parciales:')
sym.pprint(Ywp)

Las gráficas se desarrollan de acuerdo al caso del ejercicio, añadiendo instrucciones como:

# grafica de Y(w)-------------------------
a1 = 3 ; b1 = 2
muestras = 51
a_wlim = a1*4
wi  = np.linspace(-a_wlim,a_wlim,muestras)

# Y(w) real e imaginaria
Fwa = sym.lambdify(w,Yw.subs({a:a1,b:b1}))
Fwi = Fwa(wi) # evalua wi
if Yw.is_constant():
    Fwi = np.ones(len(wi))*Yw # evalua wi

# grafica Y(w)
figura, graf_Fwi = plt.subplots(2,1)

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

graf_Fwi[1].plot(wi,np.imag(Fwi),label='Imag(Y(w))',
                 color='magenta')
graf_Fwi[1].legend()
graf_Fwi[1].set_xlabel('w')
graf_Fwi[1].set_ylabel('imag(Y(w))')
graf_Fwi[1].grid()
etiqueta = ' ; a='+str(a1)+' ; b='+str(b1)
plt.suptitle(r'Y(w) = $'+str(sym.latex(Ywp))+'$'+etiqueta)
#plt.show()

# f(t)  --------------------------
# intervalo de gráfica [a_lim,b_lim]
T1 = 1
a_tlim = -2*T1
b_tlim = 2*T1
muestras = 51
ti  = np.linspace(a_tlim,b_tlim,muestras)
xt = lambda t: np.exp(-b1*t)*np.heaviside(t,0)
xti = xt(ti)
ht = lambda t: np.exp(-a1*t)*np.heaviside(t,0)
hti = ht(ti)
yt = lambda t: (1/(b1-a1))*(np.exp(-a1*t)*np.heaviside(t,0)-np.exp(-b1*t)*np.heaviside(t,0))
yti = yt(ti)

# grafica en dominio de tiempo
figura, graf_yt = plt.subplots()

graf_yt.plot(ti,xti,label='x(t)', color='blue')
graf_yt.plot(ti,hti,label='h(t)', color='magenta')
graf_yt.plot(ti,yti,label='y(t)', color='orange')
graf_yt.legend()
graf_yt.set_xlabel('t')
graf_yt.grid()
graf_yt.set_title('y(t) con h(t) y x(t)')
plt.show()

5.2 Transformada de Fourier – Señales aperiódicas continuas en tiempo con Python

Referencia: Oppenjeim 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.

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(j\omega) =\int_{-\infty}^{\infty} e^{-at} \mu (t) e^{-j \omega t} \delta t X(j\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(j\omega) =-\frac{1}{a+j\omega} e^{-(a+j\omega)(\infty)} +\frac{1}{a+j\omega} e^{-(a+j\omega)(0)} X(j\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
# integral de Fourier
# blog.espol.edu.ec/telg1001/
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
j = sym.I ; W = sym.oo
w = sym.Symbol('w',real=True)

t = sym.Symbol('t', real=True,)
a = sym.Symbol('a', real=True,positive=True)
b = sym.Symbol('b', real=True,positive=True)
T = sym.Symbol('T', real=True,positive=True)

ft = sym.exp(-a*t)*sym.Heaviside(t) # Ej 4.1 Oppenheim
[ta,tb] = [-sym.oo,sym.oo] # Intervalo Integración

#ft = sym.exp(-a*np.abs(t)) # Ej 4.2 Oppenheim
#[ta,tb] = [-sym.oo,sym.oo] # Intervalo Integración

# parametros para grafica
a1 = 2 # valor a
T1 = 1 # valor T o periodo

# intervalo de gráfica [a_lim,b_lim]
a_tlim = -2*T1
b_tlim = 2*T1
muestras = 52

# PROCEDIMIENTO
ftw = ft*sym.exp(-j*w*t)

# integral general en intervalo t
Fw_F = sym.integrate(ftw,(t,ta,tb))
if Fw_F.is_Piecewise:
    Fw = Fw_F.args[0] # primer intervalo
    Fw = Fw[0] # expresión buscada
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)


# GRAFICA valores ----------------
# f(t)
fta = sym.lambdify(t,ft.subs({a:a1,T:T1}))
ti  = np.linspace(a_tlim,b_tlim,muestras)
fti = fta(ti)

# F(w) magnitud y fase
Fw_mg = sym.lambdify(w,Fw_magn.subs({a:a1,T:T1}))
Fw_fs = sym.lambdify(w,Fw_fase.subs({a:a1,T:T1}))
a_wlim = a1*4
wi  = np.linspace(-a_wlim,a_wlim,muestras)
Fwi_mg = Fw_mg(wi) # evalua wi
Fwi_fs = Fw_fs(wi) 
if Fw_fase.is_constant():
    Fwi_fs=np.ones(len(wi))*Fw_fase

# F(w) real e imaginaria
Fwa = sym.lambdify(w,Fw.subs({a:a1,T:T1}))
a_wlim = a1*4
wi  = np.linspace(-a_wlim,a_wlim,muestras)
Fwi = Fwa(wi) # evalua wi

# f(t) dominio t
figura, graf_fti = plt.subplots()
plt.plot(ti,fti,label='f(t)',color='blue')
plt.xlabel('t')
plt.ylabel('f(t)')
etiq = str(a1)+'; T='+str(T1)
plt.title(r'f(t) = $'+str(sym.latex(ft))+'$ ; a='+etiq)
plt.grid()

# F(w) real e imaginaria
figura, 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
figura, graf_Fw = plt.subplots(2,1)

graf_Fw[0].plot(wi,Fwi_mg,label='F(w)_magn',
                color='orange')
Fwa0 = Fw_mg(0)
Fwa1 = Fw_mg(-a1)
Fwa2 = Fw_mg(a1)
graf_Fw[0].stem(-a1,Fwa1,linefmt ='--')
graf_Fw[0].stem(a1,Fwa2,linefmt ='--')
etiqueta1 = '('+str(a1)+','+str(np.round(Fwa2,4))+')'
graf_Fw[0].annotate(etiqueta1, xy=(a1,Fwa2))
etiqueta0 = '('+str(0)+','+str(np.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_fs,label='F(w)_fase',
                 color='brown')
graf_Fw[1].axhline(np.pi/2,linestyle ='--')
graf_Fw[1].axhline(-np.pi/2,linestyle ='--')
Fwa1f = Fw_fs(-a1)
Fwa2f = Fw_fs(a1)
graf_Fw[1].stem(-a1,Fwa1f,linefmt ='--')
graf_Fw[1].stem(a1,Fwa2f,linefmt ='--')
etiqueta3 = '('+str(a1)+','+str(np.round(Fwa2f,4))+')'
graf_Fw[1].annotate(etiqueta3, xy=(a1,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()

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(j\omega) = \int_{-\infty}^{\infty}e^{-a|t|} e^{-j \omega t} \delta t X(j\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(j\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*np.abs(t)) # Ej 4.2 Oppenheim
[ta,tb] = [-sym.oo,sym.oo] # Intervalo Integración

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


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(j \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(j \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
[ta,tb] = [-sym.oo,sym.oo] # Intervalo tiempo

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


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.

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}

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)

las amplitudes y fases se muestran en la tabla.

 coeficientes: 
 [n,		 ak,	   bk,		 Ck,	   fase_k ]
[[ 0.00e+00  5.0428e-01  0.0000e+00  5.0428e-01 -0.0000e+00]
 [ 1.00e+00  5.9327e-02  2.3731e-01  2.4461e-01 -1.3258e+00]
 [ 2.00e+00  1.5516e-02  1.2413e-01  1.2510e-01 -1.4464e+00]
 [ 3.00e+00  6.9556e-03  8.3467e-02  8.3756e-02 -1.4877e+00]
 [ 4.00e+00  3.9244e-03  6.2790e-02  6.2912e-02 -1.5084e+00]
 [ 5.00e+00  2.5151e-03  5.0302e-02  5.0365e-02 -1.5208e+00]
 [ 6.00e+00  1.7479e-03  4.1950e-02  4.1987e-02 -1.5292e+00]
 [ 7.00e+00  1.2848e-03  3.5974e-02  3.5997e-02 -1.5351e+00]
 [ 8.00e+00  9.8396e-04  3.1487e-02  3.1502e-02 -1.5396e+00]
 [ 9.00e+00  7.7761e-04  2.7994e-02  2.8005e-02 -1.5430e+00]
 [ 1.00e+01  6.2996e-04  2.5198e-02  2.5206e-02 -1.5458e+00]]
>>> 

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

Serie de Fourier fs(t)= 
0.504279523791290
+ 0.0593270027989753 cos(2.0 t) + 0.237308011195901 sin(2.0 t)
+ 0.0155162930397320 cos(4.0 t) + 0.124130344317856 sin(4.0 t)
+ 0.00695557963850054 cos(6.0 t) + 0.0834669556620066 sin(6.0 t)
+ 0.00392435427074932 cos(8.0 t) + 0.0627896683319894 sin(8.0 t)
+ 0.00251510984434557 cos(10.0 t) + 0.0503021968869117 sin(10.0 t)
+ 0.00174793595768210 cos(12.0 t) + 0.0419504629843707 sin(12.0 t)
+ 0.00128478859564658 cos(14.0 t) + 0.0359740806781048 sin(14.0 t)
+ 0.000983960046422013 cos(16.0 t) + 0.0314867214855049 sin(16.0 t)
+ 0.000777609134604903 cos(18.0 t) + 0.0279939288457771 sin(18.0 t)
+ 0.000629955682437573 cos(20.0 t) + 0.0251982272975036 sin(20.0 t)

El bloque de ingreso para el algoritmo es:

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

T0 = np.pi
# intervalo
at_lim = 0
bt_lim = at_lim + T0

ft = sym.exp(-t/2)

# número de coeficientes
n = 11
# Evaluación para gráfica
muestras = 101

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

Espectro de frecuencias,

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

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

T0 = 2
# intervalo
at_lim = -T0/4
bt_lim = at_lim + T0

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úmero de coeficientes
n = 9
# Evaluación para gráfica
muestras = 101

coeficientes de la serie de Fourier, magnitud y fase:

 coeficientes: 
 [ k, 		 ak, 	  bk ] 
[[ 0.0000e+00  0.0000e+00  0.0000e+00]
 [ 1.0000e+00  0.0000e+00  1.6211e+00]
 [ 2.0000e+00 -1.5593e-16  0.0000e+00]
 [ 3.0000e+00  0.0000e+00 -1.8013e-01]
 [ 4.0000e+00  1.5593e-16  0.0000e+00]
 [ 5.0000e+00  1.7764e-17  6.4846e-02]
 [ 6.0000e+00 -1.5593e-16  0.0000e+00]
 [ 7.0000e+00  1.8126e-17 -3.3084e-02]
 [ 8.0000e+00  1.5593e-16  0.0000e+00]]

Serie de Fourier en forma trigonométrica:

Serie de Fourier fs(t)= 
0
+ 0 cos(3.141592653589793 t) + 1.62113893827740 sin(3.141592653589793 t)
+ -1.55926873300775e-16 cos(6.283185307179586 t) + 0 sin(6.283185307179586 t)
+ 0 cos(9.42477796076938 t) + -0.180126548697489 sin(9.42477796076938 t)
+ 1.55926873300775e-16 cos(12.566370614359172 t) + 0 sin(12.566370614359172 t)
+ 1.77635683940025e-17 cos(15.707963267948966 t) + 0.0648455575310958 sin(15.707963267948966 t)
+ -1.55926873300775e-16 cos(18.84955592153876 t) + 0 sin(18.84955592153876 t)
+ 1.81260901979617e-17 cos(21.991148575128552 t) + -0.0330844681281103 sin(21.991148575128552 t)
+ 1.55926873300775e-16 cos(25.132741228718345 t) + 0 sin(25.132741228718345 t)

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}
# INGRESO
t = sym.Symbol('t')

T0 = 2*np.pi
# intervalo
at_lim = -T0/2
bt_lim = at_lim + T0

A = 2
ft = sym.Piecewise((0, t<-(T0/4)),
                   (1, t<(T0/4)),
                   (0, t<(T0/2)),
                   (0, True))
                   
# numero de coeficientes
n = 9
# Evaluacion para grafica
muestras = 101

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

los coeficientes de la serie de Fourier

 coeficientes: 
 [k,		 ak,	   bk,		 Ck,	   fase_k ]
[[ 0.0000e+00  5.0000e-01  0.0000e+00  5.0000e-01 -0.0000e+00]
 [ 1.0000e+00  6.3662e-01  0.0000e+00  6.3662e-01 -0.0000e+00]
 [ 2.0000e+00  3.8982e-17  0.0000e+00  3.8982e-17 -0.0000e+00]
 [ 3.0000e+00 -2.1221e-01  0.0000e+00  2.1221e-01  0.0000e+00]
 [ 4.0000e+00 -3.8982e-17  0.0000e+00  3.8982e-17  0.0000e+00]
 [ 5.0000e+00  1.2732e-01  0.0000e+00  1.2732e-01 -0.0000e+00]
 [ 6.0000e+00  3.8982e-17  0.0000e+00  3.8982e-17 -0.0000e+00]
 [ 7.0000e+00 -9.0946e-02  0.0000e+00  9.0946e-02  0.0000e+00]
 [ 8.0000e+00 -3.8982e-17  0.0000e+00  3.8982e-17  0.0000e+00]]

la serie de Fourier en forma trigonométrica

Serie de Fourier fs(t)= 
0.500000000000000
+ 0.636619772367581 cos(1.0 t) + 0 sin(1.0 t)
+ 3.89817183251938e-17 cos(2.0 t) + 0 sin(2.0 t)
+ -0.212206590789194 cos(3.0 t) + 0 sin(3.0 t)
+ -3.89817183251938e-17 cos(4.0 t) + 0 sin(4.0 t)
+ 0.127323954473516 cos(5.0 t) + 0 sin(5.0 t)
+ 3.89817183251938e-17 cos(6.0 t) + 0 sin(6.0 t)
+ -0.0909456817667973 cos(7.0 t) + 0 sin(7.0 t)
+ -3.89817183251938e-17 cos(8.0 t) + 0 sin(8.0 t)

comparando la serie de Fourier con f(t):

el espectro de frecuencias,


Algoritmo en Python

# Serie de Fourier, espectro con n coeficientes
# Lathi ej 6.1 p599

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

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

T0 = 2
# intervalo
at_lim = -T0/4
bt_lim = at_lim + T0

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úmero de coeficientes
n = 9
# Evaluación para gráfica
muestras = 101

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

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

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

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

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

# evalua los coeficientes
a0  = ak.subs(k,0)/2
b0  = bk.subs(k,0)
aki = [a0]
bki = [b0]
n_i = [0]
i = 1
while not(i>=n):
    avalor = ak.subs(k,i)
    bvalor = bk.subs(k,i)
    aki.append(avalor)
    bki.append(bvalor)
    n_i.append(i)
    i = i+1

# tabla de valores ak, Bk
tabla = np.concatenate([[np.array(n_i,dtype = int)],
                        [np.array(aki,dtype = float)],
                        [np.array(bki,dtype = float)]],
                       axis=0)
tabla = np.transpose(tabla)
np.set_printoptions(precision=4)
print('\n coeficientes: \n [ k, \t\t ak, \t  bk ] ')
print(tabla)

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

print('\nSerie de Fourier fs(t)= ')
#sym.pprint(termino)
print(a0)
for i in range(1,n,1):
    termino = '+ '+str(aki[i])+' cos('+str(i*w0)+' t)'+' + '
    termino = termino +str(bki[i])+' sin('+str(i*w0)+' t)'
    print(termino)
    


# Grafica - evaluación numérica
ti = np.linspace(at_lim,bt_lim,muestras)
ftn = sym.lambdify(t,ft)
fi = ftn(ti)

fsn = sym.lambdify(t,serieF)
fsi = fsn(ti)

tk = np.linspace(-2*T0,2*T0,4*muestras)
fk = np.zeros(len(tk))
# ajuste de intervalo de periodo
tki = (tk-at_lim)%T0+at_lim
fk = ftn(tki)

# Grafica f(t) varios periodos
figura, graf_ft = plt.subplots()
graf_ft.plot(tk,fk,label='f(t)',
             color= 'blue',linestyle='dashed')
graf_ft.plot(ti,fi,label='f(t)',color= 'blue')
graf_ft.axvline(at_lim, color ='red')
graf_ft.axvline(bt_lim, color ='red')
graf_ft.set_xlabel('t')
graf_ft.set_ylabel('x(t)')
graf_ft.legend()
graf_ft.grid()
etiqueta = '  , T0='+str(np.round(T0,4))
graf_ft.set_title('x(t) ='+str(ft)+etiqueta)

#plt.show()

# Grafica f(t) y Fourier en t
figura, graf_ftT0 = plt.subplots()
graf_ftT0.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
graf_ftT0.plot(ti,fsi,label = etiqueta)
graf_ftT0.set_xlabel('t')
graf_ftT0.set_ylabel('f(t)')
graf_ftT0.legend()
graf_ftT0.grid()
graf_ftT0.set_title('Serie de Fourier de f(t)')
#plt.show()

# espectro de frecuencia
k_i = np.arange(0,n,1,dtype=float)
ak_i = np.array(aki,dtype=float)
bk_i = np.array(bki,dtype=float)
ck_i = np.sqrt(ak_i**2 + bk_i**2)
cfs_i = np.arctan(-bk_i/ak_i)

tabla = np.concatenate([tabla,
                        np.transpose([ck_i]),
                        np.transpose([cfs_i])                        ],
                       axis=1)
np.set_printoptions(precision=4)
print('\n coeficientes: ')
print(' [k,\t\t ak,\t   bk,\t\t Ck,\t   fase_k ]')
print(tabla)

# grafica de espectro de frecuencia
figura, graf_Fw = plt.subplots(2,1)
graf_Fw[0].stem(k_i,ck_i,label='Ck_magnitud')
graf_Fw[0].set_ylabel('Ck')
graf_Fw[0].set_title('Espectro de frecuencia')
graf_Fw[0].legend()
graf_Fw[0].grid()

graf_Fw[1].stem(k_i,cfs_i,label='Ck_fase')
graf_Fw[1].legend()
graf_Fw[1].set_ylabel('Ck_fase')
graf_Fw[1].set_xlabel('k')
graf_Fw[1].grid()

plt.show()

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

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

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

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


Ejemplo

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

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 ak y bk

Para determinar las expresiones de los coeficientes, en Python se usa la libreria 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:
/1.2732*sin(1.5707*k) - 0.6366*sin(3.1415*k)  for And(k > -oo, k < oo, k != 0)
|--------------------------------------------
<                      k                                
|                                                                  
\          0                                    otherwise 

 expresión bk:
0

Nota: se han truncado los decimales a 4 para facilitar la lectura en la pantalla.

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

a_k = \begin{cases} \frac{1.2732\sin (1.5707 k) - 0.6366 \sin(3.1415 k )}{k} & \text{for}\: k > -\infty \wedge k < \infty \wedge k \neq 0 \\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 numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
T = 2*np.pi

t = sym.Symbol('t')
ft = sym.Piecewise((-1,t <-T/2),
                   (-1,t <-T/4),
                   ( 1,t < T/4),
                   (-1,t < T/2),
                   (-1,True),)
# intervalo
a = -T/2
b = T/2

# número de coeficientes
n = 4

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

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

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

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

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

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: 
[0, 1.27323954473516, 1.55926873300775e-16, 
  -0.424413181578388]

 coeficientes 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  = ak.subs(k,0)/2
b0  = bk.subs(k,0)
aki = [a0]
bki = [b0]
i = 1
while not(i>=n):
    avalor = ak.subs(k,i)
    bvalor = bk.subs(k,i)
    aki.append(avalor)
    bki.append(bvalor)
    i = i+1
print('\n coeficientes ak: ')
print(aki)
print('\n coeficientes bk: ')
print(bki)


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): 
1.27323954473516*cos(1.0*t) 
+ 1.55926873300775e-16*cos(2.0*t) 
- 0.424413181578388*cos(3.0*t)

Instrucciones en Python

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

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

print('\n serie de Fourier fs(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,

# Para evaluación numérica
fsn = sym.lambdify(t,serieF)
ftn = sym.lambdify(t,ft)

# Evaluación para gráfica
muestras = 41
ti = np.linspace(a,b,muestras)
fi = ftn(ti)
fsi = fsn(ti) 

# SALIDA
# Grafica
plt.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
plt.plot(ti,fsi,label = etiqueta)
plt.xlabel('ti')
plt.legend()
plt.title('Serie de Fourier')
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.

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

# espectro de frecuencia
k_i = np.arange(0,n,1,dtype=float)
ak_i = np.array(aki,dtype=float)
bk_i = np.array(bki,dtype=float)
ck_i = np.sqrt(ak_i**2 + bk_i**2)
cfs_i = np.arctan(-bk_i/ak_i)
print('ck: ', ck_i)
print('cfs:', cfs_i)
plt.subplot(211)
plt.stem(k_i,ck_i,label='Ck_magnitud')
plt.ylabel('ck_i')
plt.title('Especro de frecuencia')
plt.legend()
plt.grid()
plt.subplot(212)
plt.stem(k_i,cfs_i,label='Ck_fase')
plt.legend()
plt.ylabel('ck_fase')
plt.xlabel('k de kw')
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 integradas

Se incorporan mejoras en la presentación de resultados, al considerar valores casicero, dígitos a mostrar del los coeficientes, presentación de tablas de coeficientes.

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

# INGRESO
T = 2*np.pi

t = sym.Symbol('t')
ft = sym.Piecewise((-1,t <-T/2),
                   (-1,t <-T/4),
                   ( 1,t < T/4),
                   (-1,t < T/2),
                   (-1,True),)
# intervalo
a = -T/2
b = T/2

# número de coeficientes
n = 11
# Evaluación para gráfica
muestras = 101
digitos  = 4    # numpy set_print_options 
casicero = 1e-15

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

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

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

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

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

# evalua los coeficientes
a0  = ak.subs(k,0)/2
b0  = bk.subs(k,0)
aki = [a0]
bki = [b0]
n_i = [0]
i = 1
while not(i>=n):
    avalor = ak.subs(k,i)
    bvalor = bk.subs(k,i)
    aki.append(avalor)
    bki.append(bvalor)
    n_i.append(i)
    i = i+1
    
##print('\n coeficientes ak: ')
##print(aki)
##print('\n coeficientes bk: ')
##print(bki)



# tabla de valores ak, Bk
tabla = np.concatenate([[np.array(n_i,dtype = int)],
                        [np.array(aki,dtype = float)],
                        [np.array(bki,dtype = float)]],
                       axis=0)
tabla = np.transpose(tabla)
np.set_printoptions(precision=digitos)
print('\n coeficientes: \n [ k, \t\t ak, \t  bk ] ')
print(tabla)

# serie de Fourier expresión
serieF = a0 + 0*t 
i = 1
while not(i>=n):
    serieF = serieF + aki[i]*sym.cos(i*w0*t)
    serieF = serieF + bki[i]*sym.sin(i*w0*t)
    i = i+1
# serie = sym.simplify(serie)

print('\n serie de Fourier fs(t): ')
# sym.pprint(serieF)
print(a0)
for i in range(1,n,1):
    termino = '+ '
    if aki[i]<0:
        termino = '- '
    aki_str = str(abs(aki[i]))
    if abs(aki[i])<casicero:
        aki_str = '0'
    termino = termino + aki_str+' cos('+str(i*w0)+' t)'
    if bki[i]>=0:
        termino = termino + '+ '
    else:
        termino = termino + '- '
    bki_str = str(abs(bki[i]))
    if abs(bki[i])<casicero:
        bki_str = '0'
    termino = termino + bki_str+' sin('+str(i*w0)+' t)'
    print(termino)


# Grafica - Para evaluacion numerica
fsn = sym.lambdify(t,serieF)
ftn = sym.lambdify(t,ft)

# Evaluacion para grafica
ti = np.linspace(a,b,muestras)
fi = ftn(ti)
fsi = fsn(ti) 

# SALIDA
# Grafica f(t) y Fourier en t
figura, graf_ftT0 = plt.subplots()
graf_ftT0.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
graf_ftT0.plot(ti,fsi,label = etiqueta)
graf_ftT0.set_xlabel('t')
graf_ftT0.set_ylabel('f(t)')
graf_ftT0.legend()
graf_ftT0.grid()
graf_ftT0.set_title('Serie de Fourier de f(t)')
#plt.show()

# espectro de frecuencia Amplitud y fase
k_i = np.arange(0,n,1,dtype=float)
ak_i = np.array(aki,dtype=float)
bk_i = np.array(bki,dtype=float)

ck_i = np.sqrt(ak_i**2 + bk_i**2)
cfs_i = np.zeros(len(ak_i))
condicion = (abs(ak_i)>=casicero)
pendiente = -bk_i[condicion]/ak_i[condicion]
cfs_i[condicion] = np.arctan(pendiente)

tabla = np.concatenate([tabla,
                        np.transpose([ck_i]),
                        np.transpose([cfs_i])
                        ],
                       axis=1)
np.set_printoptions(precision=4)
print('\n coeficientes magnitud y fase: ')
print(' [k,\t\t ak,\t   bk,\t\t Ck,\t   fase_k ]')
print(tabla)

# grafica de espectro de frecuencia
figura, graf_Fw = plt.subplots(2,1)
graf_Fw[0].stem(k_i,ck_i,label='Ck_magnitud')
graf_Fw[0].set_ylabel('Ck')
graf_Fw[0].set_title('Espectro de frecuencia')
graf_Fw[0].legend()
graf_Fw[0].grid()

graf_Fw[1].stem(k_i,cfs_i,label='Ck_fase')
graf_Fw[1].legend()
graf_Fw[1].set_ylabel('Ck_fase')
graf_Fw[1].set_xlabel('k')
graf_Fw[1].grid()

plt.show()