welch_mtm.m

使用MATLAB编程

%%------------------------------------------------------------------------
%%功能：使用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;