%%使用周期圖法以及改進的周期法實現譜估計
%%----------------------------------------
%%使用周期圖法實現譜估計
close all;
clear all;
clc;

N=128;
n=1:128;
a1=sqrt(20);
a2=sqrt(2);
w=randn(128,1);
x=a1*sin(2*pi*0.2*n)+a2*sin(2*pi*0.213*n)+w';%數據源序列

%%-----------------方法一-----------------------
figure(2);
subplot(211);
yfft=fft(x,N);
plot(n/N,10*log10(abs(yfft)));
xlabel('f');
ylabel('dB');
title('周期圖法實現譜估計');
grid on;
axis([0,1,-20,70]);
%-------------------方法二-------------------
%%改進的周期法實現譜估計
subplot(212);
w=hanning(128)';
plot(n/N,10*log10(abs(fft(w.*x))),'r');
xlabel('f');
ylabel('dB');
legend('改進的周期圖法實現譜估計');
title('加hanning窗后的周期圖法');
grid on;
%------------------方法三--------------------
%利用自相關函數間接估計功率普
%xcorr
xxcorr=xcorr(x,'unbiased');
figure(3);
n1=1:2*N-1;
plot(n1/(2*N),10*log10(abs(fft(xxcorr))));
xlabel('f');
ylabel('dB');
title('自相關函數間接譜估計法');
grid on;