%例程5－10  利用MTM法估計功率譜
% e.g.5-10.m for example5-10.
% to test function pmtm;
 
clear all;
% Generate the signal with noise and display
N=1024;
n=0:1/(N-1):1;     %Sampling frequency N-1
f1=0.1;
f2=0.9;
wn=randn(1,N);
xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn;
subplot(211);
plot(n,xn); xlabel('n');
grid on;
title(' xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn');
% Set parameter
p=0.98;
nw=3.5;
nfft=512;
% Estimate the PSD using multitaper method(MTM)
[Pxx,Pxxc,f]=pmtm(xn,nw,nfft,N-1,p);
f=(0:length(Pxx)-1)/length(Pxx);
subplot(212);
plot(f,10*log10([Pxx Pxxc])); 
xlabel('Frequency(Hz)');
ylabel('Power Spectrum(dB)');
title('PSD with 98% confidence interval');
grid on;