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]