Transformadas de Fourier
La transformada de Fourier es una fórmula matemática que transforma una señal muestreada en tiempo o espacio con la misma señal muestreada en frecuencia temporal o espacial. En el procesamiento de señales, la transformada de Fourier puede revelar características importantes de una señal, es decir, sus componentes de frecuencia.
La transformada de Fourier se define para un vector x con puntos n muestreados uniformemente por
ω=e−2πi/n es una de n raíces complejas de unidad donde i es la unidad imaginaria. Para x e y, los índices j y k oscilan entre 0 y n−1.
La función fft
de MATLAB® utiliza un algoritmo de la transformada rápida de Fourier para calcular la transformada de Fourier de los datos. Considere una señal sinusoidal x
que es una función de tiempo t
con componentes de frecuencia de 15 Hz y 20 Hz. Utilice un vector de tiempo muestreado en incrementos de 1/50 segundos en un periodo de 10 segundos.
Ts = 1/50; t = 0:Ts:10-Ts; x = sin(2*pi*15*t) + sin(2*pi*20*t); plot(t,x) xlabel('Time (seconds)') ylabel('Amplitude')
Calcule la transformada de Fourier de la señal y cree el vector f
que se corresponde con el muestreo de la señal en el espacio de frecuencia.
y = fft(x); fs = 1/Ts; f = (0:length(y)-1)*fs/length(y);
Cuando representa la magnitud de la señal como una función de frecuencia, los picos de la magnitud se corresponden con los componentes de frecuencia de 15 Hz y 20 Hz de la señal.
plot(f,abs(y)) xlabel('Frequency (Hz)') ylabel('Magnitude') title('Magnitude')
La transformada también genera una copia reflejada de los picos, que se corresponden con las frecuencias negativas de las señales. Para visualizar mejor esta periodicidad, puede utilizar la función fftshift
, que realiza un desplazamiento circular y centrado en cero en la transformada.
n = length(x); fshift = (-n/2:n/2-1)*(fs/n); yshift = fftshift(y); plot(fshift,abs(yshift)) xlabel('Frequency (Hz)') ylabel('Magnitude')
Señales sinusoidales
En aplicaciones científicas, las señales suelen corromperse debido a ruido aleatorio, lo que enmascara sus componentes de frecuencia. La transformada de Fourier puede procesar la separación del ruido aleatorio y revelar las frecuencias. Por ejemplo, cree una nueva señal, xnoise
, inyectando ruido gaussiano en la señal original, x
.
rng('default')
xnoise = x + 2.5*randn(size(t));
La potencia de la señal como una función de frecuencia es una métrica habitual utilizada en el procesamiento de señales. La potencia es la magnitud cuadrada de la transformada de Fourier de una señal, normalizada por el número de muestras de frecuencia. Calcule y represente el espectro de potencia de la señal de ruido centrada en la frecuencia cero. A pesar del ruido, podrá distinguir las frecuencias de la señal debido a los picos de potencia.
ynoise = fft(xnoise); ynoiseshift = fftshift(ynoise); power = abs(ynoiseshift).^2/n; plot(fshift,power) title('Power') xlabel('Frequency (Hz)') ylabel('Power')
Eficiencia computacional
Utilizar la fórmula de la transformada de Fourier directamente para calcular cada uno de los elementos n de y requiere n2 operaciones en punto flotante. El algoritmo de la transformada rápida de Fourier solo requiere n log n operaciones para calcularse. Esta eficiencia computacional supone una gran ventaja al procesar datos con millones de puntos de datos. Muchas implementaciones especializadas del algoritmo de la transformada rápida de Fourier son incluso más eficientes cuando n tiene factores primos pequeños, como n es una potencia de 2.
Considere los datos de audio recogidos por micrófonos submarinos en la costa de California. Estos datos se encuentran en una biblioteca del Programa de Investigación sobre Bioacústica de la Universidad de Cornell. Cargue y especifique el formato de un subconjunto de los datos de bluewhale.au
, que contiene la vocalización de una ballena azul del Pacífico. Puesto que los sonidos de las ballenas azules son sonidos de baja frecuencia, son casi imperceptibles para los humanos. La escala de tiempo en los datos se comprime con un factor de 10 para aumentar el tono y hacer el sonido más perceptible. Puede utilizar el comando sound(x,fs)
para escuchar el archivo de audio completo.
whaleFile = 'bluewhale.au'; [x,fs] = audioread(whaleFile); whaleMoan = x(2.45e4:3.10e4); t = 10*(0:1/fs:(length(whaleMoan)-1)/fs); plot(t,whaleMoan) xlabel('Time (seconds)') ylabel('Amplitude') xlim([0 t(end)])
Especifique una nueva longitud de señal que sea la siguiente potencia de 2 superior a la longitud original. Después, utilice fft
para calcular la transformada de Fourier utilizando la nueva longitud de señal. fft
rellena automáticamente los datos con ceros para aumentar el tamaño de la muestra. Este relleno puede agilizar de forma significativa el cálculo de la transformada, especialmente para tamaños de muestras con factores primos grandes.
m = length(whaleMoan); n = pow2(nextpow2(m)); y = fft(whaleMoan,n);
Represente el espectro de potencia de la señal. La gráfica indica que el gemido se compone de una frecuencia fundamental de aproximadamente 17 Hz y una secuencia de armónicos, donde se enfatiza el segundo armónico.
f = (0:n-1)*(fs/n)/10; % frequency vector power = abs(y).^2/n; % power spectrum plot(f(1:floor(n/2)),power(1:floor(n/2))) xlabel('Frequency (Hz)') ylabel('Power')
La fase de los sinusoides
Usando la transformada de Fourier, también puede extraer el espectro de fase de la señal original. Por ejemplo, cree una señal que conste de dos sinusoides con frecuencias de 15 Hz y 40 Hz. La primera sinusoide es una onda de coseno con fase −π/4 y la segunda es una onda de coseno con fase π/2. Muestree la señal a 100 Hz por 1 segundo.
fs = 100; t = 0:1/fs:1-1/fs; x = cos(2*pi*15*t - pi/4) - sin(2*pi*40*t);
Calcule la transformada de Fourier de la señal. Represente la magnitud de la transformada como una función de frecuencia.
y = fft(x); z = fftshift(y); ly = length(y); f = (-ly/2:ly/2-1)/ly*fs; stem(f,abs(z)) xlabel("Frequency (Hz)") ylabel("|y|") grid
Calcule la fase de la transformada eliminando los valores de la transformada de magnitudes pequeñas. Represente la fase como una función de frecuencia.
tol = 1e-6; z(abs(z) < tol) = 0; theta = angle(z); stem(f,theta/pi) xlabel("Frequency (Hz)") ylabel("Phase / \pi") grid
Consulte también
fft
| fftshift
| nextpow2
| ifft
| fft2
| fftn
| fftw
No hay comentarios:
Publicar un comentario