music.m

谱估计的非参数DOA(Dorection of Arrive)方法----MUSIC(Multi Signal Classification

tic                                 %程序计时开始
clear;
clc;

n=0:1/255/1000:1/1000;              %时间长度1ms,采样频率256KHz

w1=2000*pi
s1=sin(w1*n);
subplot(3,1,1)
plot(s1)		                     %输入正弦信号1,频率f1=1KHz

w2=4000*pi
s2=sin(w2*n)
subplot(3,1,2)
plot(s2)		                     %输入正弦信号2,频率f2=2KHz

S=[s1;s2];                           %输入信号矩阵

cita1=pi/6;                          %信号1方位角30度
cita2=3*pi/8;                        %信号2方位角135度

M=[0:5]			                     %阵元数M=6
ratio=1/2                            %波长间距比

a1=exp(j*2*pi*ratio*M*cos(cita1))
a2=exp(j*2*pi*ratio*M*cos(cita2))
 A=[a1.' a2.'];		                 %方位角矩阵

N=randn(length(M),length(n))    	 %白噪声矩阵

X=A*S+1.5*N			                 %阵元输出信号

Rx=X*X'/length(s1)		             %求自相关矩阵
    
[V,D]=eig(Rx)			             %自相关矩阵的特征分解,,得到特征值矩阵和与特征值一一对应的特征向量矩阵

for i=1:length(M)                    %对特征值排序
    P(i)=D(i,i)
end

[Q,index]=sort(P)                    %将与特征值对应的特征向量也按照特征值的顺序排列

for i=1:length(M)-2                  %取排了序的特征向量的前4个来构造噪声子空间En
    En(:,i)=V(:,index(i))
end

a=[0:180]                           
m=M'
a_cita=exp(j*2*pi*m*ratio*cos(pi*a/180))

temp=En'*a_cita

for b=1:180;
    fanshu=norm(temp(:,b));
    out(b)=1/fanshu^2;		        %进行谱峰搜索
end;

subplot(3,1,3)
plot(out)

toc                                 %%程序计时结束