%MUSIC算法，7均勻直線陣列 天線間距為波長一半 3信號源
clear all
clc
d=1;
lma=2;
N=2048;                                                  %設定采樣點數
SNR=[4 8 -2];                         
fs=8;                                                    %設定采樣頻
f=[0.9 0.6 0.2];                                         %設三個信號元的頻率，單位為MHz
w=[-30 35 45]*pi/180;                                    %入射方向角向量
phi=[60 20 30]*pi/180;                                   %兩入射信號初始相位
M=7;                                                     %設定陣元數
AMP = 10.^(SNR./20);              
m=0:M-1;
n=0:N-1;
S=repmat(AMP',1,N).*exp(j*(2*pi*f.'/fs*n+phi'*ones(1,N)));             %產生三個遠場正弦信號
A=exp(j*2*pi*d*m'*(sin(w)/lma));                         %方向矩陣
X=A*S;
noise=randn(M,N)+j*randn(M,N);
X=X+noise;                                               %接受數據矩陣
Rx=X*(X')/N;                                             %計算協方差矩陣
%*********************************************************************
[V,D]=eig(Rx);			                                 %自相關矩陣的特征分解,,得到特征值矩陣和與特征值一一對應的特征向量矩陣
[Q,index]=sort(diag(D));                                 %對特征值排序 
for k=1:M-length(w)                                      %將與特征值對應的特征向量也按照特征值的順序排列     
    En(:,k)=V(:,index(k));                               %取排了序的特征向量的前M-length(theta)個來構造噪聲子空間En
end
L_lim=-90;                                               %設置測向左、右界
R_lim=90
step=0.1;                                                %設置步長
theta_axis=L_lim:step:R_lim;
MD=(R_lim-L_lim)/step+1;
a_theta=exp(j*2*pi*d/lma*m'*sin(theta_axis*pi/180));
for k=1:MD                                              %進行譜峰搜索
    d(k)=a_theta(:,k)'*En*En'*a_theta(:,k);
    k=k+1; 
end
D2=-10*log10(abs(d)./max(abs(d)));
figure(1)
polar(theta_axis,D2)
figure(2)
plot(theta_axis,D2,'k');
xlabel('DOA(theta)')
ylabel('P(dB)')
ylim([0,50])
grid