exc_2_2.m

基于AR模型的现代谱估计

function exc_2_2(Ndata, pmodel,Nfft,SNR1, SNR2);
%Ndata is the length of signal x(n) and noise signal w(n)
%pmodel is the phase of AR model
%Nfft is the length of frequency response sample of AR model
%used to generate the PSD
%SNR1 and SNR2 is the signal noise ratio in response to the two signal 
%N = input('Please input the length of data: ');
N = Ndata;
wn = gaussnoise(0,1,N); %generate the guauss distributed noise w(n)
%SNR1 = input('Please input the SNR1: ');
%SNR2 = input('Please input the SNR2: ');
A1 = sqrt(2*(10^(SNR1/10))); %Amplitude of sin signal with SNR1
A2 = sqrt(2*(10^(SNR2/10)));
f1 = 100; %Hz
f2 = 110; %Hz
fs = 1000;
w1 = 2*pi*f1/fs;
w2 = 2*pi*f2/fs;
n = 0:N-1;
xn = A1*sin (w1*n+pi/2) + A2*sin (w2*n+pi/2) + wn';

%p= input('Please input the phase of the AR model: ');
p = pmodel;
%autocorrelation method
[ap_corr, en_corr] = psd_autocorrelation(xn, p);
[H, Wn] = freqz(1, [1,ap_corr'], Nfft);
P_corr = en_corr*(abs(H)).^2;
figure(1)
subplot(2,2,1)
plot(Wn*fs/2/pi,10*log10(P_corr))
title(['自相关法(','N=',num2str(N),',','p=',num2str(p),',','SNR1=',num2str(SNR1),',','SNR2=',num2str(SNR2),')']);
ylabel('信号功率谱Pxx(dB)')
xlabel('频率f(Hz)')

%autocovariance method
[ap_cova, en_cova] = psd_autocovariance(xn, p);
[H, Wn] = freqz(1, [1,ap_cova'], Nfft);
P_cova = en_cova*(abs(H)).^2;
%figure(2)
subplot(2,2,2)
plot(Wn*fs/2/pi,10*log10(P_cova))
title(['协方差法(','N=',num2str(N),',','p=',num2str(p),',','SNR1=',num2str(SNR1),',','SNR2=',num2str(SNR2),')']);
ylabel('信号功率谱Pxx(dB)')
xlabel('频率f(Hz)')

%amendment covariance method
[ap_amcova, en_amcova] = psd_amendmentcovariance(xn, p); 
[H, Wn] = freqz(1, [1,ap_amcova'], Nfft);
P_amcova = en_amcova*(abs(H)).^2;
%figure(3)
subplot(2,2,3)
plot(Wn*fs/2/pi,10*log10(P_amcova))
title(['修正协方差法(','N=',num2str(N),',','p=',num2str(p),',','SNR1=',num2str(SNR1),',','SNR2=',num2str(SNR2),')']);
ylabel('信号功率谱Pxx(dB)')
xlabel('频率f(Hz)')

%Burg method
[ap_burg, en_burg] = psd_burg(xn, p);
[H, Wn] = freqz(1, [1,ap_burg'], Nfft);
P_burg = en_burg*(abs(H)).^2;
%figure(4)
subplot(2,2,4)
plot(Wn*fs/2/pi,10*log10(P_burg))
title(['伯格法(','N=',num2str(N),',','p=',num2str(p),',','SNR1=',num2str(SNR1),',','SNR2=',num2str(SNR2),')']);
ylabel('信号功率谱Pxx(dB)')
xlabel('频率f(Hz)')