%BEAMSPACE DOA估計
clc;
clear;
%參數設置
lookdirection=[1 5 30];%信源到達角矢量
source_state=[5673163 2348321 435937417];%生成信源的隨機數初始狀態
source_amp=[1 1 1];%信源幅度

noise_state=2374357 ;%生成噪聲的隨機數初始狀態
noise_power=-10;%噪聲功率
num_element=16; %陣元數
num_V_noise=14; %噪聲矢量個數

distance=0.5;%陣元間距，以波長為參考，
num_snapshot=5000; %快拍數
 
beamwlookdirection=[90+rad2deg(asin(1/16)) 90+rad2deg(asin(3/16))  90+rad2deg(asin(5/16))];
beamw=SteerVectGen(beamwlookdirection,num_element,distance);
%收信號
steervect=SteerVectGen(lookdirection,num_element,distance);%導引矢量
num_source=length(lookdirection);%信源數
rec_signal=0;
for index=1:num_source
source=source_amp(index)*randsrc(1,num_snapshot,[-1 1 ],source_state(index));
rec_signal=rec_signal+kron(source,steervect(:,index));
end
rec_signal=rec_signal+wgn(num_element,num_snapshot,noise_power,1,noise_state,'complex');%加入噪聲

%MUSIC算法
%求相關矩陣
R=0;
for index=1:num_snapshot
    snap=rec_signal(:,index);
    R=R+snap*snap';
end
R=R/num_snapshot;

if 1
%波束域
R=beamw'*R*beamw;
num_v_noise=1;
[v,dl]=eig(R);%求特征矢量
v_noise=v(:,1:num_v_noise);
uuh=v_noise*v_noise';
%求
theta=0:0.05:180;
scan_marix=SteerVectGen(theta,num_element,distance);
for index=1:length(theta)
    temp=scan_marix(:,index);
    temp1=beamw'*temp;
    temp2=(temp1'*temp1)/(temp1'*uuh*temp1)/num_element;
    p(index)=10*log10(abs(temp2));
    end
else
%陣元域
    [v,dl]=eig(R);%求特征矢量
    v_noise=v(:,l:num_v_noise);
uuh=v_noise*v_noise';
%求
theta=0:0.05:180;
scan_marix=SteerVectGen(theta,num_element,distance);
for index=l:length(theta)
    temp=scan_marix(:,index);
    tempt1=beamw'*temp;
    temp2=(tempt'*templ)/(tempt1'*uuh*templ)/num_element;
    p(index)=10*log10(abs(temp2));
end
end

%顯示空間譜
plot(theta,p);
grid on;