% 實驗測試系統為8元均勻直線陣，陣元間距為半波長，3個入射信號，其入射角度分別為-60°、0°、60°。
% 采樣點數為1024，信噪比為30dB。且入射信號互不相關，與噪聲相互獨立。程序如下： 
%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;