基于MATLAB的阵列信号处理仿真方法

基于MATLAB的阵列信号处理仿真方法

统 仿 真 学 报© Vol. 20 No. 13

2008年7月 Journal of System Simulation Jul., 2008

第20卷第13期 系

基于MATLAB的阵列信号处理仿真方法

刘 玲，刘晓明，曾 浩

（重庆大学 通信工程学院，重庆 400030）

摘 要：介绍如何使用MATLAB构建阵列信号处理系统模型，包括相干信号模型，幅度和相位误差模型，针对不同模型，实现协方差矩阵产生方法，波达方向估计的子空间方法，自适应波束合成器的权值求解算法和方向图、阵列增益等系统参数的仿真。这些仿真模型和方法，对于各种复杂的阵列信号处理研究，具有重要的基础作用。 关键词：阵列信号；协方差矩阵；波达方向；波束合成

中图分类号：TP391.9 文献标识码：A 文章编号：1004-731X (2008) 13-3548-05

Array Signal Processing Simulation Based on MATLAB

LIU Ling, LIU Xiao-ming, ZENG Hao

（Department of Communication Engineering, Chongqing University, Chongqing 400030, China）

Abstract: MATLAB is used to develop the system model of array signal processing, such as the correlated signals, the error

of amplitude and phase. To deal with these different models, MATLAB offers great instructions and fuctions, which make the simulating considerably convenient. The key programs were proposed, which could simulate the covariance matrix estimation, estimate the DOA (direction of arrival) of signal based on subspace, calculate the weight vector of adaptive beamformer, and get some basic system parameters. Some examples show the detail procedure and the programs are useful to the array signal processing research. Furthermore, these methods take an important function to many other complicated simulations.

Key words: array signal; covariance matrix; DOA; beamforming

引 言

阵列信号处理源于60年代相控阵天线技术[1]，目前，广泛用于雷达、通信、超声波、语音、水下探测等等不同领域。阵列信号处理的对象，是由阵列天线接收的采样快拍数据，而利用各种信号处理的算法和工具，主要实现两个目的

[2-4]

则第i个阵元接收信号相对于原点信号的时间延迟为

τi=

aTpi

(1) c

x

：一是进行空间谱估计，通过对信号分析，确定信号波

y

达方向；二是进行波束合成，利用盲和非盲算法，得到阵列加权矢量。就实现手段而言，早期相控阵中，采用模拟方式，而目前则是在数字域实现，利用FPGA或者DSP这些硬件平台，完成信号处理算法。阵列信号处理的算法研究，往往通过仿真进行。在MATLAB中，如何建立仿真模型，如何实现各种基本研究参数评价指标仿真，对于复杂的阵列信号处理具有重要意义。

图1 阵列几何结构示意图

其中，c是光速，而

sinθcosφ

(2) a= sinθsinφ

cosθ

1 阵列信号处理的信号模型

由于阵列信号处理对象都是阵列天线接收信号，所以，不同应用领域的信号模型是相同或相似的。对于一个远场窄带零均值的入射信号，其N阵元构成阵列如图1所示。如果用单位方向矢量a表示信号来向，矢量pi表示阵元坐标，

xi

pi= yi zi

i=1,2...,N (3)

如果原点位置接收信号的复数表示为

0(t)=x0(t)ejωct (4) x

则由于信号是窄带的，可以不考虑包络延时，阵列接收信号矢量可以表示为

e

jτ1 x1(t)

x(t)= ... =x0(t)ejωct ... (5)

e jτN xN(t)

3548

收稿日期：2007-08-20 修回日期：2008-02-21

基金项目：国家自然科学基金项目（30570473）

作者简介：刘玲 (1978- ) 女，四川广安人，博士生，研究方向为阵列信号处理，第四代移动通信技术；刘晓明（1963-）男，重庆长寿人，教授，研究方向为现代信号处理，无线通信技术。

基于MATLAB的阵列信号处理仿真方法

2008年7月 刘

玲，等：基于MATLAB的阵列信号处理仿真方法 Jul., 2008

通常采用前向平滑的估利用有限的K个快拍进行平滑估计。

计算法，这也就是最大似然法。

K

x=1∑x(k)xH(k)=1XXH (10) R

Kk=1K其中，数据矩阵

X= x(1)...x(K) =Xs+∑Xij+Xn

j=1J

定义波数矢量

sinθcosφ 2π (6) k= sinθsinφ λ cosθ 由此可以表示方向矢量

e v(k)= ...

e jkp

T

jkTp1

N

(7)

(11) ...i0j(K) +Xn

考虑阵列信号处理一般是在基带进行，信号中已经没有载波分量，接收信号可以表示为

=v(ks) s0(1)...s0(K) +

∑v(k) i(1)

ij

0j

j=1

J

x(t)=x0(t)v(k) (8) 由此可见，阵列接收信号的确定，一方面决定于信号的复基带表达式，另一方面则决定于方向矢量，而方向矢量又是由信号DOA和阵列几何结构确定的。

阵列接收的信号由于是由三个部分组成，数据矩阵也可以分解为三个部分。不论期望信号，干扰，还是热噪声，都假设是高斯分布的复数信号，则可以用randn（）函数模拟产生这些信号。

同样采用3.1中的假设条件，只是此时使用前向平滑进行协方差矩阵估计，MATLAB程序为

s0 = sqrt(SNR/2)*randn(1,K)+j*sqrt(SNR/2)*rand(1,K); i0 = sqrt(INR/2)*randn(1,K)+j*sqrt(INR/2)*rand(1,K); Xn = sqrt(1/2)*randn(N,K)+j*sqrt(1/2)*rand(N,K); Xs = Vs*s0; Xi = Vi*i0; X = Xs + Xi + Xn; Rx=X*X’/K;

2 阵列协方差矩阵

阵列接收信号包括三个部分，一部分是期望信号，一部分是干扰信号，一部分是热噪声。期望信号只有一个，干扰可以多个，而这些信号均值都是零。接收信号的自相关矩阵，是描绘不同信号成分间相互关系的主要统计量，而期望、干扰、噪声均值都为0，自相关矩阵就是协方差矩阵。算法仿真时，如和获得协方差矩阵统计量，是关键步骤。

2.1 协方差矩阵的公式法

在均值为零，信号和干扰相互独立前提下，阵列的接收信号协方差矩阵也就是自相关矩阵，其表达式为

Rx=E{x(t)xH(t)}

22

=σs2v(ks)vH(ks)+∑σijv(kj)vH(kj)+σnI

j=1J

2.3 相干信号的协方差矩阵估计

在阵列信号处理中，通常假设各个信号之间是不相干的。在这样的假设前提下，信号的协方差矩阵是满秩的，不会出现“病态”，协方差矩阵的逆阵也是存在的，从而，各种

(9)

DOA估计算法和波束合成算法才能正常使用。然而，由于多径或者人为干扰等情况的存在，信号间可能存在相干性。此时两个信号可以表示为

s1(t)=s0(t) αejβ (12)

=Rs+Ri+Rn

2

s

2ij

2n

这个协方差矩阵表达式是一个没有估计误差的表达式，

σ，σ，σ分别是期望信号，第j个干扰信号和热噪声

的方差，空间具有一个期望和J个干扰。MATLAB中的协方差矩阵确定，有两种方法，一种是直接使用（9）式，得到无估计误差的协方差矩阵。

比如，对于一个N阵元的标准线阵，假设空间期望信号us=cosθs，信噪比SNR=pdB,存在干扰信号

参数α和β反应了二者数值关系。为了描绘相干信号时协方差矩阵，则用randn（）产生第一个信号数据后，根据参数α和β，就可以得到相干信号数据。

s0 = sqrt(SNR/2)*randn(1,K)+j*sqrt(SNR/2)*rand(1,K); s1= s0*alpha*exp(beta);

通过2.2介绍的协方差矩阵估计模型，得到了相干情况下的Rx，但必须经过对协方差矩阵的去相干处理，获取满秩的协方差矩阵，才能用于后续阵列信号处理。一种经典的去相干处理方法是采样前后向平滑的协方差矩阵处理算法[5]。定义矩阵J，这是一个次对角线元素为1，其余为0的方阵，即

001

(13) J= 0...0

100

ui=cosθi，干扰噪声比为INR=qdB，对于热噪声，通常

假设σ=1。则MATLAB中采用（9）式可以得到协方差矩

2

n

阵，其程序为

n=[-(N-1)/2 : (N-1)/2]’; Vs=exp(-j*pi*us*n); Vi= exp(-j*pi*ui*n); SNR=10^(p/10); INR=10^(q/10);

Rx=SNR*Vs*Vs’+INR*Vi*Vi’+eye(N);

2.2 协方差矩阵的估计法

在实际的系统中，协方差矩阵只可能通过估计得到，即

3549

则协方差矩阵估计值为

基于MATLAB的阵列信号处理仿真方法

2008年7月 系 统 仿

K

=1∑Rx(k)xH(k)+Jx*(k)xT(k)Jx

2Kk=11

(14) =XXH+JX XTJ

2K

真 学 报 Jul., 2008

前提是信号个数J要是已知量。

对于一个均匀园阵，假设阵元N=8，阵列位于x-y平面，圆心位于原点，半径r为半波长，第i个阵元坐标可以表示为

2π

rcosN(i 1)

2π (20) pi=rsin(i 1) N

0

()

()

由上式可知，去相干协方差矩阵求解程序是在得到相干信号和热噪声产生的数据矩阵X后，采用

Rx=(X*X’+J*conj(X)* transpose (X)*J)/(2*K);

3 DOA估计的子空间算法

对于信号DOA估计，是阵列信号处理的一个重要研究内容。基于子空间的DOA估计和最大似然DOA估计是两种重要算法，但研究表明，这些算法却可以统一为加权的子空间拟和算法[6]。最大似然算法较子空间算法具有更好的

同时结合（6）式，可以得到方向矢量。

2π

jπsinθcos φ 1 1) N e ...v(k)= (21)

2π jπsinθcos φ N 1) N

e

C-R下界，可以克服相干信号影响，但子空间算法却在运算量上优势明显。

3.1 子空间模型

假设空间具有J个信号，他们的方差分别是σ，其中

2

j

假设空间具有三个信号，其DOA来向分别为（2600，200）、、（800，600），信噪比均为10dB，则MATLAB（1800，400）

程序首先应该得到协方差矩阵Rx，然后进行特征分解，最后是谱峰搜索。主要程序为

fi=2*pi/N*[0:1:N-1]’;

[F,D]=eigs(Rx,N-3,'SM'); xa=[0:2:360]*pi/180; ya=[0:1:90]*pi/180; for k=1:length(ya)

for m=1:length(xa)

Va=exp(j* pi*sin(ya(k))*cos(xa(m)-fi)); P(k,m)=1/(Va'*F*F'*Va); end end

mesh(x,y,10*log10(P));

j=1,..,J，而热噪声的方差为σn2，那么，协方差矩阵Rx的

N个特征值分别为

2

i=1,...,J σ2j+σn

λi= (15) 2

i=J+1,...,Nσn

对于J个较大特征值，其对应特征矢量为J个正交矢量，他们张成的空间成为信号子空间。

span{φ1...φJ} (16)

而对于N J个较小特征值，其对应特征矢量也是N J个正交矢量，而他们张成的空间就是噪声子空间。

span{φJ+1...φN} (17)

3.2 子空间MUSIC法

MUSIC算法是典型的子空间算法[7]，利用信号子空间和噪声子空间正交性，进行谱峰搜索。由于接收信号是位于信号子空间的，而信号子空间又是和噪声子空间正交的，所以，通过把接收信号方向矢量向噪声子空间投影，其值应该为零，即空间谱函数P(k)具有最大值。

P(k)=

1

(18)

v(k)Πv(k)

φJ+1

...φN] ... (19)

φN

相应仿真结果如图2所示。图的上部分为三维空间谱，下部分是空间谱的投影图。

其中投影算子

Π=[φJ+1

从空间谱表达式（18）中可以看出，进行谱峰搜索的主要问题在于方向矢量表示和协方差矩阵的特征分解。根据方向矢量的定义（7）式，其决定于信号波数和阵元坐标。阵元坐标是固定的，而谱峰搜索过程就是对波数中的DOA进行扫描过程，DOA是空间谱函数的自变量。对于噪声空间投影，虽然定义非常清楚，但在工程实现过程中，计算量却很大。对于MATLAB而言，特征值和对应特征矢量的求解相对较为简单，只需调用函数eigs（）就可以直接得到，但

图2 二维MUSIC谱

4 自适应波束合成技术

自适应波束合成技术，本质上就是空域自适应滤波。虽

3550

基于MATLAB的阵列信号处理仿真方法

2008年7月 刘 玲，等：基于MATLAB的阵列信号处理仿真方法 Jul., 2008

然目前研究前沿在于空域、时域、频域、码域的相互融合，但就阵列信号处理本身而言，自适应波束合成是重要研究内容。

w=inv(Rx)*Vs/(Vs'*inv(Rx)*Vs); u = -1:1/1000:1;

beam = 20*log10(w'*exp(j*n*pi*u)); plot(u,beam);

4.1 最优权值求解

得到空域自适应滤波的最优权矢量，是自适应波束合成的最终目标。通过对接收信号进行加权求和，可以使阵列主瓣对准期望，而零陷对准干扰，从而大大提高系统增益。最优权矢量解，对于不同的准则，有不同的表达形式。比如，对于最小功率准则的波束合成器，其最优权矢量表达式为

1

x

仿真结果如图3所示。图3中，实线表示采用公式法得到协方差矩阵，进而求解方向图。此时，主瓣方向和零陷都正常。而虚线表示在不同快拍数估计协方差矩阵时的方向图。图中可以看出，当快拍数较少，协方差矩阵误差较大，一是旁瓣电平较高，二是零陷深度较低，三是可以对期望信号也产生零陷，这就是“自零陷”效应。

[8]

wMP=µRv(ks) (22)

1

µ=H (23)

v(ks)Rx 1v(ks)

可见，权矢量决定于期望信号方向矢量和阵列接收信号的协方差矩阵。系数µ是常数，其意义在于保证主瓣期望信号方向上的增益为单位增益。

而求解最优权矢量解的过程，则是采用各种自适应的算法，无论是盲算法还是非盲算法。一种基于QR分解的逆

RLS算法可以归纳为，初始化信号协方差矩阵逆阵P和权矢量w为0，然后按如下步骤进行迭代：

第一步，构造前置矩阵

1λ-xH(n)P(n 1)

F(n)= (24) -110P(n1) λ 第二步，对上述矩阵进行QR分解，得到下三角后置矩阵

γ 1(n)0T

B(n)= (25) nγnnkP()()()

图3 阵列方向图

4.3 阵列增益

阵列增益是研究阵列性能的重要参数，其定义为阵列输出信干噪比和输入信干噪比之比[9]。

G=

SINRout

(29) SINRin

第三步，根据后置矩阵，计算第n次迭代预测误差ξ(n)和权矢量w(n)

式中，根据假设条件

SINRin=

ξ(n)=d(n) w(n 1)x(n) (26) w(n)=w(n 1)+k(n)ξ*(n) (27)

式中，d(n)为参考信号。在使用MATLAB仿真过程中，关键步骤在于由数据矢量x(n)构造前置矩阵后，无法直接进行QR分解得到后置矩阵，而是通过如下语句实现

H

σs2

∑σ

j=1

(30)

2n

2ij+σ

SINRout=

wHRsw (31) wHRi+Rnw

可见，阵列增益表示中，需要确定各个量的方差，信号协方差矩阵，以及权矢量，而这些量通过前文分析，都有各自实现方法。最小功率准则阵列增益程序为

Rs=SNR*Vs*Vs’; Ri=INR*Vi*Vi’; Rn=eye(N); Rx=Rs+Ri+Rn;

w=inv(Rx)*Vs/(Vs'*inv(Rx)*Vs); SINRin=SNR/(INR＋1)；

SINRout=w’*Rs*w/(w’*(Ri+Rn))*w); G=10*log10(SINRin/SINRout);

B=triu(qr(F’))’;

4.2 阵列方向图

通过阵列方向图，可以直观了解阵列波束是否具有合理指向。方向图是阵因子和阵元因子乘积，通常自适应波束合成信号处理算法中不考虑阵元因子，故阵列方向图阵因子为

H

Beam=wv(k) (28)

式中，权是通过自适应迭代获得，而方向矢量中的角度是自变量。但在不考虑自适应迭代算法自身情况下，若要对阵列参数进行仿真，可以直接采用式（22）权矢量。

假设空间均匀线阵N=8，空间期望信号为us=0，两个干扰信号u1=0.4，u2= 0.6，INR1=10dB，SNR=20db，

INR2=20dB则方向图的仿真首先是求解协方差矩阵，然后

同样假设一个均匀线阵，空间期望信号us=0，

SNR=20dB，同时一个干扰信号，ui=0.4，考察不同INR

下阵列增益同阵元数关系，如图4所示。图中可以看出，阵元数越多，干扰功率越小，增益越大。

是权值表达式，最后用式（28）描绘出方向图。

3551

基于MATLAB的阵列信号处理仿真方法

2008年7月 系 统 仿

真 学 报 Jul., 2008

transactions on antennas and propagation (S0096-1973), 1964, 26(5): 161-169. [2]

Yuan Q, Chen Q, Sawaya K. Accurate DOA estimation using array antenna with arbitrary geometry [J]. IEEE transactions on antennas and propagation (S0096-1973), 2005, 53(4): 1352-1357. [3]

Abdel-Samad A, Davidson T N, Gershman A B. Robust Transmit Eigen Beamforming Based on Imperfect Channel State Information [J]. IEEE transactions on signal processing (S1053-587X), 2006, 54(5): 1596-1609. [4]

Krim H, Viberg M. Two decades of array signal processing research [J]. IEEE siganl processing magazine (S1053-5888), 1996, 13(4): 67-94.

图4 阵元数同阵列增益

[5] [6]

曾浩，刘玲，覃剑. 时间平滑ML协方差矩阵估计算法及性能分析[J]. 系统仿真学报，2007，19(19)：4517-4520.

Viber M, Ottersten B, Kailath T. Detection and estimation in sensor arrays using weighted subspace fitting [J]. IEEE transactions on signal processing (S1053-587X), 1991, 39(11): 2436-2449.

5 结论

阵列信号处理内容非常丰富，各种算法不断发展和改进。但基于MATLAB的阵列信号处理仿真模型和基本方法，确始终没有改变。基本内容仍是上文所述协方差矩阵获取，谱函数搜索，权值求解，方向图和阵列增益描绘。所以，这些仿真方法对于阵列信号处理研究，有重要意义。

[7] Schmid R O. Multiple emitter location and signal parameter

estimation [J]. IEEE transactions on antennas and propagation (S0096-1973), 1986, AP-34(3): 276-280. [8]

Van Trees H L. Optimum array peocessing [M]. USA: John Wiley & Son, 2002: 451-452.

[9] 张贤达，保铮. 通信信号处理[M]. 北京：国防工业出版社，2002：

320-325.

参考文献：

[1] Chose R N. Electronically adaptive antenna systems [J]. IEEE

(上接第3532页)

能为零，则表明系统没有完成均衡。

综上所述，两余度串联舵机的伺服回路采取如图8所示的均衡线路时，当两通道的差异在一定范围内时，可通过“增弱减强”的原理对输入信号进行均衡，若不一致量在一定时间内持续超出规定的均衡权限，则报出均衡告警信号。

数学模型，研究了各种情况下串联舵机伺服系统的时域和频域特性，并在此基础上对双余度伺服回路均衡控制进行了仿真分析。结果表明，当输入Ui在线性段范围内变化时，串阻尼0.45的近似二阶系联伺服系统模型可简化为频率8Hz，

统。均衡线路通过“增弱补强”的原理，使双余度伺服通道协调输出，保证了系统工作的高精度和稳定性。该伺服控制方案可应用于同类直升机串联舵机伺服系统设计中，且等效得到的简化模型可作为伺服回路的监控模型用于飞行中实时检测。

参考文献：

[1] 姚纪文. 自动控制元件及其线路[M]. 北京: 国防工业出版社,

1980.

[2] 舒志兵, 刘峻泉, 林锦国, 等. 闭环伺服系统的数学模型研究[J].

系统仿真学报, 2002, 14(12): 1611-1613.

[3] 臧怀泉, 黄镇海, 尹汝泼, 等. 液压伺服系统建模的新方法[J]. 计

算机仿真, 2002, 19(5): 53-55.

[4] 杨云, 沈毅力, 李天石. 流体脉宽调制控制系统的建模与仿真[J].

系统仿真学报, 2003, 15(2): 171-178.

[5] 邹海峰, 孙力, 阎杰. 飞行器舵机电动伺服加载系统研究[J]. 系统

仿真学报, 2004, 16(4): 657-659.

[6] 段海滨, 于秀芬, 王道波. 基于内模PID鲁棒控制的飞行仿真伺服

系统设计[J]. 中国空间科学技术, 2004, (6): 1-6.

图10 J1变化320%时伺服系统的响应曲线

4 结论

本文建立了某型自动驾驶仪俯仰串联舵机伺服系统的

3552

本文来源：https://www.bwwdw.com/article/lfk1.html