您的当前位置：首页正文

(完整word版)MUSIC算法

来源：好兔宠物网

专业综合

课程设计报告

空间谱估计算法

一、设计任务

实现空间谱估计算法，并考察算法性能。

二、方案设计

1）由均匀线阵形式，确定阵列的导向矢量； 2）由阵列导向矢量，对接收信号进行建模仿真； 3）根据多重信号分类算法实现空间谱估计；

4）考察算法性能与信噪比，采样率，观测时间等参数的关系。

三、设计原理

3.1空间谱估计数学模型

空间谱估计就是利用空间阵列实现空间信号的参数估计的一项专门技术。整个空间谱估计系统应该由三部分组成：空间信号入射、空间阵列接收及参数估计。相应地可分为三个空间，即目标空间、观察空间及估计空间，也就是说空间谱估计系统由这三个空间组成，其框图见图1。

图1 空间谱估计的系统结构

对于上述的系统结构，作以下几点说明。

(1)目标空间是一个由信号源的参数与复杂环境参数张成的空间。对于空间谱估计系统，就是利用特定的一些方法从这个复杂的目标空间中估计出信号的未知参数。

(2)观察空间是利用空间按一定方式排列的阵元，来接收目标空间的辐射信号。由于环境的复杂性，所以接收数据中包括信号特征(方位、距离、极化等)和空间环境特征(噪声、杂波、干扰等)。另外由于空间阵元的影响，接收数据中同样也含有空间阵列的某些特征(互耦、通道不一致、频带不一致等)。这里的观察空间是一个多维空间，即系统的接收数据是由多个通道组成，而传统的时域处理方法通常只有一个通道。特别需要指出的是：通道与阵元并不是一一对应，通道是由空间的一个、几个或所有阵元合成的(可用加权或不加权)，当然空间某个特定的阵元可包含在不同的通道内。

(3)估计空间是利用空间谱估计技术(包括阵列信号处理中的一些技术，如阵列校正、空域滤波等技术)从复杂的观察数据中提取信号的特征参数。

从系统框图中可以清晰的看出，估计空间相当于是对目标空间的一个重构过程，这个重构的精度由众多因素决定，如环境的复杂性、空间阵元间的互耦、通道不一致、频带不一致等。

3.2 阵列信号处理

首先，考虑N个远场的窄带信号入射到空间某阵列上，阵列天线由M个阵元组成，这里假设阵元数等于通道数，即各阵元接收到信号后经过各自的传输信道送到处理器，也就是说处理器接收来自M个通道的数据。

si(t)ui(t)ej(0t(t))

si(t)ui(t)ej(0(t)(t)) (3.2-1) 式中，ui(t)是接受信号的幅度，(t)是接收信号的相位，0是接收信号的频率。在窄带远场信号源的假设下，有

ui(t)ui(t) (3.2-2) (t)(t)根据式(3.2-1)和式(3.2-2)，显然有下式成立：

si(t)si(t)ej0 (3.2-3)

则可以得到第L个阵元接收信号为

xl(t)glisi(tli)nl(t) l1,2,,M (3.2-4)

i1N式中，gli为第L个阵元对第i个信号的增益，nl(t)表示第L个阵元在t时刻的噪声，li表示第i个信号到达第L个阵元时相对参考阵元的时延。

将M个阵元在特定时刻接收的信号排列成一个列矢量，可得

x1(t)g11ej011x(t)j021ge221j0M1x(t)MgM1eg12ej012g22ej022gM2ej0M2g1Nej01Ns1(t)n1(t)g2Nej02Ns2(t)n2(t)(3.2-5) s(t)n(t)gMNej0MNNM在理想情况下，假设阵列中各阵元是各向同性的且不存在通道不一致、互耦

等因素的影响，则式(3.2-4)中的增益gli可以省略(即归一化1)，在此假设下式(3.2-5)可以简化为

x1(t)ej011x(t)j0212ej0M1x(t)Meej012ej022ej0M2ej01Nej02Ns1(t)n1(t)n(t)s(t)22 (3.2-6)

s(t)n(t)ej0MNNM将式(3.2-6)写成矢量形式如下：

x(t)As(t)n(t) (3.2-7)

式中，x(t)为阵列的M1维快拍数据矢量，n(t)为阵列的M1维噪声数据矢量，s(t)为空间信号的N1维矢量，A为空间阵列的MN维流型矩阵(导向矢量阵)，且

Aa1(0)a2(0)aN(0) (3.2-8)

其中导向矢量

exp(j01i)exp(j)02i i1,2,,N (3.2-9) ai(0)exp(j)0Mi式中02f2c，c为光速，为波长。

由上述的知识可知，一旦知道阵元间的延迟表达式τ，就很容易得出待定空间阵列的导向矢量或阵列流型。下面推导一下空间阵元间的延迟表达式。假设空

间任意两个阵元，其中一个为参考阵元(位于原点)，另一个阵元的坐标为(x ,y, z)，两阵元的几何关系见图，图中“×”表示阵元。

图2 空间任意两阵元的几何关系

由几何关系可以推导出两阵元的波程差为

1c这里的波程差其实就是位于x轴上两阵元间的延迟、位于y轴上两阵元间的延迟和位于z轴上两阵元间的延迟之和。

(xcoscosysincoszsin) (3.2-10)

根据式(3.2-10)的结论，下面给出实际环境中常用的几种阵列及阵元间的相互延迟表达式。

(1)平面阵设阵元的位置为(xk,yk)(k1,2,,M)，以原点为参考点，另假设信号入射参数为(i,i)(i1,2,,N)，分别表示方位角与俯仰角，其中方位角表示与x轴的夹角。

(2)线阵设阵元的位置为xk(k1,2,,M)，以原点为参考点，另假设信号入射参数为i(i1,2,,N)，表示方位角，其中方位角表示与y轴的夹角（即与线阵法线的夹角），则有

1ki(xksini) (3.2-11)

c(3)均匀圆阵设以均匀圆阵的圆心为参考点，则有

kicoscr2(k1) (3.2-12) icosiM其中方位角表示与x轴的夹角，r为圆半径。

3.3 多重信号分类算法

多重信号分类(MUSIC)算法是Schmidt RO等人在1979 年提出的。这一算法

的提出促进了特征结构分类算法的兴起和发展，该算法已成为空间谱估计理论体系中的标志性算法。此算法提出之前的有关算法都是针对阵列接收数据协方差矩阵进行直接处理，而MUSIC算法的基本思想则是将任意阵列输出数据的协方差矩阵进行特征分解，从而得到与信号分量相对应的信号子空间和与信号分量相正交的噪声子空间，然后利用这两个子空间的正交性来估计信号的参数(如入射方向、极化信息及信号强度等)。MUSIC算法在特定的条件下具有很高的分辨力、估计精度及稳定性。下面仅介绍经典MUSIC算法。

窄带远场信号的DOA 数学模型为

X(t)A()S(t)N(t) (3.3-1)

对于阵元间距2的均匀线阵，阵列的导向矢量

exp(ju) (3.3-2) a()exp(j(M1)u)显然，上式可以表示成usin函数，即1u1就对应线阵的观察范围为

9090。

1T阵列数据的协方差矩阵为

REXXH

HH2AESSAI  ARSAH2I (3.3-3) 由于信号与噪声相互独立，数据协方差矩阵可分解为与信号、噪声相关的两

部分，其中RS是信号的协方差矩阵，ARSAH是信号部分。

对R进行特征分解有

RUUH (3.3-4-a)

式中，U为特征矢量矩阵，其中由特征值组成的对角阵如下：

1上式中的特征值满足如下关系：

2 (3.3-4-b) M12定义如下两个对角矩阵：

NN1M2 (3.3-5)

1sN1N显然当空间噪声为白噪声时，有

2 (3.3-6-a) N (3.3-6-b) MN2N2I(MM)(MN) (3.3-6-c)

式中，S为大特征值组成的对角阵，N为小特征值组成的对角阵。再将特征矢量矩阵分为与特征值对应的两部分：一是与大特征值对应的信号子空间USe1e2eN；二是UNeN1eN2eM，即与小特征值对应的噪声子空间。这样，式（3.3-4-a）可以进一步写成如下形式：

ReeHiiii1NjN1eeMHjii

USHUNUSUN

H USSUSHUNNUN (3.3-7)

式中，US是由大特征值对应的特征矢量张成的子空间也即信号子空间，而UN是由小特征值对应的特征矢量张成的子空间也即噪声子空间。

理想条件下数据空间中的信号子空间与噪声子空间是相互正交的，即信号子空间中的导向矢量也与噪声子空间正交

aH()UN0 (3.3-8)

经典MUSIC算法正是基于上述这个性质提出的，但考虑到实际接收数据矩阵是有限长的，即数据协方差矩阵的最大似然估计为

L1ˆRXXH (3.3-9) Li1ˆ。由于噪声的对R进行特征分解可以计算得到噪声子空间特征矢量矩阵UNˆ并不能完全正交，存在a()与U也就是说式(3.3-8)并不成立。因此实际上求DOAN是以最小优化搜索实现的，即

ˆUˆHMUSICargminaH()UNNa() (3.3-10)

所以，MUSIC算法的谱估计公式为

PMUSIC1 (3.3-11) HHˆˆa()UNUNa()下面给出MUSIC算法的流程：

ˆ，即式(3-9)； (1)由阵列的接收数据得到数据协方差矩阵Rˆ进行特征分解； (2)对Rˆ得特征值进行信号源数判断； (3)由Rˆ与噪声子空间Uˆ； (4)确定信号子空间USN(5)根据信号参数范围由式(3-11)进行谱峰搜索；

(6)找出极大值点对应的角度就是信号入射方向。实际非理想情况下得到的协方差矩阵的特征值满足：

12NN1M (3.3-12)

而并不满足式(3.3-5)，所以判断信号源数需要用到有关信号源数估计的方法，如AIC准则、MDL准则，盖氏圆法，平滑秩序列法等。

四、结果及分析

4.1仿真源程序

clc; clear ; close all; tic

M=8;%阵列的天线数 N=3;%信源数

snap=100;%快拍数目，即时域的采样点数 L=snap; C=3e8; %光速 w0=2*pi*1100; d=20;%阵列间隔 lamda=2*d; %波长

theta1=5; theta2=30;

theta3=50;%三个信号入射方向角，是与线阵法线的夹角

snr=10;%信噪比值 fs=1000; %采样频率

ts=1/fs; %采样间隔

t=(0:snap-1)*ts;%采样时刻（是个向量）

a=[0:M-1]';%阵列矢量 %非相干信号源 s1=sin(w0*t);

%s2=sin(1.1*w0*t); %s3=sin(1.2*w0*t);

%阵列导向矢量，把计算延时的公式直接代入了

a_theta1=exp(-j*2*pi*d/lamda*a*sin(theta1/180*pi)); %a_theta2=exp(-j*2*pi*d/lamda*a*sin(theta2/180*pi)); %a_theta3=exp(-j*2*pi*d/lamda*a*sin(theta3/180*pi)); A=[a_theta1];

%A=[a_theta1 a_theta2 a_theta3];%子阵的导向矢量，每个子阵完全相同 %S=[s1;s2;s3]; S=[s1];

X0=A*S;%阵列接收到的不含噪声的信号 R_s=S*S';

%%%%%%%产生一个和X0阵列相同的噪声%%%%

randn('state',0);

real_noise=randn(size(X0)); randn('state',3);

imag_noise=randn(size(X0));

noise=(real_noise+j*imag_noise)/(2^0.5); noise=10^(-snr/20)*noise; %%%%%%%%%####传统的MUSIC算法 X=X0+noise;

R_x=X*X'/length(t);%计算协方差矩阵

[U,V]=eig(R_x); %计算协方差矩阵的特征向量和特征值 Un = U(:,1:7);

theta0=-80:0.1:80;

theta_all=(theta0/360)*2*pi; for i1=1:length(theta_all)

theta = theta_all(i1)

a_theta11=exp(-j*2*pi*d/lamda*a*sin(theta)); %a_theta21=exp(-j*2*pi*d/lamda*a*sin(theta));

%a_theta31=exp(-j*2*pi*d/lamda*a*sin(theta)); A1=[a_theta11];%阵列导向矢量

%A2=[a_theta11 a_theta21 a_theta31 ]; P1(i1)=1/abs(A1'*Un*Un'*A1); end

plot(theta_all/pi*180,P1);

4.2实验结果及分析

(1)信噪比snr=10db;采样率 fs=1000;