%%------------------------------------------------------------------------%%功能：使用MATLAB對平穩隨機信號以改進的功率譜法（welch法）以及Thomson%%      多窗估計法（Thomson Multitaper Method）做功率譜估計%%作者：李鵬%%時間：2009.2.13%%------------------------------------------------------------------------fs=1;n1=0:1/fs:1000;n=n1(:,1:end-1);ln=length(n);vn=randn(1,ln);xn=cos(0.35*pi*n)+cos(0.4*pi*n)+vn;%%************************************************************************%%以下用修正的周期圖法welch法對平穩隨機信號xn作功率譜分析%%Figuer 1 示意隨分段長度（窗函數時域長度）增加，方差特性變差figure(1);[pxx,w]=pwelch(xn,128,64);subplot(2,1,1);plot(w/pi*fs,10*log10(abs(pxx)));%橫軸需根據采樣頻率sample作相應調整%無輸出參數時繪圖的橫軸為rad/sample，對應采樣率為1時繪出正確圖形，否則需調整橫軸axis([0 1 -20 20])xlabel('Normalized Frequency(rad/sample)');ylabel('Power/Frequency(dB/rad/sample)');title('Power Spectral Density Estimate via Welch(modified periodogram method)');grid on;[pxx,w]=pwelch(xn,256,128);subplot(2,1,2);plot(w/pi*fs,10*log10(abs(pxx)));axis([0 1 -20 20])xlabel('Normalized Frequency(rad/sample)');ylabel('Power/Frequency(dB/rad/sample)');grid on;%%************************************************************************%%以下用修正的周期圖法Thomson多窗估計法（MultiTaper Method, MTM）%%對平穩隨機信號xn作功率譜分析%%Figure 2 示意隨窗函數增多，方差特性改善，但頻率分辨率下降[pxx,pxxc,w]=pmtm(xn,2);figure(2);subplot(2,1,1);plot(w/pi*fs,10*log10(abs(pxx)));axis([0 1 -20 20])xlabel('Normalized Frequency(rad/sample)');ylabel('Power/Frequency(dB/rad/sample)');title('Power Spectral Density Estimate via Thomson Multitaper');grid on;[pxx,pxxc,w]=pmtm(xn,4);subplot(2,1,2);plot(w/pi*fs,10*log10(abs(pxx)));axis([0 1 -20 20])xlabel('Normalized Frequency(rad/sample)');ylabel('Power/Frequency(dB/rad/sample)');grid on;