求根MUSIC算法源程序

%ROOT_MUSICALOGRITHM%DOAESTIMATIONBYROOT_MUSICclearall;closeall;clc;source_number=2;%信元数sensor_number=8;%阵元数N_x=1024;%信号长度snapshot_number=N_x;%快拍数w=[pi/4pi/6].';%信号频率l=((2*pi*3e8)/w(1)+(2*pi*3e8)/w(2))/2;%信号波长d=0.5*l;%阵元间距snr=0;%信噪比theta1=45;theta2=60;%两个信号的入射角度A=[exp(-j*(0:sensor_number-1)*d*2*pi*sin(theta1*pi/180)/l);exp(-j*(0:sensor_number-1)*d*2*pi*sin(theta2*pi/180)/l)].';%阵列流型s=10.^((snr/2)/10)*exp(j*w*[0:N_x-1]);%仿真信号%x=awgn(s,snr);x=A*s+(1/sqrt(2))*(randn(sensor_number,N_x)+j*randn(sensor_number,N_x));%加了高斯白噪声后的阵列接收信号R=(x*x')/N_x;%对协方差矩阵进行特征分解[V,D]=eig(R);Un=V(:,1:sensor_number-source_number);Gn=Un*Un';%找出多项式的系数，并按阶数从高至低排列a=zeros(2*sensor_number-1,1)';fori=-(sensor_number-1):(sensor_number-1)a(i+sensor_number)=sum(diag(Gn,i));end%使用ROOTS函数求出多项式的根a1=roots(a);%除掉大于1的增根（方程式在此条件下变形的）a2=a1(abs(a1)1);%挑选出最接近单位圆的N个根[lamda,I]=sort(abs(abs(a2)-1));f=a2(I(1:source_number));%计算信号到达方向角source_doa=[asin(angle(f(1))/pi)*180/piasin(angle(f(2))/pi)*180/pi];source_doa=sort(source_doa);disp('source_doa');disp(source_doa);