uca_moshi.m

阵列信号处理

function UCA_MOSHI  %========均匀圆阵模式空间算法 by 王方勇 学号：ZZ07050019  2008.4.18===========
                    %========通过模式空间法形成虚拟均匀线阵并用经典MUSIC算法进行谱峰搜索==========
                    %========本程序显示了模式空间算法性能(成功概率、估计偏差、方差)与信噪比的关系=
        
clc;
 clear;
  close all;
   tic
    format short
     SNR=0:30;
      Oknum=[];                          %记录成功次数
       for snr=1:length(SNR)             %分别在0-30dB范围内做实验
           Oknum(snr)=0;
           Estimated_Angle=[];
       for numb=1:100                    %分别做独立实验100次
%===============+参数设定+==========
M=16;                                                        %阵元数目
 c=1500;                                                     %波速度
  f=30000;                                                   %频率
   lamda=c/f;                                                %波长
    snap=100;                                                %快拍数
     ratio=5/pi;                                             %圆阵半径与波长比
      r=ratio*lamda;                                         %圆阵半径
       Beita=2*pi*r/lamda;
        K=4;                                                 %可激发的最大模式数
         fs=100;                                             %快拍采样频率
          ts=1/fs;
           t=(0:snap-1)*ts;
            N=2;                                            %信源数目
             theta_in=[30 60];                              %信号输入角度矢量                   
%==========================================================================
X=Signal(t,M,N,snap,r,lamda,theta_in,SNR(snr),f);           %阵元接收数据
%=======构造(2K-1)*(2K-1)维矩阵===========
     for ii=-K:K
         JJ(ii+K+1)=j^ii*besselj(ii,-Beita);
     end
     J=diag(JJ,0);
%=======构造F矩阵========================
     for ii=-K:K
         for jj=0:M-1
             Fh(ii+K+1,jj+1)=exp(j*2*pi/M*ii*jj);
         end
     end
     F=Fh';
%=======构造T矩阵========================
     T=1/M*inv(J)*Fh;
%=======以下为用MUSIC算法进行谱峰搜索 ===
     Ry=T*X*X'*T'/length(t);;
     [U,R]=eig(Ry);
       for iii=1:2*K+1                                        %对特征值进行排序并得噪声子空间
           b(iii)=abs(R(iii,iii));
       end
       [c,e]=sort(b);
       for iii=1:2*K+1-N
           Un(:,iii)=U(:,e(iii));
       end
Theta=0:360;                                           %360度范围内全方位搜索
  a=[-K:K]';
    for ii=1:length(Theta)
        a_theta=exp(j*a*Theta(ii)*pi/180);                    %虚拟均匀线阵的阵列流型矢量
        Pmusic(ii)=10*log10(1/abs(a_theta'*Un*Un'*a_theta));
    end
    Estimated_Angle(numb,:)=Peak_Seek(Pmusic,Theta,N);
     Error=max(abs(Peak_Seek(Pmusic,Theta,N)-theta_in));
     if Error<3
         Oknum(snr)=Oknum(snr)+1;
     end
    
 end
  Theta_mean=mean(Estimated_Angle);                     %估计均值
    Theta_error(snr)=sum(abs(Theta_mean-theta_in))/2;   %估计偏差
      Theta_var(snr)=sum(var(Estimated_Angle))/2;       %估计方差
     
end
     figure(1);
      plot(SNR,Oknum,'-s');grid on;xlabel('SNR/dB');ylabel('成功概率%');title('均匀圆阵模式间算法成功概率与信噪比的关系');
       figure(2);
      plot(SNR,Theta_error,'-*');axis([0 30 0 10]);grid on;xlabel('SNR/dB');ylabel('估计偏差(度)');title('均匀圆阵模式间算法估计偏差与信噪比的关系');
       figure(3);
      plot(SNR,Theta_var,'-o');axis([0 30 0 10]);grid on;xlabel('SNR/dB');ylabel('估计方差');title('均匀圆阵模式间算法估计方差与信噪比的关系');
%     figure(1);
%      plot(Theta,Pmusic);AXIS([0 360 -15 35]);grid on;xlabel('方位角');ylabel('峰值');title('均匀圆阵模式间算法');
toc
return
%==============角度搜索程序=============
%==============用于搜索谱峰对应的角度===
%===输出前N个峰对应的角度===============
 function out=Peak_Seek(pmusic,p_index,N)
 jj=0;
 mm=1;
 temp=pmusic(1);
 peak=[];
 index=[];
 for ii=2:length(pmusic)
     if temp<pmusic(ii)
         jj=1;
     else
         if jj==1
             peak(mm)=temp;
             index(mm)=p_index(ii-1);
             mm=mm+1;
         end
         jj=0;
     end
      temp=pmusic(ii);
 end
  [peak,in]=sort(peak);
    peak=fliplr(peak);
    in=fliplr(in);
    for zz=1:N
        out(zz)=index(in(zz));
    end 
    [out,mm]=sort(out);
 return

%========================================信号源函数--均匀圆阵=========================
%M 阵元数 ,N信源数 ,snap 快拍数 ,r 圆阵半径,lamda 波长 ,theta 信源方位角,SNR 信噪比,f0信号中心频率
%==========================================================================
function out=Signal(t,M,N,snap,r,lamda,theta,SNR,f0)   
 a=[1:M]';
   for ii=1:N
       S(ii,:)=exp(i*2*pi*(f0*t+0.5*5*2^(ii-1)*t.^2));                      %窄带信号
   end
   for jj=1:N
       A(:,jj)=exp(-i*2*pi*r/lamda*cos((2*pi*(a-1)/M)-theta(jj)/180*pi));   %阵列流型
   end
  X0=A*S;
  out=awgn(X0,SNR);                                                         %加信噪比为SNR的高斯白噪声
return
%==================================程序结束+===================================