clear all
clc
p=3;
M=7;
fc=1e9;
DOA=[20,40,60]/180*pi;
DOA1=[20,40,60];
disp('输入信号来向');
fs=3*fc;
N=512;
T=1/fs;
c=3e8;
d=0.15;
t=0:T:(N-1)*T;
s1=sqrt(2)*cos(2*pi*fc*t);
mt=sqrt(2)*cos(2*pi*5e8*t);
s2=mt.*cos(2*pi*fc*t);
mt=sqrt(2)*cos(2*pi*1e7*t+pi/8);
s3=(mt).*cos(2*pi*fc*t);
cor=corrcoef([s1;s2;s3]')
ss=[s1;s2;s3];
s=ss(1:p,:);
A=zeros(M,p);
for k=1:p
for kk=1:M
A(kk,k)=exp(-j*2*pi*fc*(kk-1)*d*sin(DOA(k))/c);
end
end
y1=A*s;
RMSE=zeros(1,31);ex=50;
for snr= -10:20;
y=awgn(y1,snr);
R=y*y'/N;%赫米特矩阵
pg=p;
[v,dd]=eig(R);
if(dd(1,1)>dd(2,2))
Un=v(:,pg+1:M);
Us=v(:,1:pg);
else
Un=v(:,1:(M-pg));
Us=v(:,(M-pg+1):M);
end
fai=pinv(Us(1:(M-1),:))*Us(2:M,:);
[v1,dd1]=eig(fai);
pha=zeros(1,3);
theta=zeros(1,3);
for i=1:p
pha(1,i)=angle(dd1(i,i));
theta(1,i)=asin(-pha(1,i)*c/2/pi/d/fc)*180/pi;
end
for k=1:(p-1)
for i=1:(p-k)
if(theta(i)>theta(i+1))
temp=theta(i);
theta(i)=theta(i+1);
theta(i+1)=temp;
end
end
end
RMSE(snr+11)=sqrt((theta-DOA1)*((theta-DOA1).'))/ex;
end
snr=-10:20;
plot(snr,RMSE,'r*-');axis([-10 20 0 0.25]);
xlabel('信噪比/dB');ylabel('均方根误差/degree');
legend('M=7,p=3,N=512,snr=-10:20,ex=50');grid on;
