基于Chirp-z变换的快速谱峰搜索方法

来源：宝玛科技网

设计与实现

基于Chirp-z变换的快速谱峰

【摘要】为减少空间谱谱峰搜索过程的计算量，文章提出了一种先采用FFT算法进行预搜索，确定峰值的大概位置，然后再采用Chirp-z变换算法在峰值附近进行精细搜索以确定峰值精确位置的方法。经仿真试验，该方法能够有效提高搜索效率，减少计算量。

【关键词】波达方向估计空间谱 Chirp-z变换谱峰搜索

1 引言

信号的波达方向（Direction of Arrival，DOA）估计是信号处理中一个重要研究内容，针对这一问题的算法大多需要对波达方向进行搜索，例如常规波束形成算法；高分辨率谱估计算法，如MUSIC算法[1]、Capon算法、线性预测算法等，由于其高分辨率的特点，这一类算法需要在所有角度范围内进行精细的谱峰搜索，其庞大的运算量成为工程实现的瓶颈。围绕这个问题，大量学者对谱峰搜索的过程进行优化研究。Liang Tao 等[2]和崔明岳等[3]在波达方向搜索过程采用FFT算法，简化了运算，但当阵元数较少时，由于其基于常规波束的算法，分辨率会有限。P. Stoica[4]提出，先采用阵列的较小子阵进行谱峰位置的预估计，再在此基础上采用整个阵列在谱峰位置附近做更精细的搜索，但其子阵的预估计过程有可能丢失峰值。本文将提出适合在一个小的区域内做精细搜索的方法，并适合与上述方法结合。

2 信号模型

考虑一个由N个全向阵元组成的线性阵列。假设M个远场窄带信号（M>N）分别从θi(i=1,2,…M)方向入射到接收阵列，且入射波信号和噪声不相关，以阵元1为参考阵元，各阵元与参考阵元间距为di(i=0,1,…N-1,d0=0)，如图1所示，则第k个阵元的输出信号可表示为：

（1）

其中si(t)为第i个信号的复包络，λ为中心波长，nk(t)为第k个阵元中的零均值高斯加性白噪声。

收稿日期：2012-03-22

图1 均匀线性阵列

2012年第14期

责任编辑：陈雍君********************

则阵列的输出信号矢量可表示为：

X(t)=[x1(t),x2(t),…xN(t)]T=A(θ)S(t)+N(t) （2）式中：

A(θ)=[a(θ1),a(θ2),…，a(θM)] （3）

令Zk=AW-k，k=0,1,…K-1。K为所要分析的复频谱点数，而A=ejθ0，为单位圆上复频谱分析的起点位置，W=e-jφ，φ为复频谱分析角度步长，则：

zk=AW−k=ejθ0ejkϕ=ej(θ0+kϕ) （9）

因此式（8）可表示为：

为N×M维阵列流形矩阵，a(θi)为对应的方向向量，且有

Y(zk)=∑y(n)A−nWkn k=0,1,…K-1 （10）

n=0

N−1

ak(θi)=e

-j

2πλdisinθi

,k=1,2,LN （4）

式（10）实际上是区间[θ0,θ0 +(K-1)φ]上的DFT变换，直接计算需要NK次复数乘法和（N-1）K次复数加法。采用Bluestein[11]提出的等式：

S(t)、N(t)分别为入射信号矢量和噪声矢量。常规波束形成算法对阵列输出式（2）直接进行搜索估计DOA，而高分辨率谱估计技术需要进一步处理，以MUSIC算法为例，由式（2）阵列输出信号自相关阵为：

nk=[n2+k2−(k−n)2] （11）

将式（11）代入式（10）：

Y(zk)=∑y(n)AWW−n

n=0N−1

n2−(k−n)2R=E[X(t)X(t)]=A(θ)RsA(θ)+σI （5）式中E为求均值运算，Rs=E[S(t)SH(t)]，当信号不相干时，Rs为满秩的M维对角阵。对R进行特征分解，由R的N-M个小特征值对应的特征向量U张成N-M维噪声子空间，采用MUSIC算法来估计信号DOA时，空间谱为：

P(θ)=

= （6）H

a(θ)UUa(θ)aH(θ)U2

-j2πHH2

Wk2=W

kN−12n=0

∑[y(n)A

−n

W]W

n2−

(k−n)2

（12）

令：

g(n)=y(n)AW

−n

n22[1]

，n=0,1,…N-1 （13）

h(n)=W

则：

−

n22

（14）

aH(θ)=[1,e

λd2sinθLe

-j

2πλdN−1sinθ](θ=−

π2

,−

π+Δθ,L,)22

π （7）

Δθ为搜索步长。观察式（6）和式（7），实际上当阵列为均匀阵列时，运算aH(θ)U可以看作对U的每一列做的频域非均匀DFT变换（因为对均匀步长Δθ做了非线性变换），因此如果采用均匀DFT变换，则实际上做的是非均匀角度步长的搜索，这也是搜索过程可以用FFT算法来实现的原因，引言部分介绍的方法基本上均基于这一点。采用FFT算法搜索是对整个角度域进行搜索，且其步长是固定的，为进行精细的搜索只有对U的每一列进行补零操作[6]。下面介绍采用Chirp-z变换仅对某一区域进行精细搜索的方法。

g(n)h(k−n)=W[g(k)∗k(k)],k=0,1LK−1 Y(zk)=W∑n=0

（15）

kN−12k2

∗表示卷积运算。式（15）表明Y(Z)的计算可以通

过卷积运算来实现，而卷积运算可以通过基于FFT的算法来快速实现，这便是Chirp-z变换的基本原理，其运算流程可以通过图2表明。

3 基于Chirp-z变换的搜索方法

设y(n)(0≤n≤N-1)为有限长序列，其Z变换为：

图2 Chirp-z变换算法流程

由式（14）知线性系统h(n)是非因果的，h(n)的取值在n=-N+1到n=K-1之间，为有限长序列，点数为N+K-1。

Y(z)=∑y(n)z

n=0

N−1

−n

（8）

输入g(n)点数为N，两者卷积输出长度为2N+K-2。由于序

责任编辑：陈雍君********************

2012年第14期

列DFT乘积对应的反变换是序列的圆周卷积，因此圆周卷积代替线性卷积不产生混叠的点数为2N+K-2，但由于只需Y(Zk)的前K个值，因此圆周卷积点数缩减到N+K-1。综上列出Chirp-z变换的实现步骤为：

（1）选择一个最小的整数L，使得L≥N+K-1，且L=2m，以便采用基2FFT算法。

（2）对g(n)进行补零使其变为L点序列。采用FFT算法计算此序列的L点DFT为G(r)。

（3）对h(n)进行周期延拓形成L点序列h(n)：⎧−n

⎪W2,0≤n≤K−1⎪ （16）h(n)=⎨0,K≤n≤L−N

⎪(L−n)2

⎪W−2,L−N+1≤n≥L−1⎩

法必需计算整个单位圆上的DFT，因此要对K’乘以2），且L’=2m，则采用FFT算法的复数乘法运算为(L’log2L’)/2。而采用Chirp-z变换来搜索，可以不改变搜索步长而根据实际需要来设计变换点数。

设阵元个数为37，阵元间距d=λ/2，λ为信号波长。四个等功率窄带信号分别从-35°、-33°、12°和32°四个方向入射到阵列。通过点FFT预估计得到的功率谱如图3所示：

采用FFT算法计算此序列的L点DFT为H(r)。（4）将G(r)与H(r)相乘，得Q(r)= G(r)H(r)。（5）采用FFT算法求Q(r)的L点IDFT为q(k)，实际上只需计算q(k)的前K个值。

（6）最后求Y(zk) =W

图3 基于FFT预估计和Chirp-z变换的功率谱搜索结果确立3个峰值位置的大概位置，然后在3个峰值附近分别采用Chirp-z变换进行精细搜索，设搜索步长为0.5°（注

q(k),0≤k≤K−1 。

意此时对应的是单位圆上0.5°），变换的点数为K=40，则L=，此时K’=360，L’=1024。图3给出了采用Chirp-z变换进行精细搜索得到的功率谱。表1给出了DFT、FFT和Chirp-z变换三种搜索方法的复数乘法量对比。

运用Chirp-z变换进行式（6）的搜索过程时，用U的每一列替代y(n)。由于Chirp-z变换适合对某一区域进行精细的搜索，因此可以先采用FFT算法对式（6）进行预搜索，确定峰值的大概位置，当阵元数较小时可以采用补零操作。然后采用Chirp-z变换算法在峰值附近进行精细的搜索以确定峰值的精确位置，同时分辨出靠得近的峰值。值得注意的是，对式（6）进行搜索时，只要搜索的步长足够小便能分辨出峰值，这与FFT算法的补零操作不提高分辨率是有区别的。

表1 DFT、FFT和Chirp-z变换搜索方法的复数乘法量

搜索方法复数乘法（次）

DFT13320

FFT5120

Chirp-z变换

1025

表1中，Chirp-z变换的复数乘法量包括了预估计带来的计算量，可见其计算量上的优势。当N很大而搜索步长很小时，其计算量上的优势会更明显。Chirp-z变换适合对某一区域进行精细的搜索，适合于DOA的跟踪，因为一

4 性能分析及仿真试验

由Chirp-z变换的实现过程，其所需复数乘法运算为

旦确定DOA的初始值，则其跟踪的搜索过程可以局限在一个小的附近区域[10]。

Llog2L+5N+L+K （17）2

由式（9）知，对K点Chirp-z变换，其计算的DFT在单位圆上的起点和步长可以任意控制，并且对阵元个数N和变换点数K没有。因此Chirp-z变换为搜索过程提供了很大的灵活性。设搜索步长为∆θ，采用DFT计算的点数为K’=π/∆θ，复数乘法运算为NK’，设L’≥2K’（由于FFT算

责任编辑：陈雍君********************

5 结论

综上所述，本文提出的基于Chirp-z变换的搜索方法适合对一定区域进行精细的搜索。通过结合FFT算法先进行预搜索，确定峰值的大概位置，然后采用Chirp-z变换算法在峰值附近进行精细搜索，能够提高搜索效率，减少计算量。而且Chirp-z变换计算的DFT在单位圆上的起点和步

2012年第14期

长可以任意控制，可以根据实际需要以任意搜索步长来设计变换点数。由于对阵元个数NK没有，因此Chirp-z变换为搜索过程提供了很大的灵活性，适合于DOA的跟踪，恰当的运用可以大大减少搜索过程计算量。另外，本文提出的搜索方法并不只限于MUSIC算法的搜索过程，事实上在需要进行谱峰搜索的场合均可适用。

1999,6(10): 273-275.

[5] A Jakobsson. Computationally efficient two-dimensional

Capon spectrum analysis[J]. IEEE Trans Sig Proc,

2000,48(9): 2651-2661.

[6] A S Paine,Qinetiq Malvern. Fast MUSIC for large 2-D

element digitised phased array radar[C]. Proceedings of the international, 2003. ★

参考文献：

[1] R O Schmidt. Multiple emitter location and signal

parameter estimation[C]. Proc. RADC Spectral Estimation Workshop, Rome, NY, 1979: 243-258.

[2] Liang Tao,H K Kwan. A novel approach to fast DOA

estimation of multiple spatial narrowband signals[C]. The

作者简介

王依恺：本科，主要从事通信技术研究与总体规划工作，重点研究方向为传输系统设计、IP网络技术的实时应用和网络时统技术应用。

45th Midwest Symposium on Circuits and Systems, 2002: I-411-I-414.

[3] 崔明岳,刘云安. 信号到达角的快速估计算法[J]. 北京邮

电大学学报, 2003,26(1): 69-72.

[4] P. Stoica. Maximum-Likelihood DOA Estimation by Data-Supported Grid Search[J]. IEEE Signal Processing Letters,

英运营商开通免费Wi-Fi区借助奥运搜索新商业模式据悉，在今年夏季伦敦奥运会开幕式即将举行之际，英国移动运营商O2已在伦敦西区（West End）开通一个免费Wi-Fi互联网接入服务区。

该免费Wi-Fi服务覆盖范围包括牛津街（Oxford Street）和摄政街（Regent Street）等购物区、莱斯特广场（Leicester Square）周围的剧院区以及博览会路（Exhibition Road）周边的各个博物馆。在本届奥运会结束后，该Wi-Fi服务区仍将得以保留。

该Wi-Fi服务由O2同西敏斯特（Westminster）及肯辛顿切尔西（Kensington & Chelsea）这两个自治市联合推出。该服务本身属免费性质，但用户需观看相应广告内容。在经过一次性登录后，用户将首先看到一个O2的主页，然后就是英国零售商House of Fraser发布的商品推销内容。

O2在一则新闻稿中表示，今后该Wi-Fi服务所发布广告内容，将有可能根据“用户所在不同位置”，而投放针对性更强的广告，从而使用户能够快速找到当前所处位置周围的商家。

西敏斯特市政委员会负责人菲利帕•罗伊（Philippa Roe）表示，此次开通的免费Wi-Fi服务，其实也是伦敦的一个长期性项目，目的是将伦敦打造成全球范围内“科技人性化程度最强的城市之一”。

从今年5月开始，英国第二大付费电视服务商维珍传媒（Virgin Media）也在部分伦敦地铁站中开通了类似服务，并计划在伦敦奥运会旅游旺季期间对该服务采取收费措施。目前尚未清楚O2的免费Wi-Fi服务究竟将持续多长时间。但O2一位发言人表示，今年奥运会结束后，该免费服务仍将持续“很长一段时间”。（腾讯科技）

责任编辑：陈雍君********************

2012年第14期

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文

全部频道

基于Chirp-z变换的快速谱峰搜索方法