rmse.m

谱估计源程序：是谱估计中的经典算法esprit,MUSIC等算法的源程序

%%%PRO-ESPRIT
clear all; close all;clc;
J=sqrt(-1);
source_number=3;
source_doa=[40 70 80];
sensor_number=10;
snapshort_number=100;
%snr=10;
 
mot=200;
k=0;
ef1=zeros(1,25);
ef2=zeros(1,25);
ef3=zeros(1,25);
for snr=1:1:25
    
k=k+1;
A=exp(-J*(0:sensor_number-1)'*pi*sin(source_doa*pi/180));
%theta=[20,40,60];
%fei=[10,30,50];
%for i=1:3
%A(:,i)=exp(-j*pi*cos(theta(i)*pi/180-2*pi*(0:9)/10)*sin(fei(i)*pi/180));
%end
for j=1:1:mot 
s=(randn(source_number,snapshort_number)+J*randn(source_number,snapshort_number));
x=A*s;
X=awgn(x,snr);


M=[exp(j*pi*sin(source_doa(1)*pi/180)),exp(j*pi*sin(source_doa(2)*pi/180)),exp(j*pi*sin(source_doa(3)*pi/180))];
M=diag(M);
y=A*M*s ;
Y=awgn(y,snr);
 
[Ux,Sx,Vx]=svd(X);
Uxs=Ux(:,1:source_number);
Sxs=Sx(1:source_number,1:source_number);
Vxs=Vx(:,1:source_number);
fy=inv(Sxs)*Uxs'*Y*Vxs;
% compute thetae
[Ufy,Sfy]=eig(fy);
for i=1:3
    doa(1,i)=asin(angle(Sfy(i,i))/pi)*180/pi;
end

f1=(doa(1,1)-source_doa(1,1))^2;
f2=(doa(1,2)-source_doa(1,2))^2;
f3=(doa(1,3)-source_doa(1,3))^2;
ef1(1,k)=f1+ef1(1,k);
ef2(1,k)=f2+ef2(1,k);
ef3(1,k)=f3+ef3(1,k);
end

ef1(1,k)=sqrt(ef1(1,k)/mot);
ef2(1,k)=sqrt(ef1(1,k)/mot);
ef3(1,k)=sqrt(ef1(1,k)/mot);
end

plot(snr,ef1,'-diamond',snr,ef2,'-square',snr,ef3,'-*');
title('Azimuth angle Estimation Performance');
xlabel('snr');
ylabel('RMSE(deg)');
legend('reflection 1','reflection 2','reflection 3');
 
axis 'auto y';