exa6_11.m

来自「新版Matlab 7辅助信号处理技术与应用」· M 代码 · 共 28 行

% exa060207_pmtm.m, for example 6.2.7 
% to estimate the PSD of x(n) by the multitaper method (MTM);

clear;
Fs=1000;% 采样频率
% 产生含有噪声的序列
var=sqrt(1/exp(1.0));
n=0:1/Fs:1;
N=length(n);
e=var*randn(1,N);
w0=100*pi;
w1=50*pi;
xn=exp(j*w0*n-j*pi)+exp(j*w1*n-j*0.7*pi)+e;

% 绘制信号波形
subplot(211)
plot(n,abs(xn))
xlabel('n')
title('x(n)=exp(j*w0*n-j*pi)+exp(j*w1*n-j*0.7*pi)+e(n)')

% 用 多窗口法估计序列的功率谱；
nw=3.5;
nfft=512;
p=0.99;
[Pxx,Pxxc,f] = pmtm(xn,nw,nfft,Fs,p);
subplot(212)
psdplot([Pxx Pxxc],f,'hz','db')
title('99% 置信度的PMTM功率谱估计')