%例程5－8  利用Bartlett法估計功率譜
% e.g.5-8.m for example5-8;
% to test function psd;
 
clear all;
% Generate the signal with noise and display
N=1024;
n=0:1/(N-1):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(311);
plot(n,xn); xlabel('n');
title(' xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn');
% Estimate the PSD
window=boxcar(N/8);
[Pxx,Pxxc,f]=psd(xn,1024,1000,window,0,0.9);
 f=(0:length(Pxx)-1)/length(Pxx);
plot_Pxx=10*log10(abs(Pxx));
subplot(312);
plot(f,plot_Pxx);
grid on
xlabel('Frequency'), 
ylabel('Power Spectrum (dB)');
%Estimate the PSD with 90% confidence interval
subplot(313);
plot(f,10*log10(Pxx-Pxxc(:,1)), f,10*log10(Pxx+Pxxc(:,1)));
xlabel('Frequency(Hz)');
ylabel('Power Spectrum(dB)');
title('PSD with 90% confidence interval');
grid;