mmusic.m

改进的music算法 通信信号处理 功率普估计 波达方向估计

%改进music算法10阵元等距阵波达方向估计（d=λ/2）,2信源，前向-后向 空间平滑 5个子阵（6阵元） 最小范数
clear
w1=pi/4;
w2=pi/3;
kp=2048;
a1=[exp(-j*pi*sin(w1)*[0:9])].';
a2=[exp(-j*pi*sin(w2)*[0:9])].';
A=[a1,a2];
k=1:kp;
s=[cos(k*0.015);cos(k*0.015)];
n=randn(10,kp)+j*randn(10,kp);
x=A*s+n;
r=x*(x)'/kp;
J=rot90(eye(6));
for m=1:5
    b=sprintf('r%d',m);
    c=sprintf('rb%d',m);
    rm=sprintf('rm%d',m);
    F=[zeros(6,m-1) eye(6) zeros(6,5-m)];
    f=F.';
    eval([b,'=','F','*r*','f;']);
    eval([c,'=J*conj(',b,')*J;']);
    eval([rm,'=(',b,'+',c,')/2;']);
end
rx=(rm1+rm2+rm3+rm4+rm5)/5;
[V,D]=eig(rx);
[bn,index]=sort(diag(D));
Vn=V(:,index);
Un=Vn(:,[1:4]);
II=Un*(Un)';
w=[-pi/2:pi/180:pi/2];
for l=1:length(w)
    ao=[exp(-j*pi*sin(w(l))*[0:5])].';
    pm(l)=(ao)'*ao*(inv((ao)'*II*ao));
     p(l)=10*log10(abs(pm(l)));
end
t=w*180/pi;
plot(t,p);grid;
title('改进music算法');
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