%構造信號源
n=-pi:0.01:pi;                           %范圍    
x1=100*cos(50*sqrt(2.0)*pi*n+pi/3);      %構造三個余弦信號
x2=100*cos(50*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];     
cys=600;   %采樣數 
SN1=20;SN2=10;SN3=12;  %信噪比
SN=[SN1;SN2;SN3];
degrad=pi/180;  %將角度制化為弧度制    

%信號和噪聲
load shy3    %載入信號源
TT=1:cys;    %采樣數矩陣
S=[x1(TT).';x2(TT).';x3(TT).']; 
nr=randn(m,cys);  %建立高斯隨機序列
ni=randn(m,cys);
E=nr+ni*j;        %建立高斯噪聲
Ps=S*(S/cys)';    %方向向量
Ps1=diag(Ps);
refp=2*10.^(SN./10);
tmp=sqrt(refp./Ps1);
S1=diag(tmp)*S;

%計算協方差矩陣，特征值分解
tmp=-i*pi*sin(angle'.*degrad);  %相位差
tmp1=[0:m-1]';        %信號到達陣元的相位差組成的向量
A=tmp1*tmp;
A1=exp(A);            %m*p維方向矩陣
X=(A*S1)+E;           %第k個陣元上的接收信號
Rxx=X*(X/cys)';       %接收信號的協方差
[E,S,V]=svd(Rxx);     %計算特征值
Vs=E(:,1:p);          
Vu=E(:,p+1:m);

%計算空間譜函數
th2=[-90:90]';         %線性信號角度范圍
tmp2=-i*pi*sin(th2'.*degrad);  %[-90 90]范圍內的相位差組成的向量
tmp3=tmp1;
A2=tmp2*tmp3;       
A3=exp(A2);            %m*p維方向矩陣
num=diag(A3'*A3);      
Ena=Vu'*A3;
den=diag(Ena'*Ena);     
doa=num./den;          %MUSCI估計值

%做出空間譜圖
semilogy(th2,doa);     %y軸對數坐標曲線
title('MUSIC算法');    %圖表標題，x軸y軸標簽及范圍
xlable('角度（弧度制）');
ylable('值');
axis([-90 90 0.1 1e4]);