1. 程式人生 > 實用技巧 >DSP - 使用 FFT 獲得功率頻譜密度估計

DSP - 使用 FFT 獲得功率頻譜密度估計

具有指定取樣率的偶數長度輸入

對於一個取樣率為 1 kHz 的偶數長度訊號,分別使用fftperiodogram獲得其週期圖。比較二者的結果。

建立一個含N(0,1) 加性噪聲的 100 Hz 正弦波訊號。取樣頻率為 1 kHz。訊號長度為 1000 個取樣點。使用隨機數生成器的預設設定以獲得可重現的結果。

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

使用fft獲取週期圖。訊號是偶數長度的實數值訊號。由於訊號是實數值訊號,您只需要對正負頻率之一進行功率估計。
為了保持總功率不變,將同時在兩組(正頻率和負頻率)中出現的所有頻率乘以因子 2。零頻率 (DC) 和 Nyquist 頻率不會出現兩次。繪製結果。
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;

計算並使用periodogram繪製週期圖。二者的結果相同。

periodogram(x,rectwin(length(x)),length(x),Fs)

具有歸一化頻率的輸入

通過fft為使用歸一化頻率的輸入生成周期圖。建立一個帶N(0,1) 加性噪聲的正弦波訊號。該正弦波的角頻率為π/4弧度/取樣點。使用隨機數生成器的預設設定以獲得可重現的結果。

rng default
n = 0:999;
x = cos(pi/4*n) + randn(size(n));

使用fft獲取週期圖。訊號是偶數長度的實數值訊號。由於訊號是實數值訊號,您只需要對正負頻率之一進行功率估計。為了保持總功率不變,將同時在兩組(正頻率和負頻率)中出現的所有頻率乘以因子 2。零頻率 (DC) 和 Nyquist 頻率不會出現兩次。繪製結果。

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;

計算並使用periodogram繪製週期圖。二者的結果相同。

periodogram(x,rectwin(length(x)),length(x))

具有歸一化頻率的複數值輸入

使用fft為具有歸一化頻率的複數值輸入生成周期圖。採用一個帶N(0,1) 復噪聲的復指數訊號,角頻率為π/4弧度/取樣點。採用隨機數生成器的預設設定,以獲得可重現的結果。

rng default
n = 0:999;
x = exp(1j*pi/4*n) + [1 1j]*randn(2,length(n))/sqrt(2);

使用fft獲得週期圖。由於輸入是複數值,此處求[0,2π)弧度/取樣點區間內的週期圖。繪製結果。

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

使用periodogram獲取並繪製週期圖。比較 PSD 估計。

periodogram(x,rectwin(length(x)),length(x),'twosided')

Reference

  1. MathWorks