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 tEjemplo
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 = 0Instrucciones 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()