Estimations de la densité spectrale de puissance avec la FFT

Cet exemple montre comment obtenir des estimations non paramétriques équivalentes de la PSD (densité spectrale de puissance) avec les fonctions periodogram et fft. Les différents cas indiquent comment mettre à l’échelle la sortie de fft pour des entrées de longueur paire, pour des fréquences normalisées et en hertz, ainsi que pour des estimations de la PSD unilatérale et bilatérale. Une fenêtre rectangulaire est utilisée dans tous les cas.

Cet exemple porte sur des signaux stationnaires dont le contenu fréquentiel est constant dans le temps. Pour les signaux non stationnaires dont le contenu fréquentiel est dépendant du temps, la transformée de Fourier à court terme est un meilleur outil d’analyse. Pour plus d’informations, veuillez consulter Spectrogram Computation with Signal Processing Toolbox.

Entrée de longueur paire avec fréquence d’échantillonnage

Obtenez le périodogramme d’un signal de longueur paire échantillonné à 1 kHz en utilisant les fonctions fft et periodogram. Comparez les résultats.

Créez un signal composé d’une onde sinusoïdale de 100 Hz et d’un bruit additif N(0,1). La fréquence d’échantillonnage est de 1 kHz. La longueur du signal est de 1 000 échantillons.

fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*100*t) + randn(size(t));

Obtenez le périodogramme avec fft. Le signal est à valeurs réelles et de longueur paire. Comme le signal est à valeurs réelles, vous n’avez besoin des estimations de puissance que pour les fréquences positives ou négatives. Pour conserver la puissance totale, multipliez toutes les fréquences présentes dans les deux ensembles (fréquences positives et négatives) par un facteur de 2. La fréquence nulle (DC) et la fréquence de Nyquist ne se produisent pas deux fois. Tracez le résultat.

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,pow2db(psdx))
grid on
title("Periodogram Using FFT")
xlabel("Frequency (Hz)")
ylabel("Power/Frequency (dB/Hz)")

Calculez et tracez le périodogramme avec periodogram. Démontrez que les deux résultats sont identiques.

periodogram(x,rectwin(N),N,fs)

mxerr = max(psdx'-periodogram(x,rectwin(N),N,fs))

mxerr = 
3.4694e-18

Entrée avec fréquence normalisée

Utilisez fft pour générer un périodogramme pour une entrée avec une fréquence normalisée. Créez un signal composé d’une onde sinusoïdale et d’un bruit additif N(0,1). L’onde sinusoïdale a une fréquence angulaire de $$\pi /4$$ rad/échantillon.

N = 1000;
n = 0:N-1;
x = cos(pi/4*n) + randn(size(n));

Obtenez le périodogramme avec fft. Le signal est à valeurs réelles et de longueur paire. Comme le signal est à valeurs réelles, vous n’avez besoin des estimations de puissance que pour les fréquences positives ou négatives. Pour conserver la puissance totale, multipliez toutes les fréquences présentes dans les deux ensembles (fréquences positives et négatives) par un facteur de 2. La fréquence nulle (DC) et la fréquence de Nyquist ne se produisent pas deux fois. Tracez le résultat.

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,pow2db(psdx))
grid on
title("Periodogram Using FFT")
xlabel("Normalized Frequency (\times\pi rad/sample)")
ylabel("Power/Frequency (dB/(rad/sample))")

Calculez et tracez le périodogramme avec periodogram. Démontrez que les deux résultats sont identiques.

periodogram(x,rectwin(N),N)

mxerr = max(psdx'-periodogram(x,rectwin(N),N))

mxerr = 
4.4409e-16

Entrée à valeurs complexes avec fréquence normalisée

Utilisez fft pour générer un périodogramme pour une entrée à valeurs complexes avec une fréquence normalisée. Le signal est une exponentielle complexe avec une fréquence angulaire de $$\pi /4$$ rad/échantillon et un bruit N(0,1) à valeurs complexes.

N = 1000;
n = 0:N-1;
x = exp(1j*pi/4*n) + [1 1j]*randn(2,N)/sqrt(2);

Utilisez fft pour obtenir le périodogramme. Comme l’entrée est à valeurs complexes, obtenez le périodogramme sur l’intervalle $$[0,2\pi )$$ rad/échantillon. Tracez le résultat.

xdft = fft(x);
psdx = (1/(2*pi*N)) * abs(xdft).^2;
freq = 0:2*pi/N:2*pi-2*pi/N;

plot(freq/pi,pow2db(psdx))
grid on
title("Periodogram Using FFT")
xlabel("Normalized Frequency (\times\pi rad/sample)")
ylabel("Power/Frequency (dB/(rad/sample))")

Utilisez periodogram pour obtenir et tracer le périodogramme. Comparez les estimations de la PSD.

periodogram(x,rectwin(N),N,"twosided")

mxerr = max(psdx'-periodogram(x,rectwin(N),N,"twosided"))

mxerr = 
2.2204e-16

Voir aussi

Applications

• Signal Analyzer

Fonctions

• fft | periodogram | pspectrum

Ressources pédagogiques

• Analyse de Fourier (Ressources pédagogiques MathWorks)