¿Cómo calculo numéricamente el desarrollo en serie de Fourier (DSF) de una función continua usando la FFT? F. J. Fraile Peláez En este documento explico cómo se usa la FFT para calcular el desarrollo en serie de Fourier (DSF) de una función continua (y periódica, claro). 1. ¿Qué es la FFT? Si quiere saberlo de verdad, búsquelo en un libro... A nuestros efectos, sólo necesitamos introducir la FFT de una forma muy breve y funcional. Concretamente: Considere un señal discreta periódica, ˜[] de  puntos. Al ser ˜[] una secuencia periódica, tiene un desarrollo en serie de Fourier que viene dado por ˜[] = −1 X 2      (1) =0 Pues bien, los  coeficientes  del desarrollo son precisamente la FFT de ˜[] Es decir, FFT [˜ []] = {0  2      −1 } 2. (2) FFT y Desarrollo en Serie de Fourier (DSF) de una función continua Veamos ahora cómo se relacionan los  de (2) (la FFT) con los coeficientes del desarrollo en serie de Fourier de una función periódica continua. Sea ˜() una función periódica de periodo  Su desarrollo en serie de Fourier es: ˜() = ∞ X  2     = =−∞ con ∞ X  0   (3) =−∞ 2  (4)  Supongamos que ˜() tiene un número finito de armónicos  (o que se puede limitar, con error despreciable, a un número finito de armónicos  ) O sea,  = 0 ||    Para representar a ˜() mediante sus muestras, debemos muestrear a más del doble de su ancho de banda; es decir, llamando  a la frecuencia de muestreo y 0 = 0 (2) 0 = 2   Supongamos que muestreamos a  veces el doble del ancho de banda. Es decir:   2( 0 ) = 1 (5) 2   ≥ 1  Llamando  = 1  la secuencia obtenida al muestrar ˜() será:  =  ¶ µ 2 ˜[] = ˜( ) =  exp    =− ¶ µ  X  2  exp   =  2  =− ¶ µ ¶ µ   X X 2 2 =  exp   =  exp    2   =− =− (6)  X  con (7)  (8)  = 2   Y por tanto   (9)  Obsérvese que  debe ser entero; de lo contrario, la secuencia no sería periódica y no podríamos hacer uso de (2). Puesto que  es entero,  debe ser entero o, al menos, un número fraccionario adecuado. Por ejemplo, si  = 7  podría ser 1 5; pero no 1 4  = La expresión (7) es casi análoga a (1). Intentemos hacerlas iguales. En primer lugar, el sumatorio de (7) se puede extender formalmente a  términos, sin más que considerar que los  sobrantes son todos cero: (−1)2 X ˜[] = =−(−1)2 ˜[] = 2 X =−2−1 µ ¶ 2  exp      impar µ ¶ 2  exp      par. (10) (11) ¡ ¢ ¡ ¢ Ahora bien, observamos que exp  2  = exp { ± } 2  Utilizaremos esta propiedad   para escribir (1) en forma idéntica a (10) o (11). Procedemos de la siguiente manera: Supongamos que  es impar. Dividamos en dos el sumatorio (1): ˜[] =  −1 X =0  2     = ( −1)2 X  2     =0 +  −1 X =(−1)2+1 2     ≡ 1 + 2  (12) ¡ ¢ ¡ ¢ 2 Pero haciendo uso de que exp  2  = exp { − }  , el segundo sumatorio se puede   escribir: 2 = −1 X 2  (−)    (13) =(−1)2+1 Llamemos 0 ≡  − en ( 13). Cuando  = ( −1)2+1 0 = −( −1)2 y cuando  =  −1 0 = −1 de manera que 2 −1 X 2 = 0 +  0 2   (14)  0 =−(−1)2 Insertando (14) en (12), queda (cambiando el nombre del índice de 2  0 → ): ˜[] = −1 X (−1)2  2   +  + X 2      (15) =0 =−(−1)2 Comparando (15) con (10), vemos que: ⎧ ⎪ ⎨    = ⎪ ⎩ +   −1 0≤≤ 2 (16)  −1 − ≤   0 2 Un resultado similar se obtiene cuando  es par (se deja como ejercicio para usted, “lector”). La conclusión es, por tanto, que los coeficientes  de la FFT de la secuencia discreta obtenida de muestrear un periodo de la función continua periódica son, convenientemente ordenados, los coeficientes  del desarrollo en serie de Fourier (DSF) de la función continua periódica. 3. Ejemplo Veamos un ejemplo. Sea la función ˜() = 3 cos3 () (17) El ancho de banda de ˜() es, claramente, 3 así que en este caso ˜() tiene exactamente 3 armónicos,  = 3. El número de muestras en un periodo ha de ser, por tanto, (18)  = 2  = 6  Tomemos, por ejemplo,  = 3 5 con lo que1  = 21 Para una ˜() tan sencilla como (17), es fácil obtener directamente la expresión analítica del desarrollo en serie de Fourier, lo cual nos servirá para comparar los resultados. Concretamente, 3 cos3 () = 0 375 −3 + 1 125 − + 1 125  + 0 375 3  (19) Obtengamos este resultado mediante la FFT. La secuencia periódica será 3 ˜[] = ˜( ) = ˜( ) = ˜ ( ) = 3 cos µ µ ¶ ¶ 2 2 3 3 = 3 cos  = 3 cos   µ   ¶ (20) Obsérvese que las muestras desde  = 0 hasta  =  − 1 cubren, efectivamente, un periodo completo de la secuencia periódica ˜[] [El periodo de la función continua ˜() es  = 2 –y no 2(3)–, a causa de la componente de frecuencia  –que aparece claramente en (19)]. El programa en MATLAB para calcular los coeficientes de Fourier sería, en este caso, el siguiente: 1 Como sabe, cuando hay que hacer cálculos intensivos es muy juicioso ajustar  a una potencia de 2 (8, 16, 32, 64, 128...), pues es así como el algoritmo FFT muestra su enorme eficiencia. 3 N = 21; x = 3*cos(2*pi/N*(0:N-1).^3; y0 = fft(x)/N; eje0 = (0:N-1); stem(eje0,abs(y0),’o’) figure y = fftshift(y0); % Estos son los coeficientes c_k eje = (-(N-1)/2:(N-1)/2); stem(eje,abs(y),’o’) La función fft proporciona el vector y0, cuyos elementos son los { } de (2). El resultado se divide por  porque la rutina fft de matlab difiere en este factor de la FFT tal como la hemos definido. La figura 1 muestra el vector y0. 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 5 10 15 20 Figura 1: FFT de ˜[]:Vector y0 formado por los   Para obtener los coeficientes  hay simplemente que reordenar los  de acuerdo con (16). Esto lo hace automáticamente la función fftshift de matlab. El resultado, vector y, se muestra en la figura 2. El eje de abcisas toma los valores del vector eje, definido para mostrar el número de armónico (es decir, ) desde −( − 1)2 = 10 hasta ( + 1)2 = 10 El resultado (eje, y) es: 4 -10.0000 -9.0000 -8.0000 -7.0000 -6.0000 -5.0000 -4.0000 -3.0000 -2.0000 -1.0000 0 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.3750 - 0.0000i 0.0000 + 0.0000i 1.1250 - 0.0000i 0.0000 1.1250 - 0.0000i 0.0000 - 0.0000i 0.3750 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i Claramente, aparecen los coeficientes correctos de las frecuencias  y 3; todos los demás son cero. 1.4 1.2 1 0.8 0.6 0.4 0.2 0 −10 −5 0 5 10 Figura 2: Coeficientes del desarrollo en serie de Fourier de ˜(): Vector y formado por los   Comparando las figuras 1 y 2 se observa cómo, efectivamente, la función fftshift ha reordenado los coeficientes de acuerdo con la fórmula (16): 5 0 1 2 3 = 0 = 0 = 1 = 1 125 = 2 = 0 = 3 = 0 375 −1 = −1+21 = 20 = 1 125 −2 = −2+21 = 19 = 0 −3 = −3+21 = 18 = 0 375 Todos los demás, 0 La elección de  = 21 en este ejemplo es arbitraria. Podíamos haber elegido un número menor; por ejemplo,  = 43 ⇒  = 8 En este caso se obtiene: -4.0000 -3.0000 -2.0000 -1.0000 0 1.0000 2.0000 3.0000 0.0000 0.3750 - 0.0000i 0.0000 + 0.0000i 1.1250 - 0.0000i 0.0000 1.1250 + 0.0000i 0.0000 - 0.0000i 0.3750 + 0.0000i, como era de esperar. c MMVI ° F. J. Fraile Peláez, Universidad de Vigo. Este documento se puede distribuir libremente sin modificaciones. Este documento ha sido escrito en LaTeX 2 mediante Scientific Workplace 5.5  http://www.com.uvigo.es/fot/jfraile/sci/ ↑ ↑↑ http://www.com.uvigo.es/fot/jfraile/ ↑↑↑ http://www.com.uvigo.es/fot/ 6