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 analizan las ventajas de utilizar las representaciones del dominio de frecuencia versus el dominio del tiempo de una señal y se ilustran los conceptos básicos mediante datos simulados y reales. El ejemplo responde a preguntas básicas como: ¿Cuál es el significado de la magnitud y la fase de una FFT? ¿Mi señal es periódica? ¿Cómo mido el poder? ¿Hay una, o más de una señal en esta banda?

El análisis del dominio de frecuencia es una herramienta de suma importancia en las aplicaciones de procesamiento de señales. El análisis del 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 del dominio del tiempo muestra cómo una señal cambia con el tiempo, el análisis del 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 cambio de fase que debe aplicarse 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 una transformación. Un ejemplo es la transformada 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 los componentes de frecuencia es la representación del dominio de frecuencia de la señal. La transformada inversa de Fourier convierte la función de dominio de frecuencia de nuevo en una función de tiempo. Las funciones de MATLAB le permiten calcular la transformada de Fourier discreta (DFT) de una señal y la inversa de esta transformación, respectivamente.fftifft

Información de magnitud y fase de la FFT

La representación del dominio de frecuencia de una señal lleva información sobre la magnitud y la fase de la señal en cada frecuencia. Esta es la razón por la salida de la computación FFT es compleja. Un número complejo, , tiene una parte real, , y una parte imaginaria, , de tal que . La magnitud de se calcula como , y la fase de se calcula como . Puede utilizar las funciones de MATLAB y, respectivamente, obtener la magnitud y la fase de cualquier número complejo.absangle

Utilice un ejemplo de audio para desarrollar una visión sobre qué información es transportada 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'); 

Se utiliza 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 de la FFT es un vector complejo que contiene información sobre el contenido de frecuencia de la señal. La magnitud le indica la fuerza de los componentes de frecuencia en relación con otros componentes. La fase le indica cómo todos los componentes de frecuencia se alinean en el tiempo.

Graficar la magnitud y los componentes de fase del espectro de frecuencias de la señal. La magnitud se traza convenientemente en una escala logarítmica (dB). La fase se desencapsula utilizando la función para que podamos ver una función continua de la 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 transformada inversa de Fourier al vector de dominio de frecuencia, Y, para recuperar la señal de tiempo. La marca ' simétrica ' indica que se está tratando con una señal de tiempo real-valorado por lo que será cero los pequeños componentes imaginarios que aparecen en la transformada inversa debido a inexactitudes 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 es también debido a las inexactitudes numéricas mencionadas anteriormente. Reproducir y escuchar 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, quite los componentes de frecuencia por encima de 1 kHz directamente de la salida FFT (haciendo las magnitudes 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') 

Vuelve a usar la señal filtrada en el dominio de tiempo.ifft

ylp = ifft(Ylp,'symmetric'); 

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

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

La fase de una señal tiene información importante acerca de Cuándo aparecen las notas de la canción en el tiempo. Para ilustrar la importancia de la fase en la señal de audio, quite la información de fase por completo tomando la magnitud de cada componente de frecuencia. Tenga en cuenta que al hacer esto, 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') 

Vuelve a obtener la señal en el dominio de tiempo y reproduce el audio. No se puede reconocer el sonido original en absoluto. La respuesta de magnitud es la misma, no se han quitado 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 representarla irreconocible.

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

Búsqueda de periodicidades de señal

La representación del 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 que no son visibles en absoluto cuando usted mira la señal en el dominio del tiempo. Por ejemplo, el análisis del dominio de frecuencia resulta útil cuando se busca el comportamiento cíclico de una señal.

Analyzing Cyclic Behavior of the Temperature in an Office Building

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. Mire 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 observando la señal del dominio del tiempo. Sin embargo, el comportamiento cíclico de la temperatura se hace evidente si miramos su representación del dominio de la frecuencia.

Obtenga la representación del 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 de temperatura controlada en un calendario de 7 días. La primera línea espectral indica que las temperaturas del edificio siguen un ciclo semanal con temperaturas más bajas durante los fines de semana y temperaturas más altas durante las semanas. 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]) 

Medir potencia

La función calcula la 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 horaria se distribuye con frecuencia, tiene unidades de vatios/Hz. Calcule el espectro de potencia integrando cada punto del PSD en el intervalo de frecuencias 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 integrarlos en 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.

Measuring Harmonics at the Output of a Non-Linear Power Amplifier

Cargue los datos medidos en la salida de un amplificador de potencia que tenga una distorsión de tercer orden de la forma Dónde es la tensión de salida y es la tensión de entrada. Los datos se capturaron con una frecuencia de muestreo de 3,6 kHz. La entrada consiste en una sinusoide 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 el segundo y tercer armónicos en 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 trama 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  

Improving Power Measurements for Noisy Signals

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 usted sólo analizó una realización corta de la señal ruidosa. 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 dividará en segmentos más pequeños de una longitud especificada, computará 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 consista en 500e3 puntos de la salida del amplificador. Mantenga el número de puntos usados para realizar las FFTs como 3600 de modo que el suelo (500e3/3600) = 138 FFTs se promedian 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 trama, elimina eficazmente todos los picos de frecuencia no esenciales causados por el ruido.pwelch El componente espectral a 180 Hz que fue enterrado en el ruido es ahora visible. El promedio elimina la varianza del espectro y esto produce mediciones de potencia más precisas de manera efectiva.

Measuring Total Average Power and Power Over a Frequency Band

Medir la potencia media total de una señal de dominio del 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 coincide 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 usted quiere medir la potencia total disponible sobre una banda de frecuencias? Puede utilizar la función para calcular la potencia sobre cualquier banda de frecuencias 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.

Calcule 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ómputo 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 el promedio 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  

Encontrar componentes espectrales

Una señal puede estar compuesta por uno o más componentes de frecuencia. La capacidad de observar todos los componentes espectrales depende de la resolución de frecuencia de su análisis. La resolución de frecuencia o el ancho de banda de resolución del espectro de potencia se define como R = 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.

Analyzing a Building's Earthquake Vibration Control System

Los sistemas de control de Active Mass driver (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, según 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 en el suelo. Las mediciones de aceleración se registraron en el primer piso de una estructura de prueba de tres pisos bajo condiciones sísmicas. Las mediciones se tomaron sin el sistema de control de masa activa (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. Utilice con segmentos de longitud 64 puntos de datos para obtener el suelo (10E3/64) = 156 FFT promedios y un ancho de banda de resolución de FS/64 = 15,625 Hz.pwelch Como se mostró antes, el promedio reduce los efectos de ruido y produce mediciones de potencia más precisas. Utilice 512 puntos FFT. El uso de NFFT > N interpolar eficazmente los puntos de frecuencia que renderizan 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 NFFT-Point FFT del vector de relleno cero).

Los espectros de potencia de aceleración de bucle abierto y de cierre 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 los datos de vibración y sabe que las vibraciones tienen un comportamiento cíclico. Entonces, ¿cómo es que las parcelas de espectro mostradas arriba no contienen ninguna línea espectral aguda típica del comportamiento cíclico? Tal vez usted está perdiendo esas líneas porque no son resolvable con la resolución obtenida con 64 longitudes de segmento de punto? Aumente la resolución de frecuencia para ver si hay líneas espectrales que no se resolvieron antes. Para ello, aumente la longitud del segmento de datos utilizada en la función a 512 puntos.pwelch Esto produce una nueva resolución de FS/512 = 1,9531 Hz. En este caso, el número de promedios FFT se reduce al suelo (10E3/512) = 19. Claramente, hay un intercambio 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 de la resolución de frecuencia le permite observar tres picos en el espectro de bucle abierto y dos en el espectro de bucle de cierre. Estos picos no eran resolvable 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 es de 5,86 Hz, y los picos de frecuencia equiespaciados sugieren que están relacionados armónicamente. Si bien 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 muesca el componente armónico a 17,58 Hz. Por lo que el sistema de control no sólo reduce la vibración, sino que también lo acerca a una sinusoide.

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. El aumento del número de puntos FFT interpolar los datos de frecuencia para darle más detalles sobre el espectro, pero no mejora la resolución.

Conclusiones

En este ejemplo aprendió a realizar el análisis del 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. Aprendió cómo calcular la potencia total o el poder sobre una banda particular de frecuencias de una señal ruidosa. Usted entendió cómo el aumento de 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.

Leer más

Para obtener más información sobre el análisis del dominio de frecuencia, consulte Signal Processing Toolbox.

Referencia: J. g. Manolakis, "procesamiento de señales digitales. Principios, algoritmos y aplicaciones ", Prentice Hall, 1996.

Apéndice

En este ejemplo se utilizan las siguientes funciones auxiliares.