1-d~music.m

一维线形阵元的music算法DOA估计

%构造信号源
n=-pi:0.01:pi;                           %范围    
x1=100*cos(10*sqrt(2.0)*pi*n+pi/3)';      %构造三个余弦信号
x2=100*cos(20*sqrt(2.0)*pi*n+pi/3)';
x3=100*cos(50*sqrt(2.0)*pi*n+pi/3)';
save shy3 x1 x2 x3                        %保存信号源

%天线阵元及接收信号
clear
close all
i=sqrt(-1);    %复数单位
j=i;
m=8;      %天线阵元数
p=3;      %信号个数
angle1=45;angle2=60;angle3=75;    %信号来向
angle=[angle1;angle2;angle3];     
nn=600;   %采样数 
SN1=20;SN2=10;SN3=12;  %信噪比
SN=[SN1;SN2;SN3];
degrad=pi/180;  %将角度制化为弧度制    

%信号和噪声
load shy3    %载入信号源
TT=1:nn;    %采样数矩阵
S=[x1(TT),x2(TT),x3(TT)].'; 
nr=randn(m,nn);  %建立高斯随机序列
ni=randn(m,nn);
E=nr+ni*j;        %建立高斯噪声
Ps=S*S'/nn;    %方向向量
Ps1=diag(Ps);
refp=2*10.^(SN/10);
tmp=sqrt(refp./Ps1);
S1=diag(tmp)*S;

%计算协方差矩阵，特征值分解
A=zeros(m,p);
k=[0:m-1]';
for t=1:p
    A(:,t)=exp(-j*pi*k*sin(angle(t)*degrad));
end
X=A*S1+E;           %第k个阵元上的接收信号
Rxx=X*X'/nn;       %接收信号的协方差
[U,s,v]=svd(Rxx);     %计算特征值
Vs=U(:,1:p);          
Vu=U(:,p+1:m);

%计算空间谱函数
for ang=0:90
    k=[0:m-1]';
    AA=exp(-j*pi*k*sin(ang*degrad));
    WW=AA'*Vu*Vu'*AA;
    Pmusic(ang+1)=abs(8./WW);
end

%做出空间谱图
aangle=[0:90];
plot(aangle,Pmusic);
title('MUSIC线形谱');    %图表标题，x轴y轴标签及范围
xlabel('角度');
ylabel('角谱');
grid on