基于MUSIC算法的信号方位估计

up目录

一、理论基础

二、核心程序

三、测试结果

一、理论基础

        在阵列信号处理的许多应用中,需要准确估计空间信号源的方向及空间分布,通常称为“空间谱”。传统的处理方法是利用天线波束作空域扫描,其缺点是分辨能力受到由阵列孔径大小决定的所谓瑞利限的限制。一旦两个信号源处于波束之内时便无法分辨。后来提出“超角分辨”技术,即可以突破瑞利限的限制,实现对处于同一波束内的信号源的分离。在这些方法中,比较典型的有基于空间协方差矩阵特征值分解的一类算法,如Ｓｃｈｍｉｄｔ提出的多重信号分类(ＭＵＳＩＣ)法、ｍｉｎ ｎｏｒｍ法以及Ｒｏｙ提出的基于旋转不变技术的信号参数估计方法(ＥＳＰＲＩＴ)等。它们利用空间协方差矩阵的特征向量来构造信号子空间与噪声子空间,由于它们相互正交或信号矢量经旋转后空间参数不变,因此可确定信号的波达方向。

       对于有相干信号存在的情况,信号矢量将有可能落入噪声子空间中,导致空间协方差矩阵发生秩亏损,因此在这种情况下,基于空间协方差矩阵特征值分解的方法将会失效。而多径相干信号在雷达、通信、声纳信号处理的应用中是很常见的。如在雷达测高或低角跟踪的应用中,目标直接回波与地面反射波是强相关的;类似的例子还有通信中基站与移动台之间的信号传输,因此,在多径信号存在的情况下,如何进行高分辨处理是一个重要的研究课题。目前已提出不少方法来解决这一问题。但这些方法大多是基于空间平滑技术来纠正协方差矩阵,然后应用ＭＵＳＩＣ等正交化方法进行处理。当然高分辨技术还有最大似然(ＭＬ)方法等,但其运算量过大,难以实时实现,故这里不作考虑。综上所述,上述方法都是基于空间协方差矩阵的特征值分解来构造信号子空间与噪声子空间的,只不过是采用空间平滑或降维等措施来解决空间协方差矩阵的秩亏损问题。针对以上情况,本文也从另一个角度出发,寻求另一种矩阵,这种矩阵同样具有空间谱的特征,但却不受信号相关性的影响,从而可对多径相关信号作出正确的分离。

       MUSIC算法是一种基于矩阵特征空间分解的方法。从几何角度讲，信号处理的观测空间可以分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间则由协方差矩阵中所有最小特征值（噪声方差）对应的特征向量组成。
        假设在接收器处有一个带有M个天线的线性阵列。在天线阵列处以入射角θ \thetaθ接收入射信号，天线之间的间距d是信号的半波长。信号在不同的天线上具有不同的传播长度，由于传播路径比天线间隔d长得多，因此相邻天线之间的路径差可以表示为d s i n ( θ ) dsin(\theta)dsin(θ) 。因此相邻天线之间引入了− 2 π f d s i n ( θ ) / c -2 \pi fdsin(\theta)/c−2πfdsin(θ)/c的相位差，其中f是信号频率，c是光速。

       多重信号分类( Multiple Signal Classification，MUSIC)算法的基本原理是把阵列输出数据的协方差矩阵进行特征分解，得到与信号分量相对应的信号子空间和与信号分量正交的噪声子空间，然后利用两个子空间的正交性实现信号的入射方向估计。

 综上所述，MUSIC算法的计算步骤：

（1）利用阵列所接收的L快拍数据，可以得到数据的协方差矩阵；
（2）对进行特征值分解，由的特征值就可以判断信号源数；
（3）得到信号子空间与噪声子空间；
（4）根据式子进行谱峰搜索；
（5）找出极大值所对应的角度就是入射信号的角度。

       MUSIC（Multiple Signal Classification）算法的方法类似于MVDR算法，只是在最后计算的时候有些许不同。该算法建立在一下假设基础上：   ①、噪声为高斯分布，每个麦克风之间的噪声分布为随机过程，之间相互独立，空间平稳，即每个麦克风噪声方差相等。 ②、要求空间信号为高斯平稳随机过程，且与单元麦噪声互不相干，相互独立。 ③、目标声源数目小于阵列元数目。 ④、阵元间隔不大于来波信号中频率最高波长的1/2。   MUSIC算法是空间谱估计技术的代表之一，它利用特征结构分析。其基本原理是将协方差矩阵进行特征值分解。它通常把空间信号分为两种，一种是信号特征向量对应的信号空间，另一种是噪声向量对应的噪声空间，利用噪声空间和信号空间的正交性原理，构造空间谱函数进行搜索，从而预估出DOA信息。

二、核心程序

function [beam_data ,time_openwindow,the_data]=beamforming(n1,n2,angle,T)[time_openwindow,raw_data]=readit(T);
Timeslice=size(raw_data,2);
for k=n1:1:n2%截取要处理的数据the_data(k-n1+1,:)=raw_data(k,:);
endj=sqrt(-1);
N=n2-n1+1;%初始化变量
f=300000;
c=1500;
d=c/f/2;
angle=angle*pi/180;
beam_data=zeros(1,Timeslice);
for t=1:1:Timeslice%在每一个时间片数上做处理s=hilbert(the_data(:,t));for k=1:N%求没个元上的插入相移v(k)=exp(-j*2*pi*d*sin(angle)*f*(n1-1)/c)*exp(-j*2*pi*d*sin(angle)*f*(k-1)/c);end%---------------加窗处理——------------------
A=0.54;
B=0.46;
C=0;
for n=1:1:N
w(n)=A-B*cos(2*pi*(n-1)/(N-1))+C*cos(4*pi*(n-1)/(N-1));
endfor k=1:N%将各个阵元上的相加beam_data(t)=beam_data(t)+w(k)*s(k)*v(k);end
end
up42

三、测试结果