%?coherent?DOA?estimation
%?using?classic?MUSIC?
%?d-D?searching?MUSIC?d-D?searching?using?big?eigen?vector?1-D?searching?
%?collus?wang?2009-11-28

clc
clear
close?all

M?=?8;??????????????%?陣元數
theta?=?[70100]*pi/180;???????%?信號入射方向
w?=?2*pi*[3030]*1e3;??????????%?入射信號頻率
Coherent_On?=?1;
d?=?length（theta）;??????????????????%?信號個數
T?=?2e6^（-1）;???????????????????????????%?采用間隔
N?=?5000;????????????????????????????????????%?快拍數

f?=?2.4e9;??????????????????????????%?載波頻率
c?=?3e8;????????????????????????????%?光速
D?=?c/f/2;??????????????????????????%?陣元間距
SNR_dB?=?10;????????????????%?SNR?in?dB

A?=?zeros（Md）;?????????????%?導向矩陣
for?index?=?1?:?d
????A（:?index）?=?[exp（?1j*2*pi*f*[0:M-1]*D*cos（theta（index））/c?）].‘;
end

s?=?zeros（d?N）;????????????%?信號
for?index?=?1?:?d
????s（index:）?=?exp（?1j*w（index）*[0:N-1]*T?）;%+?rand（1N）*2*pi
end

rho?=?ones（1d）;????????????%?相干源相關系數
rho（2:d）?=?exp（1j*rand（1d-1）*2*pi）;
if?（Coherent_On）
????for?index?=?2:d
????????s（index:）?=??s（1?:）.*rho（index）;
????end
end

sigma_noise?=?1/sqrt（10^（SNR_dB/10））;???????%?加噪?以單個信號的信噪比進行計算
X?=?zeros（M?N）;
X?=?A*s?+?sigma_noise*（randn（MN）+1j*randn（MN））/sqrt（2）;???????%?接收信號

%?-------------------?Beginning?of?Algorithms?-------------
R_xx?=?X*X‘/N;??????????????????%?calution?of?R_xx
[UDia]=eig（R_xx）;??????????????%?eigen?value?decomposition
Dia

sum?=?trace（abs（Dia））;??????????%?find?big?eigen?values
temp?=?0;
pct?=?0.90??????????????????????????????????%?if?the?sum?of?biggest?eigen?values?exceed?pct?of?trace?we?say?we?get?the?biggest?values.
for?index??=?M?:?-1?:??1
????temp?=?temp?+?abs（Dia（indexindex））;
????if?temp?>?pct*sum
????????break
????end
end
sig_index?=?index?:?M
noise_index?=?1:?index-1;

res?=?0.5;??????????????????????????????????????%?resolution?of?peak?searching?（in?deg）
theta_search?=?0:?res*pi/180:?pi;

%?-------------?classic?music---------------
P_music_classic?=?zeros（1?length（theta_search）?）;
P?=?U（:?noise_index）*U（:noise_index）‘;
for?index?=?1:length（theta_search）
????E?=?[exp（1j*2*pi*f*[0:M-1]*D*cos（theta_search（index））/c?）].‘;
????P_music_classic（index）?=?1/（E‘*P*E）;
end

figure（10）
plot（theta_search/pi*180?abs（P_music_classic））
hold?on
st