yanchi_xiangjia.m

DOA估计中的延迟相加法

function p=yanchi_xiangjia(angle)
clear
clc


%假定信号数已确定为i=3个
        %  1）入射信号为相干信号，相对时间暂定为Tau0=0、Tau1=3、Tau2=5；
        %  2）入射方向可以随机产生。此处暂定入射方向为：theta1=5度，theta2=30度，theta3=45度；
theta1=10*pi/180;theta2=40*pi/180;theta3=60*pi/180;

%空间平滑
        %  1）阵列为ULA，阵元数暂定为Narray=8；
        %  2）子阵长度暂定为p=4；


%产生发送信号和噪声序列。
%注意：此时未考虑相互之间的时延！信号之间为相干；
    
    
    N=1000;
    bit1=2*rem(unidrnd(2,1,N),2)-1;
    bit2=2*rem(unidrnd(2,1,N),2)-1;
    bit3=2*rem(unidrnd(2,1,N),2)-1;
    sig=[bit1;bit2;bit3];
   
    
    SNRdB=-5;%暂定一个信噪比
    n0=1/(10^(SNRdB/10))%*(3/8);%3个相干信号的总功率/8个阵元的总的噪声功率＝5dB    
    sigma=sqrt(n0);
    
    noise00=randn(1,N);
        noise0=sigma*noise00;
    noise11=randn(1,N);
        noise1=sigma*noise11;
    noise22=randn(1,N);
        noise2=sigma*noise22;
    noise33=randn(1,N);
        noise3=sigma*noise33;
    noise44=randn(1,N);
        noise4=sigma*noise44;
    noise55=randn(1,N);
        noise5=sigma*noise55;
    noise66=randn(1,N);
        noise6=sigma*noise66;
    noise77=randn(1,N);
        noise7=sigma*noise77;
    
        noise=[noise0;noise1;noise2;noise3;noise4;noise5;noise6;noise7];
    
    %阵列流形
    %aa=[1,                     1,                     1;
    %    exp(-j*pi*sin(theta1)),exp(-j*pi*sin(theta2)),exp(-j*pi*sin(theta3));
     %   exp(-j*2*pi*sin(theta1)),exp(-j*2*pi*sin(theta2)),exp(-j*2*pi*sin(theta3));
      %  exp(-j*3*pi*sin(theta1)),exp(-j*3*pi*sin(theta2)),exp(-j*3*pi*sin(theta3));
       % exp(-j*4*pi*sin(theta1)),exp(-j*4*pi*sin(theta2)),exp(-j*4*pi*sin(theta3));
        %exp(-j*6*pi*sin(theta1)),exp(-j*6*pi*sin(theta2)),exp(-j*6*pi*sin(theta3));
        %exp(-j*7*pi*sin(theta1)),exp(-j*7*pi*sin(theta2)),exp(-j*7*pi*sin(theta3))];
    
    a1=[1;                   
        exp(-j*pi*sin(theta1));
        exp(-j*2*pi*sin(theta1));
        exp(-j*3*pi*sin(theta1));
        exp(-j*4*pi*sin(theta1));
        exp(-j*5*pi*sin(theta1));
        exp(-j*6*pi*sin(theta1));
        exp(-j*7*pi*sin(theta1))];
   a2=[1;                   
        exp(-j*pi*sin(theta2));
        exp(-j*2*pi*sin(theta2));
        exp(-j*3*pi*sin(theta2));
        exp(-j*4*pi*sin(theta2));
        exp(-j*5*pi*sin(theta2));
        exp(-j*6*pi*sin(theta2));
        exp(-j*7*pi*sin(theta2))];
     a3=[1;                   
        exp(-j*pi*sin(theta3));
        exp(-j*2*pi*sin(theta3));
        exp(-j*3*pi*sin(theta3));
        exp(-j*4*pi*sin(theta3));
        exp(-j*5*pi*sin(theta3));
        exp(-j*6*pi*sin(theta3));
        exp(-j*7*pi*sin(theta3))];
        Rx=zeros(8,8);
        p=zeros(1,90);%绘图用
for m=1:1000
    
    %阵列的噪声协方差矩阵
    n=diag([noise0(m)*noise0(m),noise1(m)*noise1(m),noise2(m)*noise2(m),noise3(m)*noise3(m),...
        noise4(m)*noise4(m),noise5(m)*noise5(m),noise6(m)*noise6(m),noise7(m)*noise7(m)]);
    
    %阵列的协方差矩阵
            %注意：假设条件为信号和噪声是相互独立的
    
    x=a1*sig(1,m)+a2*sig(2,m)+a3*sig(3,m)+noise(:,m);
    xx=x*x';
    Rx=((m-1)*Rx+xx)/m;
end
  
  %Capon空间谱
  angle=0:0.1:90;
  for i=1:length(angle)
      theta=((i-1)/10)*pi/180;
      Atheta=[1,exp(-j*pi*sin(theta)),exp(-j*2*pi*sin(theta)),exp(-j*3*pi*sin(theta)),...
              exp(-j*4*pi*sin(theta)),exp(-j*5*pi*sin(theta)),exp(-j*6*pi*sin(theta)),exp(-j*7*pi*sin(theta))].';
     
     
       Pcapon=Atheta'*Rx*Atheta;
        p(i)=sqrt(Pcapon*Pcapon');
      % p(i)=log(p(i));
  end  
  %angle=0:0.1:90;
  plot(angle,p,'b--')%按快拍计算次数显示music谱峰
  %结果显示。正确
  title('延迟-相加法')
  xlabel('角度')
  ylabel('接收信号（db）')
  hold on