FORWARD_SMOOTHNESS_MUSIC_toeplitz_matlab

%TOEPLITZ_MUSIC ALOGRITHM %DOA ESTIMATION BY TOEPLITZ_MUSIC ALOGRITHM clear all; close all; clc; source_number=2;%信元数 sensor_number=8;%阵元数 N_x=1024; %信号长度 snapshot_number=N_x;%快拍数 w=[pi/4 pi/4].';%信号频率 l=[(2*pi*3e8)/w(1) (2*pi*3e8)/w(2)];%信号波长 d=[0.5*l(1) 0.5*l(2)];%阵元间距 snr=10;%信噪比 source_doa=[45 -60];%两个信号的入射角度 A=[exp(-j*(0:sensor_number-1)*d(1)*2*pi*sin(source_doa(1)*pi/180)/l(1));exp(-j*(0:sensor_number-1)*d(2)*2*pi*sin(source_doa(2)*pi/180)/l(2))].';%阵列流型 s=sqrt(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'/snapshot_number; %对数据协方差矩阵斜对角线上的元素进行平均,即Toeplitz化 dd=zeros(2*sensor_number-1); for i=-(sensor_number-1):(sensor_number-1) c=sum(diag(R,i))/(sensor_number-abs(i)); dd(i+sensor_number)=c; end for k=1:sensor_number R(k,k)=dd(sensor_number); end for k=1:(sensor_number-1) R(k+1,k)=dd(sensor_number-1); end for k=1:(sensor_number-2) R(k+2,k)=dd(sensor_number-2); end for k=1:(sensor_number-3) R(k+3,k)=dd(sensor_number-3); end for k=1:(sensor_number-4) R(k+4,k)=dd(sensor_number-4); end for k=1:(sensor_number-5) R(k+5,k)=dd(sensor_number-5); end for k=1:(sensor_number-6) R(k+6,k)=dd(sensor_number-6); end for k=1:(sensor_number-7) R(k+7,k)=dd(sensor_number-7); end for k=1:(sensor_number-1) R(k,k+1)=dd(sensor_number+1); end for k=1:(sensor_number-2) R(k,k+2)=dd(sensor_number+2); end for k=1:(sensor_number-3) R(k,k+3)=dd(sensor_number+3); end for k=1:(sensor_number-4) R(k,k+4)=dd(sensor_number+4); end for k=1:(sensor_number-5) R(k,k+5)=dd(sensor_number+5); end for k=1:(sensor_number-6) R(k,k+6)=dd(sensor_number+6); end for k=1:(sensor_number-7) R(k,k+7)=dd(sensor_number+7); end disp('R'); disp(R); %对Toeplitz化后的协方差矩阵进行特征分解，得到信号子空间和噪声子空间 [V,D]=eig(R); D=diag(D); disp(D); Un=V(:,1:sensor_number-source_number); Gn=Un*Un'; disp('Gn'); disp(Gn); %利用MUSIC算法进行谱峰搜索，从而得到波到达方向 searching_doa=-90:0.1:90;%线阵的搜索范围为-90~90度 for kk=1:length(searching_doa) a_theta=exp(-j*(0:sensor_number-1)'*2*pi*d*sin(pi*searching_doa(kk)/180)/l); Pmusic(kk)=1./abs((a_theta)'*Gn*a_theta); end plot(searching_doa,10*log(Pmusic)); %axis([-90 90 -90 90]); xlabel('入射角/度'); ylabel('谱峰'); legend('TOEPLITZ-MUSIC SPECTRUM'); title('拖布列兹MUSIC估计'); grid on;