music.doc

music算法的M文件

%music.m 
clear 
clear all 
%构造信号源 
n=-pi:2.0/1024*pi:pi; %范围 
s1=100*randn(1)*cos(100*randn(1)*sqrt(2.0)*pi*n+pi*randn(1)); %构造三个余弦信号 
s2=100*randn(1)*cos(100*randn(1)*sqrt(2.0)*pi*n+pi*randn(1)); 
s3=100*randn(1)*cos(100*randn(1)*sqrt(2.0)*pi*n+pi*randn(1)); 
save sig3 s1 s2 s3 
%天线阵元及接收信号 
i=sqrt(-1); 
j=i; 
m=8;%天线阵元数 
p=3; 
angle1=-60;%从-60 0 60入射的信号源 
angle2=0; 
angle3=60; 
th=[angle1;angle2;angle3]; 
nn=1024;%采样点数 
SN1=3;%信噪比 
SN2=4; 
SN3=8; 
sn=[SN1;SN2;SN3]; 
degrad=pi/180; 
load sig3 
tt=1:nn; 
S=[s1(tt);s2(tt);s3(tt)]; 
nr=randn(m,nn); 
ni=randn(m,nn); 
U=nr+j*ni; %建立高斯噪声 
Ps=S*S'/nn; 
ps=diag(Ps); 
refp=2*10.^(sn/10); 
tmp=sqrt(refp./ps); %refp,tmp,[;;] 
S2=diag(tmp)*S;%diag() then3*3 ,diag^ 
%计算协方差矩阵，特征值分解 
tmp22=-i*pi*sin(th'*degrad); %相位差 
tmp2=[0:m-1]'; %信号到达阵元的相位差组成的向量 
a2=tmp2*tmp22; 
A=exp(a2); %m*p维方向矩阵 
X=A*S2+U; %第k个阵元上的接收信号 
Rxx=X*X'/nn; %接收信号的协方差 
%计算空间谱函数 
th2=[-90:90]' %线性信号角度范围 
tmp33=-i*pi*sin(th2'*degrad); %[-90 90]范围内的相位差组成的向量 
tmp2=[0:m-1]'; 
a2=tmp2*tmp33; 
A2=exp(a2); %m*p维方向矩阵 
[V,D]=eig(Rxx) 
G=V(:,1:m-p) %信号与噪声分离 
for n=1:length(th2) 
doa(n)=1.0/(A2(:,n)'*G*G'*A2(:,n)) 
end 
%做出空间谱图 
semilogy(th2,doa); %y轴对数坐标曲线 
title('MUSIC Spectyum'); %图表标题，x轴y轴标签及范围 
xlabel('Angle(deg)'); 
ylabel('Spectrum'); 
axis([-90 90 0.1 1e3]); 
grid;