对称矩阵特征值分解的FPGA实现
刘永勤
摘 要： 针对应用于MUSIC DOA估计的数据协方差矩阵特征值分解的需要，给出一个特征值分解的硬件实现方案，并阐述了基本思想。设计采用基于CORDIC的Jacobi算法实现实对称矩阵特征值分解，并在FPGA上对5×5矩阵进行了硬件仿真，经过理论分析和实验验证，该设计可以计算出全部特征值和特征向量，为MUSIC算法的FPGA实现奠定了基础。
关键词： MUSIC算法； 特征值分解； Jacobi算法； CORDIC算法； FPGA
中图分类号： TN911?34； TN929.1 文献标识码： A 文章编号： 1004?373X（2017）12?0015?04
Abstract： Aiming at the needs of the data covariance matrix eigenvalue decomposition used in DOA estimation such as MUSIC， a hardware implementation scheme of the eigenvalue decomposition is provided and the basic idea is described in this paper. The Jacobi algorithm based on CORDIC is adopted in the design to achieve real symmetric matrix eigenvalue decomposition， and conduct the hardware emulation for 5×5 matrix in FPGA. The results of theoretical analysis and experimental verification show that the design can calculate all eigenvalues ??and eigenvectors， and has laid the foundation for FPGA implementation of MUSIC algorithm.
Keywords： MUSIC algorithm； eigenvalue decomposition； Jacobi algorithm； CORDIC algorithm； FPGA
0 引 言
多信号分类（MUSIC）[1]算法是波达方向（DOA）估计技术中最具代表性的高分辨力算法之一，因其突破了传统方法的瑞利极限而广受人们青睐。虽然MUSIC算法理论上已经成熟，但运算量大，不利于其硬件实现，对接收矢量的数据协方差矩阵进行特征值分解（EVD）占据了MUSIC算法60%的运算量，因此解决EVD在硬件上的实现是硬件实现MUSIC算法的关键所在，同时EVD作为一个相对独立的部分，可以广泛应用到人工视觉、电力电子、机械力学系统等其他众多领域中去。
Jacobi方法是计算实对称矩阵特征值的常用方法，它涉及到开方、除法等非线性运算，这对硬件实现有一定难度， Cavallaro和Luk直接利用CORDIC的两种计算模式解决Jacobi中的反正切及坐标旋转计算问题[2]。为进一步减小EVD计算量，Yang和Bohme应用了双平面旋转方法[3]，通过两个并行矩阵的并行计算实现特征值分解。目前，Jacobi算法的研究集中在对所有元素的并行计算上[4]，Brent，Luk和Van Loan提出了采用脉动阵列的結构实现矩阵特征值求解的方法[5]，充分体现了Jacobi的高度并行性，是一种经典的阵列结构。A.Ahmedsaid等人在文献[6]中对上述阵列结构进行了改进，并提出了在FPGA上的实现方案。并行计算虽然可以提高EVD计算速度，但是这种方法过程复杂，且占用资源多，不适合在单处理器系统上实现。因此本研究利用基于CORDIC的顺序Jacobi法在单片FPGA上实现EVD计算，同时针对经典Jacobi算法单次扫描时间长这一不足，简化操作流程，并通过合理设计算法结构来压缩执行时间。
1 由复数转实数
实际中经常会遇到要求解复Hermite矩阵特征值和特征向量的问题，但复矩阵的特征值分解的计算量很大，计算复杂度[7]高达，且复对称矩阵不具有实对称矩阵的任何重要性质。而经典的Jacobi方法是针对实对称矩阵的特征值分解方法，因此为了简化计算应用Jacobi方法，需要先把复Hermite矩阵转化为实对称矩阵[8?9]。
对中心有一个阵元的元均匀圆阵，将天线阵列流形按中心共轭对称形式排列，其阵列结构见图1。
2 Jacobi方法
Jacobi方法是一种经典的求解实对称矩阵特征值的计算方法之一。假设原实对称的协方差矩阵为，则Jacobi方法通过构造一系列的平面旋转矩阵，使原矩阵通过有限次的平面旋转变换转化为对角阵，对角阵的主对角线元素即为的特征值[10?11]。严格的说产生对角型所需的旋转变换是无穷的，但在实际计算中当非对角元素小到可以忽略，就可以认为矩阵已经被对角化。
对作一系列平面旋转变换，有：
3 CORDIC算法
从前文可以看出Jacobi方法涉及到矩阵旋转和反正切等复杂计算，这些计算在FPGA上是难以实现的，而CORDIC算法正是解决这一问题的有利工具。CORDIC算法是J.Voider等人于1959年提出的，它是一种循环迭代算法，通过用一系列与运算基数有关的基本角的不断偏摆来逼近所需要旋转的角度，它将很多复杂的函数运算化简为移位和加法的迭代，提供了一种适合硬件的实现方法[13?14]。其圆周坐标下的基本迭代公式如下：
其中每次旋转的基本角取，为每次旋转的方向，是根据角度累积量来判定的，即，则。CORDIC算法可以工作在旋转模式和矢量模式，旋转模式用来完成坐标旋转运算，矢量模式用来完成反正切计算得到最佳旋转角，有。在旋转模式情况下，当， 经过n次基本角旋转迭代后旋转结果为：
式中：为对结果进行修正的比例因子，。通过比较式（12）和式（6）可以很明了地看出CORDIC算法可以解决Jacobi的旋转问题。
4 FPGA实现方案
由于待分解矩阵是实对称矩阵，因此在FPGA中只需要处理其上三角元素就可以取代整个矩阵，扫描顺序为。在经典Jacobi算法中， 一次Jacobi扫描要进行1次CORDIC反正切计算和2次CORDIC旋转才能完成，单次扫描时间长，针对这一不足，具体实现时对其进行简化。
由于Jacobi中，将根据的正负性判定，从而确定符号集，符号集一旦确定就可以根据式（11）的迭代实现CORDIC旋转，省去了反正切计算的过程。
令，由可以得到：
根据式（9）可知应取初值为，，这样根据此更新公式和初值就可以迭代求得符号集。
将式（6）～式（8）展开可以得到：
则当符号集确定可以将符号集作为地址，通过查表的方式得到，从而实现对角元素式（13）的更新，而对于要扫描的非对角元素可以直接置零。式（14）中行（列）行（列）其他元素的更新则要用1次CORDIC旋转来实现，初值为。这里符号计算和CORDIC旋转可以并行进行，每算出一个值就进行一个基本角的旋转。总体结构如图2所示。
5 设计实现及结果
本设计所选用的器件是ALTERA公司的EP3C120F780C7器件，使用Verilog语言完成电路描述之后，在Quartus Ⅱ软件平台上进行编译，在Modelsim中进行仿真。以矩阵为例，采用21 b定点数，将如下矩阵作为输入矩阵进行系统仿真。
将所得定点数结果化为对应的浮点数数据在表1中给出，在Matlab中直接调用eig函数计算出的特征值和特征向量在表2中给出。表1和表2分别表示设计结果与理论结果。进行比较可以看出，本设计所得结果与理论结果基本一致；虽然第4列和第5列特征向量的符号相反，但由于各特征向量间的乘积为一个很小的接近于0的数，故其正负不影响其正交性。
假设有任意两个信号入射到五元均匀圆阵上，入射方位角分别为120°和160°，图3给出了用MUSIC算法进行DOA估计的结果，并对局部进行了放大。其中，图3中（a）为理论估计结果，其EVD计算部分采用Matkab浮点计算，图3中（b）为本文设计估计结果，其EVD计算部分采用本设计的定点操作，除了EVD计算，所有MUSIC算法中的其他计算都采用浮点操作。从图中可以看出，本设计可以有效估计出入射信號的DOA信息，证明此硬件结构可以应用于MUSIC算法的FPGA实现中。
6 结 语
本文采用CORDIC算法和Jacobi算法，在FPGA上实现了数据协方差矩阵的特征值和特征向量计算，为MUSIC算法的FPGA实现提供了保证。实对称矩阵特征分解应用广泛，其硬件实现对其他算法的工程应用有重要作用。
参考文献
[1] SCHMIDT R. Multiple emitter location and signal parameter estimation [J]. IEEE transactions on antennas and propagation， 1986， AP?34（3）： 276?280.
[2] CAVALLARO J R， LUK F T. CORDIC Arithmetic for an SVD Processor [J]. Journal of parallel and distributed computing， 1988， 5（3）： 271?290.
[3] YANG B， BOHME J F. Reducing the computations of the SVD array given by Brent and Luk [J]. SPIE advanced algorithm and architectures for signal processing IV， 1989， 1152： 92?102.
[4] GOTZE J. Parallel methods for iterative matrix decompositions [C]// IEEE International Sympoisum on Circuits and Systems. [S.l.]： IEEE， 1991： 232?235.
[5] BRENT R P， LUK F T， VAN LOAN C. Computation of the generalized singular value decomposition using mesh?connected processors [C]// Proceedings of SPIE Sympoisum on Real?Time Signal Processing， 1983， 431： 66?72.
[6] AHMEDSAID A， AMIRA A， BOURIDANE A. Improved SVD systolic array and implementation on FPGA [C]// Proceedings of 2003 IEEE International Conference on Field?Programmable Technology. [S.l.]： IEEE， 2003： 35?42.
[7] LINEBARGER D A， DEGROAT R D， DOWLING E M. Efficient direction?finding methods employing forward?backward averaging [J]. IEEE transactions on signal processing， 1994， 42（8）： 2136?2145.
[8] KIM Minseok， ICHIGE Koichi， ARAI Hiroyuki. implementation of FPGA based fast DOA estimator using unitary MUSIC algorithm [C]// Proceedings of 58th IEEE Vehicular Technology Conference. [S.l.]： IEEE， 2003： 213?217.
[9] 徐德琛，刘志文，徐友根.某测向系统中MUSIC算法的FPGA实现[J].北京理工大学学报，2010，30（9）：1107?1111.
[10] 于正文，尹庆莉.求特征值的Jacobi方法[J].山东科学，2011，24（6）：19?21.
[11] WILKINGSON J H. The algebraic eigenvalue problem [M]. UK： Oxford Science Publications， 1999.
[12] BRAVO I， JIMENEZ 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 and Applications. Madrid， Spain： IEEE， 2006： 301?316.
[13] 杨宏，李国辉，刘立新.基于FPGA的CORDIC算法的实现[J].西安邮电学院学报，2008，13（1）：75?77.
[14] KIM Minseok， ICHIGE K.， ARAI H. Design of Jacobi EVD processor based on CORDIC for DOA estimation with MUSIC algorithm [C]// Proceedings of 2002. The 13th IEEE International Symposium on Personal， Indoor and Mobile Radio Communications. Lisboa， Portugal： IEEE， 2002： 120?124.
摘 要： 针对应用于MUSIC DOA估计的数据协方差矩阵特征值分解的需要，给出一个特征值分解的硬件实现方案，并阐述了基本思想。设计采用基于CORDIC的Jacobi算法实现实对称矩阵特征值分解，并在FPGA上对5×5矩阵进行了硬件仿真，经过理论分析和实验验证，该设计可以计算出全部特征值和特征向量，为MUSIC算法的FPGA实现奠定了基础。
关键词： MUSIC算法； 特征值分解； Jacobi算法； CORDIC算法； FPGA
中图分类号： TN911?34； TN929.1 文献标识码： A 文章编号： 1004?373X（2017）12?0015?04
Abstract： Aiming at the needs of the data covariance matrix eigenvalue decomposition used in DOA estimation such as MUSIC， a hardware implementation scheme of the eigenvalue decomposition is provided and the basic idea is described in this paper. The Jacobi algorithm based on CORDIC is adopted in the design to achieve real symmetric matrix eigenvalue decomposition， and conduct the hardware emulation for 5×5 matrix in FPGA. The results of theoretical analysis and experimental verification show that the design can calculate all eigenvalues ??and eigenvectors， and has laid the foundation for FPGA implementation of MUSIC algorithm.
Keywords： MUSIC algorithm； eigenvalue decomposition； Jacobi algorithm； CORDIC algorithm； FPGA
0 引 言
多信号分类（MUSIC）[1]算法是波达方向（DOA）估计技术中最具代表性的高分辨力算法之一，因其突破了传统方法的瑞利极限而广受人们青睐。虽然MUSIC算法理论上已经成熟，但运算量大，不利于其硬件实现，对接收矢量的数据协方差矩阵进行特征值分解（EVD）占据了MUSIC算法60%的运算量，因此解决EVD在硬件上的实现是硬件实现MUSIC算法的关键所在，同时EVD作为一个相对独立的部分，可以广泛应用到人工视觉、电力电子、机械力学系统等其他众多领域中去。
Jacobi方法是计算实对称矩阵特征值的常用方法，它涉及到开方、除法等非线性运算，这对硬件实现有一定难度， Cavallaro和Luk直接利用CORDIC的两种计算模式解决Jacobi中的反正切及坐标旋转计算问题[2]。为进一步减小EVD计算量，Yang和Bohme应用了双平面旋转方法[3]，通过两个并行矩阵的并行计算实现特征值分解。目前，Jacobi算法的研究集中在对所有元素的并行计算上[4]，Brent，Luk和Van Loan提出了采用脉动阵列的結构实现矩阵特征值求解的方法[5]，充分体现了Jacobi的高度并行性，是一种经典的阵列结构。A.Ahmedsaid等人在文献[6]中对上述阵列结构进行了改进，并提出了在FPGA上的实现方案。并行计算虽然可以提高EVD计算速度，但是这种方法过程复杂，且占用资源多，不适合在单处理器系统上实现。因此本研究利用基于CORDIC的顺序Jacobi法在单片FPGA上实现EVD计算，同时针对经典Jacobi算法单次扫描时间长这一不足，简化操作流程，并通过合理设计算法结构来压缩执行时间。
1 由复数转实数
实际中经常会遇到要求解复Hermite矩阵特征值和特征向量的问题，但复矩阵的特征值分解的计算量很大，计算复杂度[7]高达，且复对称矩阵不具有实对称矩阵的任何重要性质。而经典的Jacobi方法是针对实对称矩阵的特征值分解方法，因此为了简化计算应用Jacobi方法，需要先把复Hermite矩阵转化为实对称矩阵[8?9]。
对中心有一个阵元的元均匀圆阵，将天线阵列流形按中心共轭对称形式排列，其阵列结构见图1。
2 Jacobi方法
Jacobi方法是一种经典的求解实对称矩阵特征值的计算方法之一。假设原实对称的协方差矩阵为，则Jacobi方法通过构造一系列的平面旋转矩阵，使原矩阵通过有限次的平面旋转变换转化为对角阵，对角阵的主对角线元素即为的特征值[10?11]。严格的说产生对角型所需的旋转变换是无穷的，但在实际计算中当非对角元素小到可以忽略，就可以认为矩阵已经被对角化。
对作一系列平面旋转变换，有：
3 CORDIC算法
从前文可以看出Jacobi方法涉及到矩阵旋转和反正切等复杂计算，这些计算在FPGA上是难以实现的，而CORDIC算法正是解决这一问题的有利工具。CORDIC算法是J.Voider等人于1959年提出的，它是一种循环迭代算法，通过用一系列与运算基数有关的基本角的不断偏摆来逼近所需要旋转的角度，它将很多复杂的函数运算化简为移位和加法的迭代，提供了一种适合硬件的实现方法[13?14]。其圆周坐标下的基本迭代公式如下：
其中每次旋转的基本角取，为每次旋转的方向，是根据角度累积量来判定的，即，则。CORDIC算法可以工作在旋转模式和矢量模式，旋转模式用来完成坐标旋转运算，矢量模式用来完成反正切计算得到最佳旋转角，有。在旋转模式情况下，当， 经过n次基本角旋转迭代后旋转结果为：
式中：为对结果进行修正的比例因子，。通过比较式（12）和式（6）可以很明了地看出CORDIC算法可以解决Jacobi的旋转问题。
4 FPGA实现方案
由于待分解矩阵是实对称矩阵，因此在FPGA中只需要处理其上三角元素就可以取代整个矩阵，扫描顺序为。在经典Jacobi算法中， 一次Jacobi扫描要进行1次CORDIC反正切计算和2次CORDIC旋转才能完成，单次扫描时间长，针对这一不足，具体实现时对其进行简化。
由于Jacobi中，将根据的正负性判定，从而确定符号集，符号集一旦确定就可以根据式（11）的迭代实现CORDIC旋转，省去了反正切计算的过程。
令，由可以得到：
根据式（9）可知应取初值为，，这样根据此更新公式和初值就可以迭代求得符号集。
将式（6）～式（8）展开可以得到：
则当符号集确定可以将符号集作为地址，通过查表的方式得到，从而实现对角元素式（13）的更新，而对于要扫描的非对角元素可以直接置零。式（14）中行（列）行（列）其他元素的更新则要用1次CORDIC旋转来实现，初值为。这里符号计算和CORDIC旋转可以并行进行，每算出一个值就进行一个基本角的旋转。总体结构如图2所示。
5 设计实现及结果
本设计所选用的器件是ALTERA公司的EP3C120F780C7器件，使用Verilog语言完成电路描述之后，在Quartus Ⅱ软件平台上进行编译，在Modelsim中进行仿真。以矩阵为例，采用21 b定点数，将如下矩阵作为输入矩阵进行系统仿真。
将所得定点数结果化为对应的浮点数数据在表1中给出，在Matlab中直接调用eig函数计算出的特征值和特征向量在表2中给出。表1和表2分别表示设计结果与理论结果。进行比较可以看出，本设计所得结果与理论结果基本一致；虽然第4列和第5列特征向量的符号相反，但由于各特征向量间的乘积为一个很小的接近于0的数，故其正负不影响其正交性。
假设有任意两个信号入射到五元均匀圆阵上，入射方位角分别为120°和160°，图3给出了用MUSIC算法进行DOA估计的结果，并对局部进行了放大。其中，图3中（a）为理论估计结果，其EVD计算部分采用Matkab浮点计算，图3中（b）为本文设计估计结果，其EVD计算部分采用本设计的定点操作，除了EVD计算，所有MUSIC算法中的其他计算都采用浮点操作。从图中可以看出，本设计可以有效估计出入射信號的DOA信息，证明此硬件结构可以应用于MUSIC算法的FPGA实现中。
6 结 语
本文采用CORDIC算法和Jacobi算法，在FPGA上实现了数据协方差矩阵的特征值和特征向量计算，为MUSIC算法的FPGA实现提供了保证。实对称矩阵特征分解应用广泛，其硬件实现对其他算法的工程应用有重要作用。
参考文献
[1] SCHMIDT R. Multiple emitter location and signal parameter estimation [J]. IEEE transactions on antennas and propagation， 1986， AP?34（3）： 276?280.
[2] CAVALLARO J R， LUK F T. CORDIC Arithmetic for an SVD Processor [J]. Journal of parallel and distributed computing， 1988， 5（3）： 271?290.
[3] YANG B， BOHME J F. Reducing the computations of the SVD array given by Brent and Luk [J]. SPIE advanced algorithm and architectures for signal processing IV， 1989， 1152： 92?102.
[4] GOTZE J. Parallel methods for iterative matrix decompositions [C]// IEEE International Sympoisum on Circuits and Systems. [S.l.]： IEEE， 1991： 232?235.
[5] BRENT R P， LUK F T， VAN LOAN C. Computation of the generalized singular value decomposition using mesh?connected processors [C]// Proceedings of SPIE Sympoisum on Real?Time Signal Processing， 1983， 431： 66?72.
[6] AHMEDSAID A， AMIRA A， BOURIDANE A. Improved SVD systolic array and implementation on FPGA [C]// Proceedings of 2003 IEEE International Conference on Field?Programmable Technology. [S.l.]： IEEE， 2003： 35?42.
[7] LINEBARGER D A， DEGROAT R D， DOWLING E M. Efficient direction?finding methods employing forward?backward averaging [J]. IEEE transactions on signal processing， 1994， 42（8）： 2136?2145.
[8] KIM Minseok， ICHIGE Koichi， ARAI Hiroyuki. implementation of FPGA based fast DOA estimator using unitary MUSIC algorithm [C]// Proceedings of 58th IEEE Vehicular Technology Conference. [S.l.]： IEEE， 2003： 213?217.
[9] 徐德琛，刘志文，徐友根.某测向系统中MUSIC算法的FPGA实现[J].北京理工大学学报，2010，30（9）：1107?1111.
[10] 于正文，尹庆莉.求特征值的Jacobi方法[J].山东科学，2011，24（6）：19?21.
[11] WILKINGSON J H. The algebraic eigenvalue problem [M]. UK： Oxford Science Publications， 1999.
[12] BRAVO I， JIMENEZ 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 and Applications. Madrid， Spain： IEEE， 2006： 301?316.
[13] 杨宏，李国辉，刘立新.基于FPGA的CORDIC算法的实现[J].西安邮电学院学报，2008，13（1）：75?77.
[14] KIM Minseok， ICHIGE K.， ARAI H. Design of Jacobi EVD processor based on CORDIC for DOA estimation with MUSIC algorithm [C]// Proceedings of 2002. The 13th IEEE International Symposium on Personal， Indoor and Mobile Radio Communications. Lisboa， Portugal： IEEE， 2002： 120?124.
