多子阵波束域Root-MUSIC算法并行处理实现方法

刘千里

（海军驻武汉四三八厂军事代表室，武汉，430060）

多子阵波束域Root-MUSIC算法并行处理实现方法

刘千里

（海军驻武汉四三八厂军事代表室，武汉，430060）

在对MSB-RMU算法计算量分析的基础上，研究了基于并行Jacobi方法的特征值分解算法，并设计了一种基于FPGA的采用双边CORDIC结构的并行Jacobi实现方法，为特征分解类算法提供了一种高效的单芯片解决方案。最后讨论了MSB-RMU算法在一块由6片DSP和2片FPGA组成的并行信号处理平台的实现过程。试验结果表明，设计的并行处理方法能够满足多子阵波束域Root-MUSIC算法的实时实现。

多子阵波束域；Root-MUSIC；并行Jacobi；FPGA；DSP

对于海底地形地貌声呐探测系统来说，关键的要素之一就是算法的研究与实现。因此，除了算法本身的高性能之外，考察海底探测算法另一个重要的标准就是其是否具有可实现性，即算法是否能实时实现。

多子阵波束域CAATI算法（MSB-CAATI）和多子阵波束域Root-MUSIC算法（MSB-RMU），两者都能够同时估计回波的方向和强度，非常适用于海底地形地貌探测系统［1］。但后者具有更好的估计性能，且实现过程较前者简单［2］。本文在对MSB-RMU算法计算量分析的基础上，研究了基于并行Jacobi方法的特征值分解算法，并设计了一种基于FPGA的采用双边CORDIC结构的并行Jacobi实现方法，为特征分解类算法提供了一种高效的单芯片解决方案。最后讨论了MSB-RMU算法在一块由6片DSP和2片FPGA组成的并行信号处理平台的实现过程。

1 算法计算量分析

MSB-RMU算法中计算量主要集中在多个子阵的波束形成、子阵波束域输出的相关矩阵运算、相关矩阵的特征分解以及多项式求根得到DOA，再根据阵元域协方差矩阵估计回波强度就可以进行海底成像。MSB-RMU算法较传统Root-MUSIC算法增加了多个子阵的波束形成，但相关矩阵阶数减少，从而特征分解的运算量大大减小。由此可见，对于大阵元数系统而言，总的运算量要低于传统算法。即便如此，算法还是具有相当的运算量，采用单独的处理器仍难以实时实现，因此还必须采用大规模并行处理平台进行计算。

MSB-RMU算法将MUSIC算法的高分辨性能和FT波束形成算法的鲁棒性有机结合，同时算法将阵元域处理转换到多个子阵的波束域处理，从而大大减小了对平均协方差矩阵RM的特征值分解运算量。尽管如此，对于大阵列情况，当子阵个数较多时，平均协方差矩阵RM的阶数仍然较大，也就成为了整个算法中运算量相对集中的部分。因此高效地实现对平均协方差矩阵的特征值分解是MSB-RMU算法实时实现的关键。

2 并行Jacobi方法实现特征值分解（EVD）

2.1 Jacobi方法特征值分解

Jacobi方法简单紧凑、精度高、收敛较快，是计算实对称矩阵全部特征值和相应特征向量的有效方法［3］。Jacobi算法采用平面旋转矩阵进行正交相似变换，将实对称矩阵化为对角矩阵，用于求解实对称矩阵的全部特征值和特征向量。下面给出了基于Jacobi方法的特征分解基本形式。

其中，R为N×N实对称矩阵，Λ为对角矩阵，其对角线上的元素即为矩阵R的特征值，I为单位矩阵，U矩阵的列向量为所有特征值对应的特征向量。Wpq为平面旋转矩阵：

2.2 并行Jacobi方法特征值分解

并行Jacobi算法［4］的基本思路是每次旋转消去多于一对的非对角元素。若矩阵R为N×N实对称矩阵，令m=N/2，显然R有N（N-1）/2对非对线元素需要消去，每次正交变换可同时消去的最大数为N/2。把每2m-1变换记为一次扫描，并行Jacobi算法若能满足下面两个条件，并行计算的效率将达到最高。

（1）构造正交变换矩阵Qj，使每一个变换都能消去N/2对非对角元素。

（2）每次扫描，任一对非对角元素只被消去一次，即在一次扫描中，2m-1个正交变换的每一个变换消去的元素，都不同于其他变换消去的元素。

对于Qj的选取有很多种方法，文献［5］给出一种选取方法。对每一次扫描，其2m-1个正交变换中的每一个变换Qj的非零元素按以下方法选取。

其中（p，q）的按照下面定义：

并行Jacobi算法每次正交相似变换可同时消去N/2对非对角元素，因此与经典Jacobi算法相比，在Jacobi扫描次数相同的情况下，所用时间约为经典Jacobi算法的2/N。

3 基于并行Jacobi方法的EVD在FPGA上的实现

利用CORDIC算法来实现Jacobi算法已被广泛的研究［6］，而利用双边CORDIC在FPGA上实现并行Jacobi算法可以得到更高的效率。基于2.2的分析，可以把每次正交相似变换记为一次消去过程，每一个变换消去的元素，都不同于其他变换消去的元素，因此消去过程可同时进行。每个消去过程可由一个双边CORDIC来实现，因此需要N/2个双边CORDIC组成一个并行处理网络。双边CORDIC解决如式（12）所示的双边旋转问题。

单个双边CORDIC处理器原理框图如图1所示。

图1 双边CORDIC原理框图

图中nσ为加法器的符号控制位，根据具体情况来选择，针对Jacobi正交相似变换，符号的选取原则是使待消去的元素逐渐变为零。

特征向量的计算也采用CORDIC来实现，所不同的是每次旋转过程中的符号由消去过程所确定。由于并行Jacobi算法是用于计算实对称矩阵特征分解的方法，因此对于Hermite矩阵而言，还需经过下面的变换［7］。

将N阶Hermite矩阵H写成实数矩阵和复数矩阵两部分，即

可见，H矩阵的特征分解问题转化为实对称矩阵的特征分解问题。

由于CORDIC算法仅需移位和加法运算即可完成，使得硬件实现起来非常简单。XILINX公司的IP核生成软件Core Generator中自带CORDIC算法的IP核，该处理核的结构如图2所示。

图2 CORDIC核结构

该CORDIC核具有两种可选的结构配置选项：一种是完全并行配置模式，该模式具有单周期数据输入输出，但布局布线所占用的面积较大；另一种是字串行实现模式，该模式的数据输入输出需要多个时钟周期，但布局布线所占用的面积较小。算法需要产生一个尺度因子来对输出结果的幅度进行校正，CORDIC核提供自动补偿的选项来进行幅度校正。CORDIC算法过程中，每个单位精度需要近似一次移位-加减操作。若选择字串行结构配置，N的输出宽度需要N个时钟周期的延迟，并在每N个时钟周期后产生一个新的输出；若选择并行结构配置，输出宽度需要N个时钟周期的延迟，并在每一个时钟周期后产生一个新的输出。因此在设计时应注意时钟延迟的影响，并根据具体工程需要在速度和实现面积上进行选择。

4 多子阵波束域Root-MUSIC算法的并行处理系统设计

4.1 并行处理系统组成

为了满足系统运算量需求，设计了一块由6片DSP和2片FPGA组成的CPCI结构并行处理系统，系统框图如图3所示。

图3 并行处理平台结构框图

系统主要完成对多个通道原始数据的正交变换和基于MSB-RMU算法的DOA估计，其中正交变换、多个子阵波束形成以及特征分解都由FPGA来完成，自相关矩阵的计算和最后的多项式求根由DSP阵列完成。下面按照信号流的顺序分别讨论各主要信号处理部分的实现方法。

4.2 基于CORDIC的数字下变频器

数字下变频器（DDC）的作用是将模数变换器采集到的原始数据载波频率搬移到基带进行处理。传统的DDC由正余弦查找表、乘法器和低通滤波器组成，如图4（a）所示。也可以采用CORDIC算法加以实现，从而无需乘法器和存储正余弦表。CORDIC算法由Volder［8］最早提出用于计算平面直角坐标系和极坐标系之间的变换。其基本思想是用一系列固定角度的不断偏摆来逼近所需旋转的角度，适当选取一些固定的角度值，可以使CORDIC运算只需进行移位和加减操作，适合硬件实现。所以在本系统中采用FPGA来实现DDC。利用CORDIC算法实现的DDC原理框图如图4（b）所示。

图4 DDC实现框图

4.3 基于DSP阵列的协方差矩阵计算

完成多子阵波束形成后，将所有子阵大于平均波束能量的输出波束送给4个DSP组成的DSP阵列，每个DSP从子阵中挑选出符合要求的输出波束，由13个子阵的对应方向输出组成新的信号矢量，其协方差矩阵表示为Ri，i=1，2…，M，再将M个协方差矩阵平均，得到平均协方差矩阵RM，对RM进行特征分解即可得到多子阵波束域噪声子空间和信号子空间。

由于矩阵运算的相对独立性，所以协方差矩阵可以采用并行计算的方式来实现。对于本系统而言，协方差矩阵Ri，为13×l3的Hermite阵，因此只需计算矩阵的上三角阵的元素即可。将所有上三角阵中的元素按照图5中所示分成4组，分别由4片DSP进行计算。

图5 协方差矩阵计算分配

4.4 平均协方差矩阵特征分解及方位估计

对平均协方差矩阵进行特征分解，采用前面介绍的并行Jacobi方法在FPGA2上实现。方位估计可以采用经典的MUSIC算法空间搜索方法来实现，也可以采用Root-MUSIC算法求根实现。对于空间扫描方法而言，由于前端的波束形成已经初步确定了回波的方位，因此不需要对整个空间进行搜索，而只需对特定角度范围进行搜索即可，用两片DSP实现。对于求根方法，采用两片DSP轮流实现求根过程。

5 总结

从MSB-RMU算法的并行实现的角度入手，针对小规模系统设计的基于FPGA的双边CORDIC结构的并行Jacobi实现方法，为特征分解类算法提供了一种高效的单芯片解决方案，对于类似算法的并行实现提供了一种思路。设计的由6片DSP和2片FPGA组成的并行信号处理平台上较好的实现了MSB-RMU算法的实时处理，但对于大规模的并行处理平台及算法的实时实现，如算法所用资源，工作频率等因素，如何评估需要什么样的硬件平台，在FPGA上完成哪些功能等方面还需要进行更深入的研究。

［1］周天.超宽覆盖海底地形地貌高分辨探测技术研究［D］.哈尔滨工程大学，2005:43-50.

［2］周天，李海森，朱志德，等.多波束测深系统中多子阵检测法的改进及其性能分析［J］.声学技术，2005，24（3）:138-143.

［3］任春丽，徐甲同，王俊平.实对称三对角矩阵特征值的一种并行算法及实现［J］.西安电子科技大学学报，1999，26（2）:217-221.

［4］张丽君，金绥更.向量算法与并行算法［M］.北京:国防工业出版社，1993:123-130.

［5］袁生光.对称矩阵特征值分解的硬件实现研究［D］.浙江大学，2008:30-37.

［6］MINSEOK KIM，KOICHI ICHIGE，HIROYUKI ARAI.Design of Jacobi EVD Processor based on CORDIC for DOA estimation with MUSIC algorithm ［C］.Personal，Indoor and Mobile Radio Communications，2002.The 13th IEEE International Symposium on，Lisboa，Portugal，2002:120-124

［7］李小波，薛王伟，孙志勇.一种求解复Hermite矩阵特征值的方法［J］.数据采集与处理，2005，20（4）: 403-406.

［8］谢伦国.并行计算机互联网络技术［M］.北京: 电子工业出版社，2004:154-169.