mmusic.m

DOA经典算法的MATLAB程序

clear
 clc
dd=1; %天线阵元的间距
a=2; %信号中心波长
b=10; %信号幅度
SNR=-10; %噪声的方差
q11=90*(pi/180);
q12=90*(pi/180);
q21=70*(pi/180);
q22=100*(pi/180);
q31=100*(pi/180);
q32=60*(pi/180);
q41=70*(pi/180);
q42=90*(pi/180);  %四个输入信号的方向（下标中第一个数字表示第几个信号，第二个数字表示是哪一个角）
A1=zeros(9,1);%求阵因子
A2=zeros(9,1);
A3=zeros(9,1);
A4=zeros(9,1);
for m=0:2
    for i=0:2
        w=m*3+i+1;
        A1(w,1)=[exp(j*(2*pi*dd/a)*(m*cos(q11)*cos(q12)+i*sin(q11)*cos(q12)))]';
        A2(w,1)=[exp(j*(2*pi*dd/a)*(m*cos(q21)*cos(q22)+i*sin(q21)*cos(q22)))]';
        A3(w,1)=[exp(j*(2*pi*dd/a)*(m*cos(q31)*cos(q32)+i*sin(q31)*cos(q32)))]';
        A4(w,1)=[exp(j*(2*pi*dd/a)*(m*cos(q41)*cos(q42)+i*sin(q41)*cos(q42)))]';
    end
end
A=[A1,A2,A3,A4];%得出A 矩阵
n=1:2000;
v1=.05;%四信号的频率
v2=.05;
v3=.05;
v4=.05;
d=[b*cos(v1*n);b*sin(v2*n);b*sin(v3*n);b*sin(v4*n)];%输入信号矢量
randn('state',1);
en=sqrt(10^SNR/10)*randn(9,2000);%高斯白噪声
U=A*d+en;%总的输入信号
%MUSIC算法部分
U1=(U)';
c=cov(U*U1);%总输入信号的协方差矩阵
[s,h]=eig(c);%求协方差矩阵的特征矢量及特征值
Vn=s(:,[5:9]);%取出与最小特征值对应的特征矢量
ch=Vn(1,:);Vnn=Vn([2:9],:);
dn=Vnn*ch'*inv(ch*ch');
DN=zeros(9,1);
DN(1,1)=1;DN([2:9],:)=dn
q1b=[pi/180:pi/180:pi];
num=0;
for n1=1:length(q1b)
    num=num+1
    q1a(n1)=q1b(n1);
    for n2=1:length(q1b)
        q1aa(n2)=q1b(n2);
        for qa=0:2
            for qb=0:2
                x=qa*3+qb+1;
                A1a(x,1)=[exp(j*(2*pi*dd/a)*(qa*cos(n1)*cos(n2)+qb*sin(n1)*cos(n2)))]'; 
            end
        end 
        Pmusic(n1,n2)=inv((A1a)'*DN*(DN)'*A1a);%应用MUSIC 法估计输出
        P1=abs(Pmusic);    
        Tb(n1,n2)=q1aa(n2);
        T(n1,n2)=q1a(n1);
    end
end
%绘出应用MUSIC 法估计的波达方向图
clf;
t1=T*180/pi;
t2=Tb*180/pi;
surf(t1,t2,P1);grid
xlabel('MN algorithm')