root-MUSIC算法

%2017.6.15
%ROOT_MUSIC ALOGRITHM
%DOA ESTIMATION BY ROOT_MUSIC
clc;
clear all;
close all;

K=2;                                                             %信源数
M=8;                                                             %阵元数
L=200;                                                           %信号长度
w=[pi/4 pi/6].';                                                 %信号频率
lamda=((2*pi*3e8)/w(1)+(2*pi*3e8)/w(2))/2;%信号波长  
d_lamda=0.5;                                                   　%阵元间距
snr=20;                                                          %信噪比
theta1=[45,60];                                                 %信号入射角
for k=1:K
A(:,k)=exp(-1j*[0:M-1]'*d_lamda*2*pi*sin(theta1(k)*pi/180));　　　%阵列流型
end
for kk=1:L
s(:,kk)=sqrt(10.^((snr/2)/10))*exp(1j*w*(kk-1));%仿真信号
end
x=A*s+(1/sqrt(2))*(randn(M,L)+1j*randn(M,L));%加入高斯白噪声
R=(x*x')/L;                                                       %协方差矩阵
%%%%%%第一种方法%%%%%%%%%%
[V,D]=eig(R);                                                     %对协方差矩阵进行特征分解
Un=V(:,1:M-K);                                                  　%取噪声子空间
Gn=Un*Un';
a = zeros(2*M-1,1)';        　　　                                 %找出多项式的系数，并按阶数从高至低排列
for i=-(M-1):(M-1)
    a(i+M) = sum( diag(Gn,i) );
end
a1=roots(a);                                                      %使用ROOTS函数求出多项式的根                            
a2=a1(abs(a1)<1);                 　　　                          %找出在单位圆里且最接近单位圆的N个根
[lamda,I]=sort(abs(abs(a2)-1));         　　　　　　             %挑选出最接近单位圆的N个根
f=a2(I(1:K));                                                    %计算信号到达方向角
source_doa=[asin(angle(f(1))/pi)*180/pi asin(angle(f(2))/pi)*180/pi];
source_doa=sort(source_doa);
disp('source_doa');
disp(source_doa);
%%%%%%%第二种方法%%%%%%%%%
% [V,D]=eig(R);
% Un=V(:,1:M-K);                                            　　　%(4.5.7）
% Un1=Un(1:M,:);
% Un2=Un(K+1:M,:);
% T=[1 0 0 0 0 0]';                                           　　 %（阵元数-信源数）*1
% c=Un1*inv(Un2)*T;                                     　　　　  % (K*(M-K))*((M-K)*(M-K))*((M-K)*1)=K*1
% c=[1,c(2,1),c(1,1)];                                  　　　     %（4.5.8）
% f=roots(c);　　　　　　　　　　　　　　　　　　　　　　　　
% source_doa=[asin(angle(f(1))/pi)*180/pi asin(angle(f(2))/pi)*180/pi];
% source_doa=sort(source_doa);
% disp('source_doa');
% disp(source_doa);