%图 均方误差随快拍数的渐变趋势
%成功概率随快拍数的渐变趋势
clear all
close all
clc
tic;
%参数设定
lamda=1;
kelm=8;
dd=0.5*lamda;
d=0:dd:(kelm-1)*dd;
theta=[-10 15 30];
doanum=length(theta);
SNR=10;
montenum=50;
SNAPM=50:50:1000;
snapnum=length(SNAPM);
Narg=2;%曲线的种类
data=zeros(snapnum,montenum,Narg,doanum);
for snapindex=1:snapnum
snap=SNAPM(snapindex);
snap
for monindex=1:montenum
monindex
S=randn(doanum,snap);
A=exp(-j*2*pi*d.'*sin(theta*pi/180)/lamda);
X=A*S;
X=awgn(X,SNR,'measured');
Rxx=X*X'/snap;
[V,D]=eig(Rxx);
Ud=diag(D)';
[Ud,I]=sort(Ud);
Ud=fliplr(Ud);
EV=fliplr(V(:,I));
% MUSIC
rxv=inv(Rxx);
En=EV(:,doanum+1:kelm);
theta0=[-90:0.05:90];
doa0num=length(theta0);
for n=1:doa0num
a=exp(-j*2*pi*d.'*sin(theta0(n)*pi/180)/lamda);
Pmusic(n)=1/(a'*En*En'*a);
end
Mpeaks=findpeak(Pmusic,theta0,doanum);
%LS-ESPRIT
Ex=EV(1:kelm-1,1:doanum);
Ey=EV(2:kelm,1:doanum);
[pes,qes]=eig(inv(Ex'*Ex)*Ex'*Ey);
ESDOA=zeros(1,doanum);
for i=1:doanum;
ESDOA(i)=real(asin(-j*(log(qes(i,i)))*lamda/(-2*pi*dd))*180/pi);
end;
% 结果整理
data(snapindex,monindex,1,:)=sort(Mpeaks);
data(snapindex,monindex,2,:)=sort(ESDOA);
end
end
%MSE
mse=zeros(snapnum,doanum,Narg);
for snapindex=1:snapnum
for monindex=1:montenum
for doaindex=1:doanum
%MUSIC
mse(snapindex,doaindex,1)=mse(snapindex,doaindex,1)+(data(snapindex,monindex,1,doaindex)-theta(doaindex))^2;
%Esprit
mse(snapindex,doaindex,2)=mse(snapindex,doaindex,2)+(data(snapindex,monindex,2,doaindex)-theta(doaindex))^2;
end
end
end
%MMSE
mmse=sqrt(mse/montenum);
mmseM=(mmse(:,1,1)+mmse(:,2,1)+mmse(:,3,1))/3;
mmseE=(mmse(:,1,2)+mmse(:,2,2)+mmse(:,3,2))/3;
figure(2)
semilogy(SNAPM,mmseM,'-b^')
hold on
semilogy(SNAPM,mmseE,'-kp')
grid on;
xlabel('快拍数')
ylabel('均方误差（°）')
legend('MUSIC算法','ESPRIT算法');
toc;
评论0