%例程5－9  利用Welch法估計功率譜
% e.g.5-9.m for example5-9;
% to test function pwelch;
 
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(211);
plot(n,xn); xlabel('n');
grid on;
title(' xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn');

% Set parameter
nfft=1024;
window=blackman(256);
noverlap=128;
range='oneside';
% Estimate the PSD 
[Pxx,f]=pwelch(xn,window,noverlap,nfft,N-1,range);
f=(0:length(Pxx)-1)/length(Pxx);
subplot(212);
plot(f,10*log10(Pxx)); 
xlabel('Frequency(Hz)');
ylabel('Power Spectrum(dB)');
title('Use Pwelch Function(Pxx(f)/(N-1)');
grid;