Transformada de Fourier discreta en Python con SciPy

Introducción

En este artículo vamos a ver cómo calcular la transformada de Fourier discreta (o DFT) de una señal en Python utilizando la transformada rápida de Fourier (o FFT) implementada en SciPy. El análisis de Fourier es la herramienta fundamental en procesamiento de señales y resulta útil en otras áreas como en la resolución de ecuaciones diferenciales o en el tratamiento de imágenes.

Nota: Puedes ver el artículo en forma de notebook de IPython mediante la herramienta nbviewer.

La transformada rápida de Fourier o FFT es en realidad una familia de algoritmos muy eficientes (\( O(n \log n)\)) para calcular la DFT de una señal discreta, de los cuales el más utilizado es el algoritmo de Cooley-Tukey. Este es el que está implementado en SciPy a través de las subrutinas FFTPACK, escritas en FORTRAN 77, y es el que vamos a utilizar. También podríamos haber utilizado la biblioteca FFTW, más moderna y escrita en C, o la implementación presente en NumPy, que es una traducción a C de FFTPACK y que funciona igual que la de SciPy, pero no lo vamos a hacer.

Nótese que estos métodos no nos permiten calcular ni la serie de Fourier de una función periódica ni la transformada de Fourier de una función no periódica.

Estas operaciones forman parte del cálculo simbólico y deben llevarse a cabo con otro tipo de programas, como SymPy o Sage. En Pybonacci puedes leer una introducción a SymPy y cómo calcular series con SymPy, así como una reseña sobre Sage.

En este artículo nos vamos a centrar en el caso unidimensional, aunque el código es fácilmente extrapolable al caso de dos y más dimensiones.

Definiciones

Antes de empezar de lleno con las FFTs vamos a hacer una brevísima introducción matemática. He intentado que esto no se convierta en una larga introducción al procesamiento de señales digitales, pero si te interesa el tema y quieres más rigor puedes acudir a cualquiera de las referencias que doy más abajo además de la documentación de scipy.fftpack.

En este artículo vamos a manejar señales discretas reales, bien creadas por nosotros o procedentes de muestrear una señal analógica a una frecuencia de muestreo fs. La transformada discreta de Fourier o DFT de la señal de entrada, que diremos que está en el dominio del tiempo, será otra señal discreta, n general compleja, cuyos coeficientes vienen definidos en SciPy por la expresión

\( \displaystyle A_k = \sum_{m=0}^{n-1} a_m \exp\left(-2\pi j{mk \over n}\right) \qquad k = 0,\ldots,n-1,\)

siendo \( a_n\) los coeficientes de la señal de entrada y \( A_k\) los de la señal obtenida. Diremos que \( A_k\) está en el dominio de la frecuencia y que es el espectro o transformada de \( a_n\). La resolución en frecuencia será mayor cuanto mayor sea el número de intervalos \( n\) de la señal de entrada. La DFT se calcula con el algoritmo de la transformada rápida de Fourier o FFT, que es más eficiente si \( n\) es una potencia entera de 2.

Cada uno de los componentes de frecuencia se puede expresar en forma compleja como \( A \exp(j \omega t)\). El espectro complejo de una función seno, por ejemplo, se obtiene directamente de su expresión en forma compleja

\( \displaystyle A \sin(\omega t) = A \frac{e^{j \omega t} – e^{-j \omega t}}{2 j} = -j \frac{A}{2} e^{j \omega t} + j \frac{A}{2} e^{-j \omega t}\)

y tendrá, por tanto, una raya espectral a la frecuencia \( \omega\) y otra a la frecuencia \( -\omega\), como ya comprobaremos.

La DFT de una señal siempre asume que la señal de entrada es un período de una secuencia periódica. Esto será importante por lo que luego veremos. Sería, por tanto, el equivalente en tiempo discreto a las series de Fourier.

FFT de una función senoidal sencilla

Vamos a calcular la FFT de una función senoidal en Python utilizando SciPy. Para que sea verdaderamente sencilla tenemos que preparar un poco los datos, o si no acabaremos observando efectos adversos. Estos los estudiaremos más adelante.

Para manejar números redondos, vamos a recrear una señal con un armónico de \( f = 400~Hz\) y otro del doble de frecuencia:

\( \displaystyle y(t) = \sin(2 \pi f) – \frac{1}{2} \sin(2 \pi \cdot 2 f)\)

donde \( 2 \pi f = \omega\). Nótese que, por lo que hemos visto en el apartado anterior, veremos cuatro rayas en el espectro de esta función: dos para la frecuencia \( f\) y otras dos para la \( 2 f\). Vamos a imponer a priori el número de intervalos y la distancia entre ellos y de ahí vamos a calcular el intervalo de tiempo. Evidentemente así no se funciona cuando muestreamos un archivo de audio, pongo por caso, pero como digo así obtendremos el resultado exacto. El código y la salida serían las siguientes:

import matplotlib.pyplot as plt
import numpy as np
from numpy import pi

n = 2 ** 6 # Número de intervalos
f = 400.0 # Hz
dt = 1 / (f * 16) # Espaciado, 16 puntos por período
t = np.linspace(0, (n - 1) * dt, n) # Intervalo de tiempo en segundos
y = np.sin(2 * pi * f * t) - 0.5 * np.sin(2 * pi * 2 * f * t) # Señal
plt.plot(t, y)
plt.plot(t, y, 'ko')
plt.xlabel('Tiempo (s)')
plt.ylabel('$y(t)$')

Esta es la señal que vamos a transformar. Fijáos en el último punto representado. Como el intervalo va desde 0 hasta n - 1, los trozos «empalman» perfectamente. Ahora vamos a hallar la DFT de la señal y. Para ello necesitamos dos funciones, que se importan del paquete scipy.fftpack:

• La función fft, que devuelve la DFT de la señal de entrada. Si queremos que esté normalizada, tenemos que dividir después por el número de intervalos N.
• La función fftfreq, que devuelve las frecuencias a las que corresponde cada punto de la DFT. Se usa conjuntamente con la función anterior y acepta dos argumentos: el número de elementos n y el espacio entre intervalos dt, que ya hemos calculado antes.

Como en este caso solo hay dos frecuencias fundamentales, vamos a representarlas utilizando la función vlines de matplotlib. El código y la salida quedarían así:

from scipy.fftpack import fft, fftfreq

Y = fft(y) / n # Normalizada
frq = fftfreq(n, dt) # Recuperamos las frecuencias
plt.vlines(frq, 0, Y.imag) # Representamos la parte imaginaria
plt.annotate(s=u'f = 400 Hz', xy=(400.0, -0.5), xytext=(400.0 + 1000.0, -0.5 - 0.35), arrowprops=dict(arrowstyle = "->"))
plt.annotate(s=u'f = -400 Hz', xy=(-400.0, 0.5), xytext=(-400.0 - 2000.0, 0.5 + 0.15), arrowprops=dict(arrowstyle = "->"))
plt.annotate(s=u'f = 800 Hz', xy=(800.0, 0.25), xytext=(800.0 + 600.0, 0.25 + 0.35), arrowprops=dict(arrowstyle = "->"))
plt.annotate(s=u'f = -800 Hz', xy=(-800.0, -0.25), xytext=(-800.0 - 1000.0, -0.25 - 0.35), arrowprops=dict(arrowstyle = "->"))
plt.ylim(-1, 1)
plt.xlabel('Frecuencia (Hz)')
plt.ylabel('Im($Y$)')

Nota: Si quieres leer más sobre añadir texto y anotaciones a una gráfica en matplotlib, puedes leer la parte VIII del nuestro magnífico tutorial de matplotlib por Kiko 🙂

¡Y ya hemos hecho nuestra primera transformada discreta de Fourier con el algoritmo FFT! Pero las cosas en la realidad son bastante más complicadas. Vamos a ver algunos de los problemas que aparecerán con frecuencia.

Fuga, «zero-padding» y funciones ventana

Vamos a analizar el ejemplo anterior, pero en esta ocasión vamos a poner menos cuidado con la preparación de los datos. Veamos qué ocurre:

n2 = 2 ** 5
t2 = np.linspace(0, 0.012, n2) # Intervalo de tiempo en segundos
dt2 = t2[1] - t2[0]
y2 = np.sin(2 * pi * f * t2) - 0.5 * np.sin(2 * pi * 2 * f * t2) # Señal
Y2 = fft(y2) / n2 # Transformada normalizada
frq2 = fftfreq(n2, dt2)
fig = plt.figure(figsize=(6, 8))
ax1 = fig.add_subplot(211)
ax1.plot(t2, y2)
ax1.set_xlabel('Tiempo (s)')
ax1.set_ylabel('$y_2(t)$')
ax2 = fig.add_subplot(212)
ax2.vlines(frq2, 0, Y2.imag)
plt.xlabel('Frecuencia (Hz)')
plt.ylabel('Im($Y_2$)')

Como se puede ver, el fenómeno no es extremadamente importante pero han aparecido otras rayas espectrales que no esperábamos. Esto se conoce como «leaking» (y yo lo voy a traducir por fuga) y es debido a que, en este caso, los trozos «no empalman exactamente». Recuerda que la DFT, y por extensión la FFT asume que estamos transformando un período de una señal periódica. Si utilizamos más puntos y extendemos la señal con ceros (esto se conoce como «zero-padding») la DFT da más resolución en frecuencia pero la fuga se magnifica:

t3 = np.linspace(0, 0.012 + 9 * dt2, 10 * n2) # Intervalo de tiempo en segundos
y3 = np.append(y2, np.zeros(9 * n2)) # Señal
Y3 = fft(y3) / (10 * n2) # Transformada normalizada
frq3 = fftfreq(10 * n2, dt2)
fig = plt.figure(figsize=(6, 8))
ax1 = fig.add_subplot(211)
ax1.plot(t3, y3)
ax1.set_xlabel('Tiempo (s)')
ax1.set_ylabel('$y_3(t)$')
ax2 = fig.add_subplot(212)
ax2.vlines(frq3, 0, Y3.imag)
plt.xlabel('Frecuencia (Hz)')
plt.ylabel('Im($Y_3$)')

Existe una manera de reducir la fuga y es mediante el uso de funciones ventana. Las funciones ventana no son más que funciones que valen cero fuera de un cierto intervalo, y que en procesamiento de señales digitales se utilizan para «suavizar» o filtrar una determinada señal. NumPy trae unas cuantas funciones ventana por defecto; por ejemplo, la ventana de Blackman tiene este aspecto:

Como se puede ver, en los extremos del intervalo es nula. Las funciones ventana reciben un único argumento que es el número de puntos. Si multiplicamos la ventana por la señal, obtenemos una nueva señal que vale cero en los extremos. Comprobemos el resultado, representando ahora el espectro de amplitud y comparando cómo es el resultado si aplicamos o no la ventana de Blackman:

n4 = 2 ** 8
t4 = np.linspace(0, 0.05, n4)
dt4 = t4[1] - t4[0]
y4 = np.sin(2 * pi * f * t4) - 0.5 * np.sin(2 * pi * 2 * f * t4)
y5 = y4 * np.blackman(n4)
t4 = np.linspace(0, 0.12 + 4 * dt4, 5 * n4)
y4 = np.append(y4, np.zeros(4 * n4))
y5 = np.append(y5, np.zeros(4 * n4))
Y4 = fft(y4) / (5 * n4)
Y5 = fft(y5) / (5 * n4)
frq4 = fftfreq(5 * n4, dt4)
fig = plt.figure(figsize=(6, 8))
ax1 = fig.add_subplot(411)
ax1.plot(t4, y4)
plt.xlabel('Frecuencia (Hz)')
plt.ylabel('$y_4(t)$')
ax2 = fig.add_subplot(412)
ax2.vlines(frq4, 0, abs(Y4)) # Espectro de amplitud
plt.xlabel('Frecuencia (Hz)')
plt.ylabel('Abs($Y_4$)')
ax3 = fig.add_subplot(413)
ax3.plot(t4, y5)
plt.xlabel('Frecuencia (Hz)')
plt.ylabel('$y_5(t)$')
ax4 = fig.add_subplot(414)
ax4.vlines(frq4, 0, abs(Y5)) # Espectro de amplitud
plt.xlabel('Frecuencia (Hz)')
plt.ylabel('Abs($Y_5$)')

Esto ya ha sido más interesante, ¿no? 🙂

Otros ejemplos: procesamiento de audio, espectrogramas

AVISO (KIKO): EL CÓDIGO EN ESTA PARTE FINAL NO FUNCIONA EN NUEVAS VERSIONES DE PYTHON DEBIDO A QUE SCIKIT.AUDIOLAB NO SE ACTUALIZA DESDE 2010 PERO HE QUERIDO DEJARLO TAL CUAL POR MOTIVOS HISTÓRICOS.

Para acabar de manera triunfal, vamos a ver cómo representar el espectrograma de un fragmento de audio, utilizando la función specgram de matplotlib. El espectrograma se obtiene de aplicar la transformada corta de Fourier (STFT) a una señal, que es como aplicar una DFT a pequeños fragmentos de la misma, para ver la evolución de las frecuencias a lo largo del tiempo. Para ello, vamos a usar este fragmento de una melodía tocada con un violín, disponible libremente en Wikimedia Commons y que está muestreado a 44,1 kHz:

Al tratarse de un archivo OGG, vamos a utilizar el SciKit audiolab.

import matplotlib.pyplot as plt
import scikits.audiolab as audiolab
sound = audiolab.sndfile('Violin_for_spectrogram.ogg', 'read')
y = sound.read_frames(sound.get_nframes())
Pxx, freqs, bins, im = plt.specgram(y, NFFT=512, Fs=44100)
plt.xlim(0, len(y) / 44100.0)
plt.ylim(0, 22050.0)
plt.colorbar(im).set_label(u'Intensidad (dB)')
plt.xlabel(u'Tiempo (s)')
plt.ylabel(u'Frecuencia (Hz)')

Se puede apreciar cómo va cambiando la frecuencia fundamental (barras horizontales inferiores) y los armónicos (todas las que están encima, cada vez más débiles) a medida que el intérprete toca una nota distinta. En definitiva, una preciosidad 😛

Espero que te haya gustado el artículo. Más abajo tienes las referencias en las que me he basado y unos cuantos enlaces que me han sido muy útiles para aprender procesamiento de señales digitales, porque cuando empecé a escribir el artículo no sabía nada sobre el tema. Sobre todo he aprendido la teoría con [II] y los ejemplos están inspirados en el capítulo «Example Applications of the DFT» de [I]. Si ves algún error o quieres apuntar alguna sugerencia, no dudes en decirlo en los comentarios.

¡Un saludo desde Milano! 🙂

Referencias

1. SMITH, Julius O. Mathematics of the discrete Fourier transform (DFT) with audio applications [en línea]. Segunda edición. Disponible en Web: <https://ccrma.stanford.edu/~jos/mdft/mdft.html>. [Consulta: 28 de septiembre de 2012]
2. LOPEZ-LEZCANO, Fernando; CACERES, Juan-Pablo. Introducción al Procesamiento Digital de Señales y a la Música por Computador [en línea]. Disponible en Web: <https://ccrma.stanford.edu/workshops/cm2007/index.shtml>. [Consulta: 28 de septiembre de 2012]

Enlaces

• Análisis espectral de una sinusoide https://ccrma.stanford.edu/~jos/mdft/Spectrum_Analysis_Sinusoid_Windowing.html
• Teoría general http://paulbourke.net/miscellaneous/dft/
• NumPy no es adecuado para series de Fourier http://stackoverflow.com/questions/4258106/how-to-calculate-a-fourier-series-in-numpy
• Reconstrucción de la señal a partir de algunos términos calculados con fft http://stackoverflow.com/questions/5843085/fourier-series-in-numpy-question-about-previous-answer
• Uso de fftfreq http://stackoverflow.com/questions/9456037/scipy-numpy-fft-frequency-analysis
• Unidades de fftfreq http://stackoverflow.com/questions/8573702/units-of-frequency-when-using-fft-in-numpy
• Frecuencias de fftfreq http://stackoverflow.com/questions/3694918/how-to-extract-frequency-associated-with-fft-values-in-python
• Artículo de The Glowing Python http://glowingpython.blogspot.com.es/2011/08/how-to-plot-frequency-spectrum-with.html
• Demo de MathWorks http://www.mathworks.es/products/matlab/examples.html?file=/products/demos/shipping/matlab/fftdemo.html
• Documentación de MATLAB http://www.mathworks.es/help/matlab/ref/fft.html
• Diferencias entre NumPy y SciPy http://stackoverflow.com/questions/6363154/what-is-the-difference-between-numpy-fft-and-scipy-fftpack
• Zero-padding of a signal http://dsp.stackexchange.com/questions/741/why-should-i-zero-pad-a-signal-before-taking-the-fourier-transform
• Uso de funciones ventana https://docs.google.com/open?id=1FliaqX0ModFkrPEvCtconY4vNlzU0r3luLkSl9Zk3WB48YyyAd12oD-Pw9IQ
• Para ampliar: DSP con SciPy http://docs.scipy.org/doc/scipy/reference/signal.html
• Python para procesamiento de audio http://lac.linuxaudio.org/2011/papers/40.pdf
• Definición de intervalos de frecuencia http://www.maximintegrated.com/glossary/definitions.mvp/term/frequency_bin/gpk/136