經典功率譜估計

1、直接法：
直接法又稱周期圖法，它是把隨機序列x(n)的N個觀測數據視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，得X(k)，然后再取其幅值的平方，并除以N，作為序列x(n)真實功率譜的估計。
Matlab代碼示例：
clear;
Fs=1000; %采樣頻率
n=0:1/Fs:1;
%產生含有噪聲的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));

2、間接法：
間接法先由序列x(n)估計出自相關函數R(n)，然后對R(n)進行傅立葉變換，便得到x(n)的功率譜估計。
Matlab代碼示例：
clear;
Fs=1000; %采樣頻率
n=0:1/Fs:1;
%產生含有噪聲的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
cxn=xcorr(xn,'unbiased'); %計算序列的自相關函數
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);

3、改進的直接法：
對于直接法的功率譜估計，當數據長度N太大時，譜曲線起伏加劇，若N太小，譜的分辨率又不好，因此需要改進。
3.1、Bartlett法
Bartlett平均周期圖的方法是將N點的有限長序列x(n)分段求周期圖再平均。
Matlab代碼示例：
clear；
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(length(n)); %矩形窗
noverlap=0; %數據無重疊
p=0.9; %置信概率
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);
pause;
figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

3.2、Welch法
Welch法對Bartlett法進行了兩方面的修正，一是選擇適當的窗函數w(n)，并再周期圖計算前直接加進去，加窗的優點是無論什么樣的窗函數均可使譜估計非負。二是在分段時，可使各段之間有重疊，這樣會使方差減小。
Matlab代碼示例：
clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %數據無重疊
range='half'; %頻率間隔為[0 Fs/2]，只計算一半的頻率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);

figure(1)
plot(f,plot_Pxx);

pause;

figure(2)
plot(f,plot_Pxx1);

pause;

figure(3)
plot(f,plot_Pxx2);