超分辨测向中奇异值分解的FPGA实现*

谈文韬,朱业腾,黎仁刚,林 明

(1.江苏科技大学 电子信息学院,江苏 镇江 212003;2.中船重工集团723研究所,江苏 扬州 225101)



超分辨测向中奇异值分解的FPGA实现*

谈文韬1,朱业腾2,黎仁刚2,林 明**1

(1.江苏科技大学 电子信息学院,江苏 镇江 212003;2.中船重工集团723研究所,江苏 扬州 225101)

奇异值分解是超分辨测向技术的核心组成部分，现有的并行实现方案适用范围窄，运算量大，迭代时间长。为了满足测向接收机系统的高实时性需求，结合双边Jacobi算法的交换策略和单边Jacobi算法的求角结构，提出了一种改进的实现方法。该实现方法修正了脉动阵列的收敛性问题，提高了复数矩阵的收敛速度。同时，给出了算法的现场可编程门阵列(FPGA)实现结构。仿真结果证明该方案耗时在百微秒以内，能够应用于电子侦察设备。

电子侦察；超分辨测向；奇异值分解；脉动阵列

1 引 言

由于现代战场电磁环境的日益复杂，超分辨测向技术以其优越的算法性能得到了越来越多的关注[1-2]，奇异值分解(Singular Value Decomposition,SVD)成为了超分辨测向技术实现与应用的掣肘。与在其他方面的应用不同[3-4]，应用于电子对抗的奇异值分解[5]对精度敏感，具有高实时性、高可靠性的严格要求。因为分析问题角度不同，一般应用于该领域的奇异值分解更注重实时性，运行时间通常约束在100 μs以内；而常规实现方法更注重资源消耗，运行时间通常为1～100 ms级[6-7]，与电子侦察系统的支援应用需求相距甚远。并行化脉动(Brent-Lu-Van，BLV)阵列在20世纪80年代被应用于SVD的计算[8]，大大缩短了超分辨测向算法的整体运行时间。Bravo等人[9]于2006年对算法提出了改进，但BLV阵列存在非对称矩阵收敛慢甚至不收敛的情况，且仅适用偶数维度实数矩阵。在工程应用中，有限位宽数据的舍入操作，导致计算过程中的矩阵数值并非完全对称，传统BLV结构会放大该误差的影响，影响算法的计算精度。

本文以多重分类(Multiple Signal Classification,MUSIC)算法中奇异值分解的应用为出发点，在双边Jacobi算法BLV并行化时间结构的基础上，综合考虑系统实时性、可靠性以及资源等需求，给出了一种适用超分辨测向算法的SVD并行化实现结构。

2 MUSIC算法

假设M元任意天线阵列的所有阵元均位于坐标系XOY平面内，第k个阵元坐标为(xk,yk,0)，第i个窄带信号波长为λi，来波方向为(θi,φi)，如图1，则第k个阵元到圆心(即原点)的相位差Δrik为

(1)

图1 天线阵列Fig.1 Antenna array

天线阵列的接收公式为

(2)

x(t)=As(t)+n(t)。

(3)

式中：P表示信号源数目，s(t)表示P个入射信号的集合，x(t)表示M个阵元处接收信号的组合，A表示天线阵的阵列流形。

MUSIC算法基本原理为

Rxx=E[x(t)x(t)H]=ARssAH+σ2In=UsΣsVs+UnΣnVn，

(4)

(5)

(6)

式中：Rxx表示接收数据的协方差矩阵;Rss表示入射信号协方差矩阵;σ2表示噪声功率;In表示单位矩阵;Us、Un分别为信号与噪声协方差矩阵的左奇异矢量，也即信号子空间与噪声子空间;Vs、Vn为右奇异矢量，对于方阵，同样可作为子空间，通常只计算左奇异矢量;上标“^”表示极大似然估计;Pmusic为空间谱。

3 SVD并行实现算法

奇异值分解算法具有较为明确的物理意义，设Rxx为M阵元天线阵列接收信号的自相关矩阵，则Rxx为一个M×M的Hermite矩阵。因为Hermite矩阵是正规矩阵，所以Rxx的特征向量中可找出一组正交基，且特征值均为实数；U与V均为由Rxx特征向量构成的M×M酉矩阵，包含对应信号在不同阵元间的相位差信息。同时，由于噪声与计算误差的存在，Rxx的特征值非零，为非奇异矩阵。

目前应用最为广泛的是双边Jacobi算法[10]与单边Jacobi算法[11]。

双边Jacobi算法是一种通过平面旋转完成消元计算，使矩阵收敛于由奇异值构成的对角矩阵的可并行化算法。

以Rxx为例，双边Jacobi算法可以表示为

(7)

式中：S为矩阵Rxx对应的奇异值对角矩阵，Gk为对应m、n行列交点数值的Givens变换矩阵的一种变形。

应用于FPGA的双边Jacobi算法并行化实现方法，是一种基于BLV脉动阵列的实现结构。该方法由两种2×2的运算单元构成，分别为主对角线上角求取主运算单元与其他位置的旋转乘法单元。

在BLV阵列结构中，每次旋转后进行相邻行列的阵列交换，该交换策略的元素遍历性较差，如图3所示。图中，U12与U21单元的b、c元素无法传递到主对角线上的旋转单元内，仅能通过传递的旋转角进行迭代消元，且计算误差会随着迭代的增加在类似位置的值内累积，导致算法收敛慢，计算结果误差较大。

图2 传统BLV阵列结构Fig.2 Structure of traditional BLV array

单边算法旋转角度基于行列正交性计算得出，这与双边算法的旋转消元方法有一定差异。本文引入了单边Jacobi算法的旋转角计算方法,其第m列与第n列元素间旋转角的计算公式为

(8)

(9)

式中:Rcol(*)表示矩阵对应列，sign(*)表示取符号。式(8)与式(9)对复数矩阵同样适用。

4 实现方案选择

针对双边Jacobi算法中BLV阵列只适用奇数维度实数矩阵、元素遍历性差,以及单边Jacobi算法收敛速度慢的问题，结合两种算法的计算特点，下面给出符合应用背景的实现方案。

4.1 双边算法交换策略

传统BLV阵列具有元素遍历性差的特点。在运算过程中，两个主运算单元副对角线上所有乘法单元内的元素收敛速度慢，甚至出现发散的情况;而双边算法的行列交换策略具有最优的元素遍历性，算法收敛速度快，结果可靠。

两两组合主对角线上元素，共有M(M-1)/2种组合；每次向BLV中输入⎣M/2」个不含重复列的组合，并按组合位置交换对应行列；按此方式进行一次计算，称作一次向量正交化，完成所有正交化组合称作“一扫”，完成一扫需要2⎣(M-1)/2」+1次并行向量正交化计算。将生成的交换策略转变为矩阵的地址，每个数据有行地址与列地址两个坐标信息，所以完成一扫所需的地址数据量L可由式(10)计算：

L=2M2(2⎣(M-1)/2」+1)。

(10)

每个地址的位宽为「lb(M+1)⎤，「*⎤表示向上取整。

组合的行列排列信息可以离线计算，并按矩阵元素地址格式量化，以常数的形式存储，不会增加计算量；或者，通过设计一个地址发生器，进行行列的遍历。

4.2 复数计算方案

常规的复数矩阵分解方法是将维度为M的复数矩阵化为维度2M的实对称矩阵，采用实数矩阵方法进行分解;而奇异值分解的计算复杂度为O(M3)，算法的复杂度变为原来的8倍以上，尤其是在需要计算特征向量时，大大增加了计算成本与实现难度。同时使用实数化处理时，增加了每一扫需要的向量正交化次数，也增加了算法的迭代时间。

此处采用单边Jacobi算法的旋转角计算方式，根据复数计算与现场可编程逻辑门阵列(Field Programmable Gate Array,FPGA)结构对公式进行变形优化。将主运算单元中的2×2矩阵进行旋转消元，需要满足

(11)

因为Rxx为Hermite复数矩阵，所以a、d为实数，b、c相互共轭。根据式(8)和式(9)，复数旋转角求取公式化简为

(12)

式中:“*′”表示复数的共轭。式(12)结构更适于FPGA实现。同时不再向同行同列传递旋转角，而是经过复数修正的旋转矩阵，即先乘以实数化矩阵的共轭转置将矩阵还原：

(13)

4.3 奇数行(列)乘法单元

在实际运算过程中，由于奇数末位行位没有配对组合，不必计算旋转角，但在旋转相乘时，必须连接旋转矩阵传递相乘，来保证行列变换的一致性。

本文通过增加单边末行(列)单边乘法运算模块，如图3虚线框中部分，对奇数余行(列)的运算进行补充:

(14)

图3 修正的BLV阵列Fig.3 Modified BLV array

最终的奇异值分解并行实现方案如图2所示。该算法结构中，误差累积因素只有三角函数正、余弦匹配误差，其他过程误差仅影响收敛速度，与计算精度无关，保证了数值的稳定性。

4.4 实现方案性能分析

此处分别对传统BLV实现方案与本文实现方案进行分析对比。

(15)

以9元天线阵接收数据实部的协方差矩阵作为BLV阵列的输入数据，则ek随正交化迭代次数的变化曲线如图4和图5所示。

图4 传统BLV方案的曲线Fig.4 Curve of traditional BLV solution

图5 本文方案的曲线Fig.5 Curve of the proposed solution

从图4中可以看出，即便算法的输入数据维度M相对较小(M<10)，BLV阵列的部分副对角线上的值无法消元，导致奇异值的均方误差不能收敛于0值，与第3节中的描述一致,因此,该实现方法无法满足超分辨测向算法的应用需求。

以9元天线阵接收复数矩阵作为本文算法的输入数据，在图5中，可以发现，均方误差值的下降曲线较为平滑，且在放大框图中可以看出，到第3扫(第27次正交变换)完成时，均方误差已经接近10-5。

试验结果说明了本文的方案选择，收敛速度更快，稳健性更好。该性能分析的过程中部分数据为双精度浮点数，比起FPGA内单精度块浮点的实现方法存在可控误差。

5 硬件实现方案

传统BLV阵列的主运算模块至少需要3个坐标旋转计算法(Coordinate Rotation Digital Computer,CORDIC)模块，CORDIC算法模块周期长，大大增加了算法输出的延时。此处使用组合逻辑与查找表方式分别设计了2个CORDIC模块，大大缩减了旋转角求取时间。旋转角与sp矩阵计算结构相似，求取顺序不同，通过一次复用减少了算法的资源消耗。旋转矩阵求取模块实现结构如图6所示。

图6 旋转角求取模块结构Fig.6 Structure of rotation angle calculation module

基于双边算法交换策略的奇异值分解迭代流程如图7所示。

图7 奇异值分解迭代流程Fig.7 Iterative process of SVD

本文以200 MHz时钟、块浮点数据结构(符号位1,指数位宽8,溢出保护位1，数据位宽24)进行算法硬件实现。每个旋转角计算单元需要 7 421个查找表(Look-up Table,LUT)、790个触发器(Flip-flop,FF)、16个DSP48E1(一种FPGA乘法器)，每个乘法单元需要858个LUT、438个FF、40个DSP48E1。

以ε=10-6为向量正交精度，完成一次9行9列的奇异值分解计算仅需3 172个周期，约15.86 μs。验证结果证明该实现方案优于文献[4-5]中的毫秒级实现方案。

6 结 论

本文以双边Jacobi算法的BLV实现方案为基础，结合单边算法与奇数维度矩阵实现方法的优点，完成了一种遍历性、稳健性、收敛速度均较优的实现方案。该方案耗时被约束在百μs以内，远低于原有的毫秒级实现方案，使超分辨测向算法在高实时性电子侦察系统中的应用成为可能。本文的实现方法消耗的片上资源较多，仍有进一步优化的空间。

[1] 胡子扬,任渊.一种最小冗余线阵的目标DOA估计方法[J].电讯技术,2014,54(11):1493-1498. HU Ziyang,REN Yuan. A DOA estimation method based on sub-minimum redundancy linear array[J].Telecommunication Engineering,2014,54(11):1493-1498.(in Chinese)

[2] 付淑娟,景小荣,张祖凡,等.基于虚拟阵列改进MUSIC算法的相干信源DoA估计[J].电讯技术,2011,51(11):63-67. FU Shujuan,JING Xiaorong,ZHANG Zufan,et al.DoA estimation of coherent sources by using virtual array-based improved MUSIC algorithm[J].Telecommunication Engineering,2011,51(11):63-67.(in Chinese)

[3] MOHANTY R,ANIRUDH G,PRADHAN T. Design and performance analysis of fixed-point Jacobi SVD algorithm on reconfigurable system[J].IERI Procedia,2014(7):21-27.

[4] QIAO H L. New SVD based initialization strategy for non-negative matrix factorization[J].Pattern Recognition Letters,2015,63(C):71-77.

[5] LIU Y,CUI H Y. Antenna array signal direction of arrival estimation on digital signal processor(DSP)[J].Procedia Computer Science,2015,55(7):782-791.

[6] MILFORD D,SANDELL M. Singular value decomposition using an array of CORDIC processors[J].Signal Processing,2014,102(9):163-170.

[7] TAI Y G,DAN C T,PSARRIS K. Scalable matrix decompositions with multiple cores on FPGAs[J].Microprocessors and Microsystems,2013,37(8):887-898.

[8] LUK F T,BRENT R P. The solution of singular-value and symmetric eigenvalue problems on multiprocessor arrays[J].SIAM Journal on Scientific and Statistical Computing,1985,6(1)：69-84．

[9] BRAVO I,JIMNEZ P,MAZO M,et al.Implementation in FPGAs of Jacobi method to solve the eigenvalue and eigenvector problem[C]//Proceedings of 2006 International Conference on Field Programmable Logic & Applications.Madrid,Spain:IEEE,2006:1-4.

[10] JAMES D,KRESIMIR V. Jacobi's method is more accurate than QR[J].SIAM Journal on Matrix Analysis and Applications,1992,13(4):1204-1245.

[11] 郭强,赵雷. 一种基于动态序列的单边Jacobi方法[J].苏州大学学报(工科版),2011,31(4):16-22. GUO Qiang,ZHAO Lei. A new one-side Jacobi basedon dynamic ordering[J].Journal of Soochow University(Engineering Science Edition),2011,31(4):16-22.(in Chinese)

FPGA Implementation of Singular Value Decomposition Algorithm in Super-resolution Direction-finding

TAN Wentao1,ZHU Yeteng2,LI Rengang2,LIN Ming1

(1.School of Electronics Information,Jiangsu University of Science and Technology,Zhenjiang 212003,China;2.No.723 Research Institute of China Shipbuilding Industry Corporation,Yangzhou 225101,China)

Singular value decomposition(SVD) is the core of super-resolution direction-finding(DF) technique.The existing parallel implementation methods have narrow applicability,intensive computation and lengthy iteration. In order to meet the application requirements of DF receiver system,a modified implementation method is proposed according to the exchange policy of bilateral Jacobi algorithm and the angle calculation structure of unilateral Jacobi algorithm. In this implementation method,the ergodic problem in Brent-Lu-Van(BLV) array is modified,and the convergence rate in complex matrix is improved.An algorithm structure on field programmable gate array(FPGA) is presented.Simulation results show that the method consumes less than 100 ms and can be applied to electronic reconnaissance equipment.

electronic reconnaissance;super-resolution direction-finding;sigular value decomposition(SVD);Brent-Lu-Van(BLV) array

10.3969/j.issn.1001-893x.2017.07.012

谈文韬,朱业腾,黎仁刚,等.超分辨测向中奇异值分解的FPGA实现[J].电讯技术，2017,57(7):801-805.[TAN Wentao,ZHU Yeteng,LI Rengang,et al.FPGA implementation of singular value decomposition algorithm in super-resolution direction-finding[J].Telecommunication Engineering,2017,57(7):801-805.]

2016-10-26；

2017-02-22 Received date:2016-10-26；Revised date：2017-02-22

国家自然科学基金青年科学基金项目(61401179)

TN971

A

1001-893X(2017)07-0801-05

谈文韬(1991—),男,江苏连云港人,硕士研究生,主要研究方向为雷达信号与信息处理理论与技术;

Email:DorusTan@163.com

朱业腾(1988—),男,江苏扬州人,工程师,主要研究方向为数字信号处理;

黎仁刚(1978—),男,江苏扬州人,博士,研究员,主要研究方向为阵列信号处理;

林 明(1960—),男,辽宁大连人,教授、硕士生导师,主要研究方向为雷达信号处理。

Email:jskdlm@qq.com

**通信作者：jskdlm@qq.com Corresponding author：jskdlm@qq.com