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