pjh2.m

PSD using BT/AR/ method,get the auto-correlation funtion

%get the simulate data
close all;
clear all;
N1=10000;	%总点数
N=1000;		%采样点数
x1=zeros(1,N1);
x=zeros(1,N);
var_v=0.1;
V=sqrt(var_v)*randn(1,N1);
w0=1.2*pi;
for n=1:N1
x1(n)=exp(i*pi*n-i*pi)+exp(i*w0*n-i*0.7*pi)+V(n);
end
x=x1(N1-N+1:N1); %选取后面1000个点
clear x1 V;
%-------------------------------------------------------------------------------------------------------------------
%get the auto-correlation funtion
for m=1:N
Pre_r(m)=0;
for n=m:N
Pre_r(m)=x(n)*conj(x(n-m+1))+Pre_r(m);
end
Pre_r(m)=Pre_r(m)/N;
end
N2=500;
r=zeros(1,N2);
r0=Pre_r(1);    %get the r(0)
r=Pre_r(2:501); %get the r(m),m=1,2,...,500
%-------------------------------------------------------------------------------------------------------------------
%PSD using BT method with estimated data
w=linspace(0,2*pi,1000);
Pbtsum=0;
for m=1:500
Pbtsum=r(m)*cos(w*m)+Pbtsum;
end
P_bt=Pbtsum*2+r0;
subplot(1,3,1),plot(w,abs(P_bt));
grid on;
title('BT method');
%-------------------------------------------------------------------------------------------------------------------
%Let P=3, construct a 4th order correlation matrix R
[X,R]=corrmtx(x,3);
[V,D]=eig(R);
En=V(1:4,1:2);
%PSD estimation using MUSIC and MVDR method
for l=1:1000
    a=[1,exp(-j*w(l)),exp(-j*w(l)*2),exp(-j*w(l)*3)];
    P_music(l)=1/(abs(conj(a)*V(1:4,1))^2+abs(conj(a)*V(1:4,2))^2);
    P_mvdr(l)=1/(conj(a)*R^(-1)*a');
end
subplot(1,3,2),plot(w,P_music);
grid on;
title('MUSIC method');
subplot(1,3,3),plot(w,abs(P_mvdr));
grid on;
title('MVDR method');
end