%MUSIC ALGORITHM
%DOA Estimation by ULA
%信號子空間法。
clear all; 
close all;
clc
N_x=1024; % Length of Signal 
N=8; % Size of Rx Matrix 
A=[4 3 3];
M=3; % Number of Signals 
w=[pi/6 pi/10 pi/8]';%信號頻率 與2不同之處
K=0:N-1;
xx=pi/3;yy=0*pi/4;zz=-pi/6; %三個信號的入射角,小于5度就難以辨別了
B=[exp(-j*pi*K'*sin(xx)) exp(-j*pi*K'*sin(yy)) exp(-j*pi*K'*sin(zz))] ; %陣列流型,信號源決定行數，陣元數決定列數
xxx=2*exp(j*w*[0:N_x-1]);%仿真信號
x=B*xxx+randn(8,N_x)+j*randn(8,N_x); %加噪聲
R=x*x';
[V D]=eig(R);
[lambda,index] = sort((diag(D)));
UU=V(:,index(6:8));
theta=-90:0.1:90;%ULA估計角度變化的范圍和頻率選擇
for i = 1:length(theta)
AA=exp(-j*pi*K*sin(theta(i)*pi/180));%方向矢量
WW=AA*(eye(size(R))-UU*UU')*AA';
Pmusic(i)=abs(1/WW);%角譜
% PCap(i) = abs(1/(AA*inv(R)*AA'));
end
Pmusic = 20*log10(Pmusic);
sita=-90:0.1:90;
plot(sita,Pmusic);
grid
xlabel('rad')