%周期圖法功率譜估計(jì)
%exa6-7_period.m, for example 6-7
%to test period.m

clear all;
Fs=1000; %采樣頻率

%產(chǎn)生含有噪聲的序列
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;

%繪制信號(hào)波形
subplot(311)
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)')

%計(jì)算序列的DFT
nfft=1024;
Xk=fft(xn,nfft);
%計(jì)算序列的PSD
Pxx1=abs(Xk).^2/N;
%繪制功率譜圖形
index=0:round(nfft/2-1);
k=index*N/nfft;
plot_Pxx1=10*log10(Pxx1(index+1));
subplot(312)
plot(k,plot_Pxx1);
ylabel('公式直接計(jì)算的功率譜')

%periodogram函數(shù)計(jì)算的功率譜
window=boxcar(length(xn));
[Pxx2,f]=periodogram(xn,window,nfft,Fs);
plot_Pxx2=10*log10(Pxx2(index+1));
subplot(313)
plot(k,plot_Pxx2)
xlabel('periodogram 函數(shù)計(jì)算的功率譜')