music_3_26.m

信号DOA的估计算法,学习空间谱估计的同学可以参考

%--------------------------------------------------------------------------------------------
%/////////////////////         MUSIC Algorithm           /////////////////////////////////////
%---------------------------------------------------------------------------------------------
%assumptions:*elements is the number of array elements,
%            *targets is the number of incident signals ,and elements is biger than targets.
%            *The vector a(w)corresponding to different values of w are linearly independent.
%            *a(w) depends on the array elements , ans its response to a signal incident from
%             the direction of the signals.
%            *The noise is mean-zero Gaussian White noise and its variance is b2.
%            %And the incident signals and the noise are uncorrelated.
%----------------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------------
%At the begining of the algorithm, creat some essential data, 
%and the data are created by the function produce_data. 
%Senondly, calculateing the covariance R of the created data,
%and sometimes the covariance R is estimated by snapshots.
%after we get the data of the covariance we can know the eigen-structure of R
%e_value is the eigenvalues of R, and e_vector is the eigenvecter corresponding to different 
%eigenvalues. The number of eigenvalues is just the same as elements.
%At last we structure the spectrum-function and estimate the DOAs for the spectrum-function.
%-----------------------------------------------------------------------------------------------
%input the parameters of the function produce_date, in order to get the signals
tic

%targets is the number of the wanted signal
targets=5
%elements is the number of array element
elements=8
%snapshots is the the snapshot of signal received by array
snapshots=1000
%dl is the distance between two elements of the array
dl=0.5

%ang is the direction-of-arrival of the wanted signal
ang=[-45,-42,0,15,45]

%vn is the signal-noise-ratio of the signal received by array
%vn is given with dB
vn=20
snr1=ones(1,targets)
snr=snr1.*vn


%produce_data is the function which produce the data as if the 
% data is received by array.
ul=produce_data(targets,elements,snapshots,dl,snr,ang);

%to get the covariance matrix R by the original data
R=ul*ul';
%cover=cov(ul');
R=R/snapshots;

%Counting the eigen-structure of the covariance R
[e_vector,e_value]=eig(R) 

%Pick up the diagonal elements of the matrix e_value in order to compare the big and small
e_value=diag(e_value)

%Compare the big and the small between the eigenvalues.
%And at the same time adjust the eigenvector corresponding to the eigenvalue.
for num1=1:elements
    for num2=num1+1:elements
        if(e_value(num1)<e_value(num2))
         chan_value=e_value(num1);
         e_value(num1)=e_value(num2);
         e_value(num2)=chan_value;
         chan_vector=e_vector(:,num1);
         e_vector(:,num1)=e_vector(:,num2);
         e_vector(:,num2)=chan_vector;
      end
   end
end

%Form the noise subspace
if (targets<elements)
   noise_vector=e_vector(:,(targets+1):elements);
else 
      noise_vector=0

end

%Then construct the spectrum function 
for num1=1:3001
    num_of_angle1(num1)=((-pi/2)/3000)*(3001-num1);
   %num_of_phase(num1)=2*pi*dl*sin(num_of_angle(num1));
end

for num2=1:3000
   num_of_angle2(num2)=(pi/2/3000)*(num2);
   %num_of_phase(num2)=2*pi*dl*sin(num_of_angle(num2));
end

num_of_angle=[num_of_angle1 num_of_angle2];
num_of_phase=2*pi*dl*sin(num_of_angle);

for num=1:elements
    dir_vector(num,:)=exp((num-1)*num_of_phase*j);
end

for num1=1:6001
    p_music(num1)=1/norm((noise_vector'*dir_vector(:,num1)),2);
end

%Solve the direction-of-arrival,and the angle corresponding to the central maximum of the 
%spectrum function is the DOA.
%The count_num counts the number of the incident signal
%The DOA is the angle of the incident signal
count_num=0;
for num=2:6000
   if(p_music(num)>p_music(num-1)&p_music(num)>p_music(num+1)&(20*log(p_music(num))>0))
      count_num=count_num+1;
      DOA(count_num)=num_of_angle(num)*180/pi
      DOA_num(count_num)=num
   end
end

%Draw the graph of the angle and the specturm function.
plot(num_of_angle*180/pi,20*log(p_music),'k');
xlabel('方向角（度）') 
ylabel('功率谱(dB)')

%count the Cramer-Rao Bound of the estimation
%count the direction Matrix
%DOA=DOA*pi/180
%for num1=1:targets
 %  for num2=1:elements
 %     A(num1,num2)=exp((num2-1)*2*pi*dl*sin(DOA(num1))*j);
  % end
%end

%eff=VARmu/VARcr
%Aa=A*A'
%Az=inv(Aa)
%for num=1:targets
%   eff(num)=1+Az(num,num)/snr(num)
%end

toc
t=toc