toeplitz.m

含有前向空间平滑、经典MUSIC算法、TOEPLITZ算法进行DOA估计的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;