Estimaciones de densidad del espectro de potencia utilizando FFT
Este ejemplo muestra cómo obtener estimaciones de densidad del espectro de potencia (PSD) no paramétricas equivalentes al periodograma utilizando fft
. Los ejemplos muestran cómo adaptar correctamente el resultado de fft
para entradas de longitud par, para frecuencia normalizada y hercios, y para estimaciones PSD unilaterales y bilaterales.
Entrada de longitud par con tasa de muestreo
Obtenga el periodograma para una señal de longitud par muestreada a 1 kHz utilizando fft
y periodogram
. Compare los resultados.
Cree una señal compuesta por una onda sinusoidal de 100 Hz en N(0,1) ruido aditivo. La frecuencia de muestreo es 1 kHz. La longitud de la señal es de 1000 muestras. Utilice los valores predeterminados del generador de números aleatorios para obtener resultados reproducibles.
rng default
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));
Obtenga el periodograma utilizando fft
. La señal tiene valor real y longitud par. Puesto que la señal tiene valor real, solo necesita estimaciones de potencia para las frecuencias positivas o negativas. Para conservar la potencia total, multiplique todas las frecuencias en ambos conjuntos (las frecuencias positivas y negativas) por un factor de dos. La frecuencia cero (DC) y la frecuencia de Nyquist no se producen dos veces. Represente el resultado.
N = length(x); xdft = fft(x); xdft = xdft(1:N/2+1); psdx = (1/(Fs*N)) * abs(xdft).^2; psdx(2:end-1) = 2*psdx(2:end-1); freq = 0:Fs/length(x):Fs/2; plot(freq,10*log10(psdx)) grid on title('Periodogram Using FFT') xlabel('Frequency (Hz)') ylabel('Power/Frequency (dB/Hz)')
Calcule y represente el periodograma utilizando periodogram
. Muestre que los dos resultados son idénticos.
periodogram(x,rectwin(length(x)),length(x),Fs)
mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x),Fs))
mxerr = 3.4694e-18
Entrada con frecuencia normalizada
Utilice fft
para generar un periodograma para una entrada utilizando frecuencia normalizada. Cree una señal compuesta por una onda sinusoidal en N(0,1) ruido aditivo. La onda sinusoidal tiene una frecuencia angular de rad/muestra. Utilice los valores predeterminados del generador de números aleatorios para obtener resultados reproducibles.
rng default
n = 0:999;
x = cos(pi/4*n) + randn(size(n));
Obtenga el periodograma utilizando fft
. La señal tiene valor real y longitud par. Puesto que la señal tiene valor real, solo necesita estimaciones de potencia para las frecuencias positivas o negativas. Para conservar la potencia total, multiplique todas las frecuencias en ambos conjuntos (las frecuencias positivas y negativas) por un factor de dos. La frecuencia cero (DC) y la frecuencia de Nyquist no se producen dos veces. Represente el resultado.
N = length(x); xdft = fft(x); xdft = xdft(1:N/2+1); psdx = (1/(2*pi*N)) * abs(xdft).^2; psdx(2:end-1) = 2*psdx(2:end-1); freq = 0:(2*pi)/N:pi; plot(freq/pi,10*log10(psdx)) grid on title('Periodogram Using FFT') xlabel('Normalized Frequency (\times\pi rad/sample)') ylabel('Power/Frequency (dB/rad/sample)')
Calcule y represente el periodograma utilizando periodogram
. Muestre que los dos resultados son idénticos.
periodogram(x,rectwin(length(x)),length(x))
mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x)))
mxerr = 1.4211e-14
Entrada de valor complejo con frecuencia normalizada
Utilice fft
para generar un periodograma para una entrada de valor complejo con frecuencia normalizada. La señal es una exponencial compleja con una frecuencia angular de rad/muestra en ruido N(0,1) de valor complejo. Establezca la configuración predeterminada del generador de números aleatorios para obtener resultados reproducibles.
rng default
n = 0:999;
x = exp(1j*pi/4*n) + [1 1j]*randn(2,length(n))/sqrt(2);
Utilice fft
para obtener el periodograma. Puesto que la entrada es de valor complejo, obtenga el periodograma a partir de rad/muestra. Represente el resultado.
N = length(x); xdft = fft(x); psdx = (1/(2*pi*N)) * abs(xdft).^2; freq = 0:(2*pi)/N:2*pi-(2*pi)/N; plot(freq/pi,10*log10(psdx)) grid on title('Periodogram Using FFT') xlabel('Normalized Frequency (\times\pi rad/sample)') ylabel('Power/Frequency (dB/rad/sample)')
Utilice periodogram
para obtener y representar el periodograma. Compare las estimaciones de PSD.
periodogram(x,rectwin(length(x)),length(x),'twosided')
mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x),'twosided'))
mxerr = 2.8422e-14
Consulte también
Apps
Funciones
fft
|periodogram
|pspectrum