%例程5－7  功率譜估計的周期圖法

% e.g.5-7.m for example5-7;
% to test function of periodogram;
 
%Estimate the PSD with periodogram method(By formula directly)
clear;
% Case 1: N=256
N=256;
n=0: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;
Pxx=10*log10(abs(fft(xn).^2)/N);
f=(0:length(Pxx)-1)/length(Pxx);
subplot(311);
plot(f,Pxx);
xlabel('Frequency(Hz)');
ylabel('Power Spectrum(dB)');
title('Periodogram  N=256');
grid;
%===========================================
%Case 2:N=1024
N2=1024;
n=0:N2-1;
f1=0.1;
f2=0.9;
wn2=randn(1,N2);
xn2=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn2;
Pxx2=10*log10(abs(fft(xn2).^2)/N2);
f=(0:length(Pxx2)-1)/length(Pxx2);
subplot(312);
plot(f,Pxx2);
xlabel('Frequency(Hz)');
ylabel('Power Spectrum(dB)');
title('Periodogram  N=1024');
grid;
%Estimate the PSD with function 'periodogram'(For Case 2 only)
window=boxcar(length(xn2));
[Pxx3,f]=periodogram(xn2,window,N2,[]);
plot_Pxx3=10*log10(Pxx3);
f=(0:length(Pxx3)-1)/length(Pxx3);
subplot(313);
plot(f,plot_Pxx3);
xlabel('Frequency(Hz)');
ylabel('Power Spectrum(dB)');
title('Use Periodogram Function, N=1024');
grid;