clear;close all;clc;
L=1024;%采樣點數為1024
SNR1=-10;SNR2=-10;%SNR1和SNR2分別為兩路正弦信號的信噪比，單位為dB
A1=sqrt(2*10^(SNR1/10));%根據dB單位的信噪比，換算出兩路信號的振幅
A2=sqrt(2*10^(SNR2/10));
f1=100;f2=110;fs=1000;%兩路信號的頻率
Num_fft=1024;%進行FFT變換時的點數
N=40;%譜估計的階數
randn('state',0);
w=randn(1,L);%產生噪聲序列
t=0:1/fs:(L-1)/fs;%根據采樣率劃分采樣時間
x=A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+w;%產生觀測信號
%xlim([0,500]);

f_begin=0;
f_end=200;
figure(1)
[Pxx1,f]=pyulear(x,N,Num_fft,fs);%用Levinson法進行功率譜估計
Pxx1=Pxx1/min(Pxx1);%找出譜密度中的最小值，作歸一化處理
Pxx1=10*log10(Pxx1);%產生譜密度的dB數據
subplot(2,2,1);
plot(f,Pxx1);grid on;xlim([f_begin,f_end]);
title('Levinson法估的功率譜密度曲線');
xlabel('頻率f');
ylabel('功率譜密度 (dB)');

[Pxx2,f]=pburg(x,N,Num_fft,fs);  %用Burg法進行功率譜估計
Pxx2=Pxx2/min(Pxx2);%找出譜密度中的最小值，作歸一化處理
Pxx2=10*log10(Pxx2);%產生譜密度的dB數據
subplot(2,2,2);
plot(f,Pxx2);grid on;xlim([f_begin,f_end]);
title('Bug法估的功率譜密度曲線');
xlabel('頻率f');
ylabel('功率譜密度 (dB)');

%figure(2)
[Pxx3,f]=pcov(x,N,Num_fft,fs);%用協方差法進行功率譜估計
Pxx3=Pxx3/min(Pxx3);%找出譜密度中的最小值，作歸一化處理
Pxx3=10*log10(Pxx3);%產生譜密度的dB數據
subplot(2,2,3);
plot(f,Pxx3);grid on;xlim([f_begin,f_end]);
title('協方差法估的功率譜密度曲線');
xlabel('頻率f');
ylabel('功率譜密度 (dB)');

[Pxx4,f]=pmcov(x,N,Num_fft,fs);  %用修正協方差法進行功率譜估計
Pxx4=Pxx4/min(Pxx4);%找出譜密度中的最小值，作歸一化處理
Pxx4=10*log10(Pxx4);%產生譜密度的dB數據
subplot(2,2,4);
plot(f,Pxx4);grid on;xlim([f_begin,f_end]);
title('修正協方差法估的功率譜密度曲線');
xlabel('頻率f');
ylabel('功率譜密度 (dB)');