Matlab　關于功率譜的計算2008-03-11 17:19關于功率譜密度的方法的討論及其實現
直接法： 
直接法又稱周期圖法，它是把隨機序列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));

間接法： 
間接法先由序列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);

改進的直接法： 
對于直接法的功率譜估計，當數據長度N太大時，譜曲線起伏加劇，若N太小，譜的分辨率又不好，因此需要改進。

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]);

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);