隨機信號功率譜密度估計
--By xzd1575
一、實驗目的
1．深入理解隨機信號功率譜密度估計
2．掌握在Matlab平台上進行信號功率譜密度估計的基本方法
二、實驗原理
1. 隨機信號功率譜密度定義
定義隨機信號信號的功率譜
為
其中
為隨機信號的自相關函數。
功率譜反映了信號的功率在頻域隨頻率
分布，因此
又稱為功率譜密度。[1]
2. 經典譜估計（非參數譜估計）方法簡介
經典譜估計的方法主要包括兩種方法：周期圖法和自相關法。
周期圖法[1]（直接法）
周期圖法又稱為直接法，它是把隨機信號
的N點觀察數據
視為一個能量有限信號，直接取
的傅里葉變換，得
，然后再取其幅值的平方，並除於N，作為對
真實功率譜
的估計。以
表示用周期圖法估計的功率譜，則
自相關法[1]（間接法）
此方法的理論基礎是維納-辛欽定理。1958年Blackman和Tukey 給出了這一種方法的具體實現，即由
估計出自相關函數
，然后對
的功率譜，記之為
，並以此作為對
的估計，即
因為這種方法求出的功率譜是通過自相關函數間接得到的，所以稱為間接法，又稱自相關法或BT法。當M較小時，上式計算量不是很大，因此，該方法是在FFT問世之前（即周期圖法被廣泛應用之前）常用的譜估計方法。
3. 參數模型譜估計方法簡介[1]
參數模型法是現代譜估計的主要內容，參數模型法的思路如下。
三、實驗步驟
1. 構造模擬信號
構造模擬信號
2. 使用經典譜估計方法對信號進行譜估計
1) 周期圖法
Matlab中Prierdogram()函數就是運用周期圖法進行譜估計。調用格式如下：
[psdestx,Fxx] = periodogram(xn,rectwin(length(xn)),length(xn),Fs);
其中輸入參數xn為待估計的離散信號，rectwin(length(xn))表示窗長為xn點的矩形窗（rectangle window）,Fs表示采樣頻率。
輸出參數Fxx表示頻率，psdestx為對應Fxx頻率的功率譜密度。
為了使周期圖法得到的功率譜密度更為平滑，提出了許多改進的方法，Welch平均周期圖法就是其中一種，在matlab中pwelch()函數就是使用該方法進行功率譜估計，pwelch()函數的調用格式如下：
pwelch(xn,hamming(256),128,1024,Fs)
輸入參數xn為輸入信號，hamming(256)為窗長為256的漢明窗，Fs為信號采樣頻率。調用后可繪制得到信號功率譜密度圖，如需要觀察得到的功率譜密度數值，可以添加相應的輸出參數，相應可以參閱matlab幫助文檔。
2) 相關函數法
相關函數法是先求信號是自相關函數，再根據維納辛欣定理，功率譜密度就是自相關函數的傅里葉變換，對自相關函數求傅里葉變換，得到功率譜密度。
需要用到matlab中xcorr()函數，其調用格式如下：
cx=xcorr(xn,'biased');
其中輸入參數xn為待求自相關函數的信號，'biased'表示使用有偏差的自相關函數求法。
輸出參數cx即為信號xn的自相關函數。
3. 使用現代譜估計方法對信號進行譜估計
伯格（Brug）譜估計是一種AR譜估計方法，可調用matalb中pburg函數，其調用格式如下：
pburg(xn,5,1024,Fs)
輸入參數xn為信號，Fs為采樣頻率。調用后可繪制得到信號功率譜密度圖，如需要觀察得到的功率譜密度數值，可以添加相應的輸出參數，相應可以參閱matlab幫助文檔。
四、實驗結果與分析
1. 經典譜估計方法和現代譜估計方法比較
圖4.1 不同功率譜估計方法比較
如圖4.1所示，對比周期圖法（periodogram）和平均周期圖法（Welch），驗證了Welch法得到的圖要比周期圖法得到的功率譜密度圖光滑。自相關法和周期圖法得到的功率譜估計在140Hz和150Hz處鋒比較尖銳，頻率分辨率要比Welch平均周期圖法高。現代AR譜估計Brug方法同樣可以在140Hz和150Hz處得到尖銳的譜峰，同時其估計的功率譜密度圖也很平滑。
2. AR譜估計中模型階數對譜估計結果的影響
(a) (b) (c)
圖4.2 AR模型階數對譜估計的影響
(a)5階 (b)14階 (c)20階
如圖4.2，對比不同AR模型階數對功率譜估計的影響，發現階數較低時，在140Hz-150Hz頻率范圍左右，只出現一個譜峰，沒有得到實際的兩個譜峰，頻率分辨率不夠，隨着模型階數的增加，到階數達到14時，可以有效地區分140Hz和150Hz處的兩個譜峰，有較好的頻率分辨率，隨着模型階數的繼續增加，在真峰（140Hz和150Hz）附近的假峰會隨着增多。
五、實驗結論
通過對比經典和現代不同譜估計方法，可以發現，現代譜估計方法既有較好的頻率分辨率，又是能使功率譜密度較為平滑，可以很到的得到信號譜峰。
現代AR譜估計中，模型的階數選擇是一個很重要的問題，選擇合適的階數，可以有效的檢查出有效信號的譜峰，如果模型階數過低，則頻率分辨率不夠，可能會丟失有效信號譜峰，如果模型階數過高，則可能出現假峰。
六、參考文獻
[1] 胡廣書． 數字信號處理：理論、算法與實現（第三 版）[M]. 北京：清華大學出版社，2012．
七、附錄
Matlab 實驗程序
%% 構造模擬信號
close all
Fs = 1000; % Sampling frequency
t = (0:Fs)/Fs; % One second worth of samples
A = [1 2 3]; % Sinusoid amplitudes
f = [150;140]; % Sinusoid frequencies
x = A*sin(2*pi*f*t);%1*sin(2*pi*150*t)+2*sin(2*pi*140*t)
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
figure(1)
plot(t(1:200),xn(1:200));
xlabel('Time(s)')
ylabel('x(t)')
title('Signal x(t) with Additional Gussian White Noise')
set(gcf,'color',[1 1 1]);
%% 周期圖法
figure(2);
subplot(2,2,1)
[psdestx,Fxx] = periodogram(xn,rectwin(length(xn)),length(xn),Fs);
plot(Fxx,10*log10(psdestx)); grid on;
xlabel('Frequency(Hz)'); ylabel('Power/frequency (dB/Hz)');
title('Periodogram Power Spectral Density Estimate');
axis([0 500 -60 10])
%% 自相關函數法（BT法）
subplot(2,2,2)
cx=xcorr(xn,'biased'); %計算自相關函數
cxdft = fft(cx);
psdx = abs(cxdft)/Fs;
freq = 0:Fs/length(psdx):Fs/2;
plot(freq,10*log10(psdx(1:length(freq)))); grid on;
title('AutoCorrelation Power Spectral Density Estimate');
xlabel('Frequency (Hz)'); ylabel('Power/frequency (dB/Hz)');
axis([0 500 -60 10])
%% Welch 平均周期法
subplot(2,2,3)
pwelch(xn,hamming(256),128,1024,Fs);
axis([0 500 -60 10])
%% Burg法 AR參數譜估計
figure(2)
subplot(2,2,4)
pburg(xn,14,1024,Fs)
axis([0 500 -60 10])
set(gcf,'color',[1 1 1]);
%% 討論不同的AR階數對Brug法的影響
figure(3)
subplot(1,3,1)
pburg(xn,5,1024,Fs)
axis([0 500 -60 10])
title('Order 5')
subplot(1,3,2)
pburg(xn,14,1024,Fs)
axis([0 500 -60 10])
title('Order 14')
subplot(1,3,3)
pburg(xn,20,1024,Fs)
axis([0 500 -60 10])
title('Order 20')
set(gcf,'color',[1 1 1]);