Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Introducción práctica al análisis de dominio de frecuencia

Este ejemplo muestra cómo realizar e interpretar el análisis básico de la señal del dominio de frecuencia. En el ejemplo se describen las ventajas de usar representaciones de dominio de frecuencia frente a dominio de tiempo de una señal e ilustra conceptos básicos mediante datos simulados y reales. El ejemplo responde a preguntas básicas tales como: ¿cuál es el significado de la magnitud y la fase de un FFT? ¿Mi señal es periódica? ¿Cómo mido la potencia? ¿Hay una o más de una señal en esta banda?

El análisis de dominio de frecuencia es una herramienta de suma importancia en las aplicaciones de procesamiento de señales. El análisis de dominio de frecuencia se utiliza ampliamente en áreas tales como comunicaciones, geología, teledetección y procesamiento de imágenes. Mientras que el análisis de dominio de tiempo muestra cómo cambia una señal con el tiempo, el análisis de dominio de frecuencia muestra cómo se distribuye la energía de la señal en un rango de frecuencias. Una representación de dominio de frecuencia también incluye información sobre el desplazamiento de fase que se debe aplicar a cada componente de frecuencia para recuperar la señal de tiempo original con una combinación de todos los componentes de frecuencia individuales.

Una señal se puede convertir entre los dominios de tiempo y frecuencia con un par de operadores matemáticos llamados transformación. Un ejemplo es la transformación de Fourier, que descompone una función en la suma de un número (potencialmente infinito) de componentes de frecuencia de onda sinusoidal. El "espectro" de componentes de frecuencia es la representación del dominio de frecuencia de la señal. La transformación inversa de Fourier convierte la función de dominio de frecuencia en una función de tiempo. Las funciones y de MATLAB le permiten calcular la transformación Discrete Fourier (DFT) de una señal y la inversa de esta transformación respectivamente.fftifft

Magnitud y Fase Información de la FFT

La representación de dominio de frecuencia de una señal contiene información sobre la magnitud y la fase de la señal en cada frecuencia. Esta es la razón por la que la salida del cálculo FFT es compleja. Un número complejo, $x$, tiene una parte real, $x_r$, y una parte imaginaria, $x_i$, de forma que $x = x_r + ix_i$. La magnitud de $x$ se calcula como $\sqrt{(x_r^2+x_i^2)}$, y la fase de $x$ se calcula como $\arctan{(x_i/x_r)}$. Puede utilizar funciones de MATLAB y obtener respectivamente la magnitud y la fase de cualquier número complejo.absangle

Utilice un ejemplo de audio para desarrollar información sobre qué información se lleva por la magnitud y la fase de una señal. Para ello, cargue un archivo de audio que contenga 15 segundos de música acústica de guitarra. La frecuencia de muestreo de la señal de audio es de 44,1 kHz.

Fs = 44100; y = audioread('guitartune.wav'); 

Utilícelo para observar el contenido de frecuencia de la señal.fft

NFFT = length(y); Y = fft(y,NFFT); F = ((0:1/NFFT:1-1/NFFT)*Fs).'; 

La salida del FFT es un vector complejo que contiene información sobre el contenido de frecuencia de la señal. La magnitud indica la intensidad de los componentes de frecuencia en relación con otros componentes. La fase le indica cómo se alinean todos los componentes de frecuencia en el tiempo.

Trazar la magnitud y los componentes de fase del espectro de frecuencia de la señal. La magnitud se traza convenientemente en una escala logarítmica (dB). La fase se desenvuelve utilizando la función para que podamos ver una función continua de frecuencia.unwrap

magnitudeY = abs(Y);        % Magnitude of the FFT phaseY = unwrap(angle(Y));  % Phase of the FFT  helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT) 

Puede aplicar una transformación inversa de Fourier al vector de dominio de frecuencia, Y, para recuperar la señal de tiempo. El indicador "simétrico" indica que está tratando con una señal de tiempo de valor real por lo que pondrá a cero los pequeños componentes imaginarios que aparecen en la transformación inversa debido a imprecisiones numéricas en los cálculos.ifft Observe que la señal de tiempo original, y, y la señal recuperada, y1, son prácticamente las mismas (la norma de su diferencia está en el orden de 1e-14). La diferencia muy pequeña entre los dos también se debe a las imprecisiones numéricas mencionadas anteriormente. Reproduce y escucha la señal no transformada y1.

y1 = ifft(Y,NFFT,'symmetric'); norm(y-y1) 
 ans =     3.7885e-14  
hplayer = audioplayer(y1, Fs); play(hplayer); 

Para ver los efectos de cambiar la respuesta de magnitud de la señal, elimine los componentes de frecuencia por encima de 1 kHz directamente de la salida FFT (haciendo que las magnitudes sean iguales a cero) y escuche el efecto que esto tiene en el sonido del archivo de audio. La eliminación de componentes de alta frecuencia de una señal se conoce como filtrado de paso bajo.

Ylp = Y; Ylp(F>=1000 & F<=Fs-1000) = 0;  helperFrequencyAnalysisPlot1(F,abs(Ylp),unwrap(angle(Ylp)),NFFT,...   'Frequency components above 1 kHz have been zeroed') 

Vuelva a obtener la señal filtrada en el dominio del tiempo mediante .ifft

ylp = ifft(Ylp,'symmetric'); 

Reproduce la señal. Todavía se puede escuchar la melodía, pero suena como si hubiera cubierto sus oídos (filtra sonidos de alta frecuencia cuando hace esto). A pesar de que las guitarras producen notas que están entre 400 y 1 kHz, a medida que toca una nota en una cuerda, la cuerda también vibra en múltiplos de la frecuencia base. Estos componentes de mayor frecuencia, conocidos como armónicos, son los que dan a la guitarra su tono particular. Cuando los quitas, haces que el sonido parezca "opaco".

hplayer = audioplayer(ylp, Fs); play(hplayer); 

La fase de una señal tiene información importante sobre cuándo aparecen las notas de la canción. Para ilustrar la importancia de la fase en la señal de audio, elimine la información de fase por completo tomando la magnitud de cada componente de frecuencia. Tenga en cuenta que al hacer esto se mantiene la respuesta de magnitud sin cambios.

% Take the magnitude of each FFT component of the signal Yzp = abs(Y); helperFrequencyAnalysisPlot1(F,abs(Yzp),unwrap(angle(Yzp)),NFFT,[],...   'Phase has been set to zero') 

Recupere la señal en el dominio de tiempo y reproduzca el audio. No se puede reconocer el sonido original en absoluto. La respuesta de magnitud es la misma, no se han eliminado componentes de frecuencia esta vez, pero el orden de las notas ha desaparecido por completo. La señal ahora consiste en un grupo de sinusoides todos alineados en el tiempo igual a cero. En general, las distorsiones de fase causadas por el filtrado pueden dañar una señal hasta el punto de hacerla irreconocible.

yzp = ifft(Yzp,'symmetric'); hplayer = audioplayer(yzp, Fs); play(hplayer); 

Búsqueda de periodicidades de señal

La representación de dominio de frecuencia de una señal le permite observar varias características de la señal que no son fáciles de ver, o no son visibles en absoluto cuando se mira la señal en el dominio de tiempo. Por ejemplo, el análisis de dominio de frecuencia se vuelve útil cuando se busca un comportamiento cíclico de una señal.

Análisis del comportamiento cíclico de la temperatura en un edificio de oficinas

Considere un conjunto de mediciones de temperatura en un edificio de oficinas durante la temporada de invierno. Las mediciones se tomaron cada 30 minutos durante aproximadamente 16,5 semanas. Observe los datos del dominio de tiempo con el eje de tiempo escalado a semanas. ¿Podría haber algún comportamiento periódico en estos datos?

load officetemp.mat Fs = 1/(60*30);   % Sample rate is 1 sample every 30 minutes t = (0:length(temp)-1)/Fs;  helperFrequencyAnalysisPlot2(t/(60*60*24*7),temp,...   'Time in weeks','Temperature (Fahrenheit)') 

Es casi imposible saber si hay algún comportamiento cíclico en las temperaturas de la oficina mirando la señal de dominio del tiempo. Sin embargo, el comportamiento cíclico de la temperatura se hace evidente si nos fijamos en su representación de dominio de frecuencia.

Obtenga la representación de dominio de frecuencia de la señal. Si traza la magnitud de la salida FFT con un eje de frecuencia escalado a ciclos/semana, puede ver que hay dos líneas espectrales que son claramente más grandes que cualquier otro componente de frecuencia. Una línea espectral se encuentra en 1 ciclo/semana, la otra se encuentra en 7 ciclos/semana. Esto tiene sentido dado que los datos provienen de un edificio con temperatura controlada en un calendario de 7 días. La primera línea espectral indica que las temperaturas de construcción siguen un ciclo semanal con temperaturas más bajas los fines de semana y temperaturas más altas durante la semana. La segunda línea indica que también hay un ciclo diario con temperaturas más bajas durante la noche y temperaturas más altas durante el día.

NFFT = length(temp);              % Number of FFT points F = (0 : 1/NFFT : 1/2-1/NFFT)*Fs; % Frequency vector  TEMP = fft(temp,NFFT); TEMP(1) = 0; % remove the DC component for better visualization  helperFrequencyAnalysisPlot2(F*60*60*24*7,abs(TEMP(1:NFFT/2)),...   'Frequency (cycles/week)','Magnitude',[],[],[0 10]) 

Potencia de medición

La función calcula el FFT de la señal y normaliza la salida para obtener una densidad espectral de potencia, PSD o un espectro de potencia desde el que se puede medir la potencia.periodogram El PSD describe cómo la potencia de una señal de tiempo se distribuye con frecuencia, tiene unidades de vatios/Hz. El espectro de potencia se calcula integrando cada punto del PSD durante el intervalo de frecuencia en el que se define ese punto (es decir, sobre el ancho de banda de resolución del PSD). Las unidades del espectro de potencia son vatios. Puede leer los valores de potencia directamente desde el espectro de potencia sin tener que integrarse a lo largo de un intervalo. Tenga en cuenta que el PSD y el espectro de potencia son reales, por lo que no contienen ninguna información de fase.

Medición de armónicos en la salida de un amplificador de potencia no lineal

Cargue los datos medidos en la salida de un amplificador de potencia que tenga distorsión de tercer orden del formulario $v_o = v_i + 0.75 v_i^2 + 0.5 v_i^3$Dónde $v_o$ es la tensión de salida y $v_i$ es la tensión de entrada. Los datos se capturaron con una frecuencia de muestreo de 3,6 kHz. La entrada $v_i$ consiste en un sinusoides de 60 Hz con amplitud de unidad. Debido a la naturaleza de la distorsión no lineal, debe esperar que la señal de salida del amplificador contenga un componente de CC, un componente de 60 Hz y armónicos segundo y tercero a 120 y 180 Hz.

Cargue 3600 muestras de la salida del amplificador, calcule el espectro de potencia y trace el resultado en una escala logarítmica (decibelios-vatios o dBW).

load ampoutput1.mat Fs = 3600; NFFT = length(y);  % Power spectrum is computed when you pass a 'power' flag input [P,F] = periodogram(y,[],NFFT,Fs,'power');  helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...   'Power spectrum (dBW)',[],[],[-0.5 200]) 

La gráfica del espectro de potencia muestra tres de los cuatro picos esperados en DC, 60 y 120 Hz. También muestra varios picos más espurios que deben ser causados por el ruido en la señal. Tenga en cuenta que el armónico de 180 Hz está completamente enterrado en el ruido.

Mida la potencia de los picos esperados visibles:

PdBW = 10*log10(P); power_at_DC_dBW = PdBW(F==0)   % dBW  [peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,'minpeakheight',-11); peakFreqs_Hz = F(peakFreqIdx) peakPowers_dBW 
 power_at_DC_dBW =     -7.8873   peakFreqs_Hz =      60    120   peakPowers_dBW =     -0.3175   -10.2547  

Mejora de las mediciones de potencia para señales de noventa

Como se ve en la trama anterior, el periodograma muestra varios picos de frecuencia que no están relacionados con la señal de interés. El espectro se ve muy ruidoso. La razón de esto es que sólo analizó una breve realización de la señal de ruidoso. Repetir el experimento varias veces y promediar eliminaría los picos espectrales espurios y produciría mediciones de potencia más precisas. Puede lograr este promedio utilizando la función.pwelch Esta función tomará un vector de datos grande, lo dividirá en segmentos más pequeños de una longitud especificada, calculará tantos periodogramas como segmentos y los promediará. A medida que aumenta el número de segmentos disponibles, la función producirá un espectro de potencia más suave (menos varianza) con valores de potencia más cercanos a los valores esperados.pwelch

Cargue una observación más grande que consta de 500e3 puntos de la salida del amplificador. Mantenga el número de puntos utilizados para realizar los FFT como 3600 de modo que el piso(500e3/3600) a 138 FFT se promedia para obtener el espectro de potencia.

load ampoutput2.mat SegmentLength = NFFT;  % Power spectrum is computed when you pass a 'power' flag input [P,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,'power');  helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...   'Power spectrum (dBW)',[],[],[-0.5 200]) 

Como se ve en la parcela, elimina eficazmente todos los picos de frecuencia espurios causados por el ruido.pwelch El componente espectral a 180 Hz que estaba enterrado en el ruido es ahora visible. El promedio elimina la varianza del espectro y esto produce efectivamente mediciones de potencia más precisas.

Medición de la potencia media total y la potencia sobre una banda de frecuencia

Medir la potencia media total de una señal de dominio de tiempo es una tarea fácil y común. Para la señal de salida del amplificador, y, la potencia media total se calcula en el dominio de tiempo como:

pwr = sum(y.^2)/length(y) % in watts 
 pwr =      8.1697  

En el dominio de frecuencia, la potencia media total se calcula como la suma de la potencia de todos los componentes de frecuencia de la señal. El valor de pwr1 consiste en la suma de todos los componentes de frecuencia disponibles en el espectro de potencia de la señal. El valor está de acuerdo con el valor de pwr calculado anteriormente utilizando la señal de dominio de tiempo:

pwr1 = sum(P) % in watts 
 pwr1 =      8.1698  

Pero, ¿qué pasa si quieres medir la potencia total disponible sobre una banda de frecuencias? Puede utilizar la función para calcular la potencia en cualquier banda de frecuencia deseada.bandpower Puede pasar la señal de dominio de tiempo directamente como entrada a esta función para obtener la potencia sobre una banda especificada. En este caso, la función estimará el espectro de potencia con el método de periodograma.

Calcular la potencia sobre la banda de 50 Hz a 70 Hz. El resultado incluirá la potencia de 60 Hz más la potencia de ruido sobre la banda de interés:

pwr_band = bandpower(y,Fs,[50 70]); pwr_band_dBW = 10*log10(pwr_band) % dBW 
 pwr_band_dBW =      0.0341  

Si desea controlar el cálculo del espectro de potencia utilizado para medir la potencia en una banda, puede pasar un vector PSD a la función.bandpower Por ejemplo, puede utilizar la función como lo hizo antes para calcular el PSD y garantizar la promediación de los efectos de ruido:pwelch

% Power spectral density is computed when you specify the 'psd' option [PSD,F]  = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,'psd'); pwr_band1 = bandpower(PSD,F,[50 70],'psd'); pwr_band_dBW1 = 10*log10(pwr_band1) % dBW 
 pwr_band_dBW1 =      0.0798  

Búsqueda de componentes espectrales

Una señal puede estar compuesta de uno o más componentes de frecuencia. La capacidad de observar todos los componentes espectrales depende de la resolución de frecuencia del análisis. La resolución de frecuencia o el ancho de banda de resolución del espectro de potencia se define como R a Fs/N, donde N es la longitud de la observación de la señal. Solo se resolverán los componentes espectrales separados por una frecuencia mayor que la resolución de frecuencia.

Análisis del sistema de control de vibraciones sísmicas de un edificio

Los sistemas de control de controlador de masa activa (AMD) se utilizan para reducir la vibración en un edificio bajo un terremoto. Un conductor de masa activo se coloca en la planta superior del edificio y, en función de las mediciones de desplazamiento y aceleración de los pisos del edificio, un sistema de control envía señales al conductor para que la masa se mueva para atenuar las perturbaciones del suelo. Las mediciones de aceleración se registraron en el primer piso de una estructura de prueba de tres pisos en condiciones de terremoto. Las mediciones se tomaron sin el sistema de control de controlador de masa activo (condición de bucle abierto) y con el sistema de control activo (condición de bucle cerrado).

Cargue los datos de aceleración y calcule el espectro de potencia para la aceleración del primer piso. La longitud de los vectores de datos es 10e3 y la frecuencia de muestreo es de 1 kHz. Se utiliza con segmentos de longitud 64 puntos de datos para obtener puntos de suelo(10e3/64) a 156 promedios FFT y un ancho de banda de resolución de Fs/64 a 15.625 Hz.pwelch Como se ha demostrado anteriormente, el promedio reduce los efectos de ruido y produce mediciones de potencia más precisas. Utilice 512 puntos FFT. El uso de NFFT > N interpola eficazmente puntos de frecuencia que representan una gráfica de espectro más detallada (esto se logra añadiendo ceros NFFT-N al final de la señal de tiempo y tomando el punto NFFT FFT del vector acolchado cero).

Los espectros de potencia de aceleración de bucle abierto y bucle cerrado muestran que cuando el sistema de control está activo, el espectro de potencia de aceleración disminuye entre 4 y 11 dB. La atenuación máxima se produce a unos 23,44 kHz. Una reducción de 11 dB significa que la potencia de vibración se reduce en un factor de 12,6. La potencia total se reduce de 0,1670 a 0,059 vatios, un factor de 2,83.

load quakevibration.mat  Fs = 1e3;                 % sample rate NFFT = 512;               % number of FFT points segmentLength = 64;       % segment length  % open loop acceleration power spectrum [P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power');  % closed loop acceleration power spectrum P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');  helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...   'Frequency in Hz','Acceleration Power Spectrum in dB',...   'Resolution bandwidth = 15.625 Hz',{'Open loop', 'Closed loop'},[0 100]) 

Está analizando datos de vibración y sabe que las vibraciones tienen un comportamiento cíclico. Entonces, ¿cómo es que las gráficas de espectro mostradas anteriormente no contienen líneas espectrales nítidas típicas del comportamiento cíclico? Tal vez le faltan esas líneas porque no se pueden resolver con la resolución obtenida con longitudes de segmento de 64 puntos? Aumente la resolución de frecuencia para ver si hay líneas espectrales que antes no se pudieran resolver. Para ello, aumente la longitud del segmento de datos utilizado en la función a 512 puntos.pwelch Esto produce una nueva resolución de Fs/512 a 1.9531 Hz. En este caso, el número de promedios ffT se reduce a suelo(10e3/512) a 19. Claramente, hay un equilibrio entre el número de promedios y la resolución de frecuencia cuando se utiliza .pwelch Mantenga el número de puntos FFT igual a 512.

NFFT = 512;            % number of FFT points segmentLength = 512;   % segment length  [P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power'); P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');  helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...   'Frequency in Hz','Acceleration Power Spectrum in dB',...   'Resolution bandwidth = 1.95 Hz',{'Open loop', 'Closed loop'},[0 100]) 

Observe cómo el aumento en la resolución de frecuencia le permite observar tres picos en el espectro de bucle abierto y dos en el espectro de bucle cerrado. Estos picos no se resolvían antes. La separación entre los picos en el espectro de bucle abierto es de aproximadamente 11 Hz, que es menor que la resolución de frecuencia obtenida con segmentos de longitud 64 pero mayor que la resolución obtenida con segmentos de longitud 512. El comportamiento cíclico de las vibraciones es ahora visible. La frecuencia de vibración principal está a 5,86 Hz, y los picos de frecuencia equiespacidos sugieren que están relacionados armónicamente. Aunque ya se ha observado que el sistema de control reduce la potencia general de las vibraciones, los espectros de mayor resolución muestran que otro efecto del sistema de control es entallar el componente armónico a 17,58 Hz. Así que el sistema de control no sólo reduce la vibración, sino que también la acerca a un sinusoides.

Es importante tener en cuenta que la resolución de frecuencia está determinada por el número de puntos de señal, no por el número de puntos FFT. Al aumentar el número de puntos FFT, se interpolan los datos de frecuencia para darle más detalles sobre el espectro, pero no mejora la resolución.

Conclusiones

En este ejemplo ha aprendido a realizar análisis de dominio de frecuencia de una señal mediante las funciones , , , . y .fftifftperiodogrampwelchbandpower Usted entendió la naturaleza compleja de la FFT y cuál es la información contenida en la magnitud y la fase del espectro de frecuencias. Ha visto las ventajas de utilizar datos de dominio de frecuencia al analizar la periodicidad de una señal. Aprendiste cómo calcular la potencia total o la potencia sobre una banda particular de frecuencias de una señal ruidoso. Entendió cómo aumentar la resolución de frecuencia del espectro le permite observar componentes de frecuencia estrechamente espaciados y aprendió sobre el equilibrio entre la resolución de frecuencia y el promedio espectral.

Lectura adicional

Para obtener más información sobre el análisis de dominio de frecuencia, consulte el cuadro de herramientas de procesamiento de señales.

Referencia: J.G. Proakis y D. G. Manolakis, "Procesamiento de señaldigital. Principios, algoritmos y aplicaciones", Prentice Hall, 1996.

Apéndice

En este ejemplo se utilizan las siguientes funciones auxiliares.