direct_kim.m

Fast DOA Estimation Algorithm using Pseudo Covariance Matrix的程序

%DIRECT DATA DOMAIN
clear all;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%初始条件
source_number=3;%信元数
sensor_number=51;%阵元数
snapshot=5;%快拍数
w=[3*pi/5 pi/5 2*pi/5].';%信号频率
m=(sensor_number+1)/2;
snr=20;%信噪比
source_doa=[2.3 8.7 14.8];%信号的入射角度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%信号产生
A=[exp(-j*(0:sensor_number-1)*pi*sin(source_doa(1)*pi/180));exp(-j*(0:sensor_number-1)*pi*sin(source_doa(2)*pi/180));exp(-j*(0:sensor_number-1)*pi*sin(source_doa(3)*pi/180))].';%阵列流型
s=sqrt(10.^(snr/10))*exp(j*w*(0:snapshot-1));%仿真信号
x=A*s+(1/sqrt(2))*(randn(sensor_number,snapshot)+j*randn(sensor_number,snapshot));%加了高斯白噪声后的阵列接收信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%接收数据重构
searching_doa=0:0.5:20;%观测方向的变化
a_pseudo=zeros(m);
   xx=x(1:m,:);
for tt=1:length(searching_doa)   %搜索循环
    a=exp(-j*(0:sensor_number-1)*pi*sin(searching_doa(tt)*pi/180));
  for ff=1:m  
    for gg=1:m
        a_pseudo(ff,gg)=a(1,ff+gg-1);
    end
end
for ss=1:snapshot       %快拍循环
x1=zeros(sensor_number,1);
for kk=1:sensor_number
    x1(kk,1)=x(kk,ss);
end  
x_pseudo=zeros(m);
for ii=1:m
    for jj=1:m
        x_pseudo(ii,jj)=x1(ii+jj-1,1);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 奇异值分解，并按大小排序
      [v,u] = eig(x_pseudo,a_pseudo);
      [val,index] = min(abs(diag(u)));
      v_min = v(:,index);    %最小特征值对应的特征向量
        a1=exp(-j*(0:m-1)*pi*sin(searching_doa(tt)*pi/180));  %操纵矢量重构
        w=v_min./(a1*v_min);  
         ww(:,tt)=w  %加权
      y=w.'*x_pseudo(:,ss);
      yy(ss,tt)=abs(y*y.');
end
Pb(tt)=10*log10(sum(yy(:,tt))/snapshot);
end

      [val1,index1] = min(abs(Pb));
w1=ww(:,index1);
searching_doa_j=0:0.1:20;
for ee=1:length(searching_doa_j)
    a2=exp(-j*(0:m-1)*pi*sin(searching_doa_j(ee)*pi/180));
    Pd(ee)=10*log10(1./abs((a2*w1)*(a2*w1).'))
end
    
    
figure(1)
plot(searching_doa,Pb,'r');
%axis([-90 90 -90 90]);
xlabel('入射角/degree');
ylabel('空间谱/dB');
legend('近似方向谱 ');
title('直接数据域谱估计');
grid on;
figure(2)
plot(searching_doa_j,Pd,'r');
%axis([-90 90 -90 90]);
xlabel('入射角/degree');
ylabel('空间谱/dB');
legend('精确方向谱 ');
title('直接数据域谱估计');
grid on;