mn.m

最小范数算法 波达方向估计 通信信号处理

%MN最小范数算法10阵元等距阵波达方向估计（d=λ/2）
clear
w1=pi/4;
w2=-pi/3;
w3=pi/4.5;
w4=-pi/3.5;
kp=2048;
a1=[exp(-j*pi*sin(w1)*[0:9])].';
a2=[exp(-j*pi*sin(w2)*[0:9])].';
a3=[exp(-j*pi*sin(w3)*[0:9])].';
a4=[exp(-j*pi*sin(w4)*[0:9])].';
A=[a1,a2,a3,a4];
k=1:kp;
s=[1.3*cos(k*0.015);sin(k*0.05);cos(k*0.02);sin(k*0.035)];
n=randn(10,kp)+j*randn(10,kp);
x=A*s+n;
r=x*(x)'/kp;
[V,D]=eig(r);
[b,index]=sort(diag(D));
Vn=V(:,index);
Un=Vn(:,[1:6]);
II=Un*(Un)';
e1=eye(10,1);
w=[-pi/2:pi/180:pi/2];
for l=1:length(w)
    ao=[exp(-j*pi*sin(w(l))*[0:9])].';
    pm(l)=abs((ao)'*II*e1)^2;
    p(l)=10*log10(1/pm(l));
end
t=w*180/pi;
plot(t,p);grid
title('最小范数算法')
ylabel('P(\theta) (dB)');
xlabel('\theta (deg)');
[Wn Pn]=ginput;
hold on
plot(Wn,Pn,'or');
for k=1:length(Wn)
    text(Wn(k),Pn(k),['  \theta=',num2str(Wn(k)),'P=',num2str(Pn(k))]);
end
hold off