ula_azimuth_fre.m

基于均匀线阵的

function ULA_Azimuth_Fre  %用MUSIC算法和MNM算法对均匀线阵方位角与频率进行联合估计
                          %%%方位角是信号与线阵的夹角
clc;
clear;
close all;
M=16;                             %阵元数目
N=4;                              %信源数目
c=1500;                           %波速度
f=30000;                           %工作频率
lamda=c/f;                        %波长
d=0.5*lamda;                      %阵元间距
loop=3;                           %延迟级数
snap=1000;                        %快拍数目
format short
tao=1/3/f;
ratio=0.5;
fs=1000;
ts=1/fs;
t=(0:snap-1)/fs;
theta_in=[40,0.3;50,0.3;80,0.2;80,0.1]; %输入信号

X=Signal(t,M,N,snap,ratio,theta_in,20,f,loop,tao);
       Rx=X*X'/snap;
      [U,R]=eig(Rx);
      size(Rx);
for iii=1:M*loop
    b(iii)=abs(R(iii,iii));
end
[c,e]=sort(b);
for iii=1:M*loop-N
U_n(:,iii)=U(:,e(iii));
end
c_1=U_n(1,:);      %1X44
E_n=U_n(2:48,:);   %47X44
d_1=[1 ;(E_n*c_1')/(c_1*c_1')];
%theta=0:1:180;
Theta=0:180;
 omig=0:0.01:0.5;
  %a=[1:M]';
  aa=[0:(loop-1)]';
for ii=1:length(Theta)
    for jj=1:length(omig)
        a_theta=[];
         for jjj=1:M
         Temp=exp(j*[2*pi*ratio*omig(jj)*(jjj-1)*cos(Theta(ii)/180*pi)+aa*2*pi*f*omig(jj)*tao]);
         a_theta=[a_theta;Temp];
         end 
       P_MUSIC(ii,jj)=10*log10(1/abs(a_theta'*U_n*U_n'*a_theta));
       P_MNM(ii,jj)=10*log10(1/abs(a_theta'*d_1*d_1'*a_theta));
 end
end
[X,Y]=meshgrid(omig,Theta);
figure(1);
mesh(X,Y,P_MUSIC);grid on;xlabel('归一化频率');ylabel('方向角');title('均匀线阵频率与方位角联合估计—MUSIC算法');
figure(2);
mesh(X,Y,P_MNM);grid on;xlabel('归一化频率');ylabel('方向角');title('均匀线阵频率与方位角联合估计—MNM算法');
return
% %========================================信号源====均匀线阵====================
% %M阵元数 ,N信源数 ,snap快拍数 ,R阵元间距与波长比 ,theta 信源方位角  ,SNR 信噪比
function out=Signal(t,M,N,snap,R,theta,SNR,f0,loop,tao)   
 %a=[0:(M-1)]';
 b=[0:(loop-1)]';
for ii=1:N
 S(ii,:)=exp(j*2*pi*(f0*theta(ii,2)*t+1*2^(ii-1)*t.^2));
end
 for ii=1:N
     A1=[];
     for jj=1:M
         Temp=exp(j*[2*pi*R*theta(ii,2)*(jj-1)*cos(theta(ii,1)/180*pi)+b*2*pi*f0*theta(ii,2)*tao]);
         A1=[A1;Temp];
     end 
       A(:,ii)=A1;
 end
  X0=A*S;
   randn('state',0);
    real_noise=randn(size(X0));
     randn('state',3);
      imag_noise=randn(size(X0));
       noise0=(real_noise+j*imag_noise)/2^0.5;
         noise=10^(-SNR/20)*noise0;
    out=X0+noise;
return