對(duì)稱矩陣特征值分解的FPGA實(shí)現(xiàn)

劉永勤

摘 要： 針對(duì)應(yīng)用于MUSIC DOA估計(jì)的數(shù)據(jù)協(xié)方差矩陣特征值分解的需要，給出一個(gè)特征值分解的硬件實(shí)現(xiàn)方案，并闡述了基本思想。設(shè)計(jì)采用基于CORDIC的Jacobi算法實(shí)現(xiàn)實(shí)對(duì)稱矩陣特征值分解，并在FPGA上對(duì)5×5矩陣進(jìn)行了硬件仿真，經(jīng)過理論分析和實(shí)驗(yàn)驗(yàn)證，該設(shè)計(jì)可以計(jì)算出全部特征值和特征向量，為MUSIC算法的FPGA實(shí)現(xiàn)奠定了基礎(chǔ)。

關(guān)鍵詞： MUSIC算法； 特征值分解； Jacobi算法； CORDIC算法； FPGA

中圖分類號(hào)： TN911?34； TN929.1 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 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 引 言

多信號(hào)分類（MUSIC）[1]算法是波達(dá)方向（DOA）估計(jì)技術(shù)中最具代表性的高分辨力算法之一，因其突破了傳統(tǒng)方法的瑞利極限而廣受人們青睞。雖然MUSIC算法理論上已經(jīng)成熟，但運(yùn)算量大，不利于其硬件實(shí)現(xiàn)，對(duì)接收矢量的數(shù)據(jù)協(xié)方差矩陣進(jìn)行特征值分解（EVD）占據(jù)了MUSIC算法60%的運(yùn)算量，因此解決EVD在硬件上的實(shí)現(xiàn)是硬件實(shí)現(xiàn)MUSIC算法的關(guān)鍵所在，同時(shí)EVD作為一個(gè)相對(duì)獨(dú)立的部分，可以廣泛應(yīng)用到人工視覺、電力電子、機(jī)械力學(xué)系統(tǒng)等其他眾多領(lǐng)域中去。

Jacobi方法是計(jì)算實(shí)對(duì)稱矩陣特征值的常用方法，它涉及到開方、除法等非線性運(yùn)算，這對(duì)硬件實(shí)現(xiàn)有一定難度， Cavallaro和Luk直接利用CORDIC的兩種計(jì)算模式解決Jacobi中的反正切及坐標(biāo)旋轉(zhuǎn)計(jì)算問題[2]。為進(jìn)一步減小EVD計(jì)算量，Yang和Bohme應(yīng)用了雙平面旋轉(zhuǎn)方法[3]，通過兩個(gè)并行矩陣的并行計(jì)算實(shí)現(xiàn)特征值分解。目前，Jacobi算法的研究集中在對(duì)所有元素的并行計(jì)算上[4]，Brent，Luk和Van Loan提出了采用脈動(dòng)陣列的結(jié)構(gòu)實(shí)現(xiàn)矩陣特征值求解的方法[5]，充分體現(xiàn)了Jacobi的高度并行性，是一種經(jīng)典的陣列結(jié)構(gòu)。A.Ahmedsaid等人在文獻(xiàn)[6]中對(duì)上述陣列結(jié)構(gòu)進(jìn)行了改進(jìn)，并提出了在FPGA上的實(shí)現(xiàn)方案。并行計(jì)算雖然可以提高EVD計(jì)算速度，但是這種方法過程復(fù)雜，且占用資源多，不適合在單處理器系統(tǒng)上實(shí)現(xiàn)。因此本研究利用基于CORDIC的順序Jacobi法在單片F(xiàn)PGA上實(shí)現(xiàn)EVD計(jì)算，同時(shí)針對(duì)經(jīng)典Jacobi算法單次掃描時(shí)間長這一不足，簡化操作流程，并通過合理設(shè)計(jì)算法結(jié)構(gòu)來壓縮執(zhí)行時(shí)間。

1 由復(fù)數(shù)轉(zhuǎn)實(shí)數(shù)

實(shí)際中經(jīng)常會(huì)遇到要求解復(fù)Hermite矩陣特征值和特征向量的問題，但復(fù)矩陣的特征值分解的計(jì)算量很大，計(jì)算復(fù)雜度[7]高達(dá)，且復(fù)對(duì)稱矩陣不具有實(shí)對(duì)稱矩陣的任何重要性質(zhì)。而經(jīng)典的Jacobi方法是針對(duì)實(shí)對(duì)稱矩陣的特征值分解方法，因此為了簡化計(jì)算應(yīng)用Jacobi方法，需要先把復(fù)Hermite矩陣轉(zhuǎn)化為實(shí)對(duì)稱矩陣[8?9]。

對(duì)中心有一個(gè)陣元的元均勻圓陣，將天線陣列流形按中心共軛對(duì)稱形式排列，其陣列結(jié)構(gòu)見圖1。

2 Jacobi方法

Jacobi方法是一種經(jīng)典的求解實(shí)對(duì)稱矩陣特征值的計(jì)算方法之一。假設(shè)原實(shí)對(duì)稱的協(xié)方差矩陣為，則Jacobi方法通過構(gòu)造一系列的平面旋轉(zhuǎn)矩陣，使原矩陣通過有限次的平面旋轉(zhuǎn)變換轉(zhuǎn)化為對(duì)角陣，對(duì)角陣的主對(duì)角線元素即為的特征值[10?11]。嚴(yán)格的說產(chǎn)生對(duì)角型所需的旋轉(zhuǎn)變換是無窮的，但在實(shí)際計(jì)算中當(dāng)非對(duì)角元素小到可以忽略，就可以認(rèn)為矩陣已經(jīng)被對(duì)角化。

對(duì)作一系列平面旋轉(zhuǎn)變換，有：

3 CORDIC算法

從前文可以看出Jacobi方法涉及到矩陣旋轉(zhuǎn)和反正切等復(fù)雜計(jì)算，這些計(jì)算在FPGA上是難以實(shí)現(xiàn)的，而CORDIC算法正是解決這一問題的有利工具。CORDIC算法是J.Voider等人于1959年提出的，它是一種循環(huán)迭代算法，通過用一系列與運(yùn)算基數(shù)有關(guān)的基本角的不斷偏擺來逼近所需要旋轉(zhuǎn)的角度，它將很多復(fù)雜的函數(shù)運(yùn)算化簡為移位和加法的迭代，提供了一種適合硬件的實(shí)現(xiàn)方法[13?14]。其圓周坐標(biāo)下的基本迭代公式如下：

其中每次旋轉(zhuǎn)的基本角取，為每次旋轉(zhuǎn)的方向，是根據(jù)角度累積量來判定的，即，則。CORDIC算法可以工作在旋轉(zhuǎn)模式和矢量模式，旋轉(zhuǎn)模式用來完成坐標(biāo)旋轉(zhuǎn)運(yùn)算，矢量模式用來完成反正切計(jì)算得到最佳旋轉(zhuǎn)角，有。在旋轉(zhuǎn)模式情況下，當(dāng)， 經(jīng)過n次基本角旋轉(zhuǎn)迭代后旋轉(zhuǎn)結(jié)果為：

式中：為對(duì)結(jié)果進(jìn)行修正的比例因子，。通過比較式（12）和式（6）可以很明了地看出CORDIC算法可以解決Jacobi的旋轉(zhuǎn)問題。

4 FPGA實(shí)現(xiàn)方案

由于待分解矩陣是實(shí)對(duì)稱矩陣，因此在FPGA中只需要處理其上三角元素就可以取代整個(gè)矩陣，掃描順序?yàn)椤Ｔ诮?jīng)典Jacobi算法中， 一次Jacobi掃描要進(jìn)行1次CORDIC反正切計(jì)算和2次CORDIC旋轉(zhuǎn)才能完成，單次掃描時(shí)間長，針對(duì)這一不足，具體實(shí)現(xiàn)時(shí)對(duì)其進(jìn)行簡化。

由于Jacobi中，將根據(jù)的正負(fù)性判定，從而確定符號(hào)集，符號(hào)集一旦確定就可以根據(jù)式（11）的迭代實(shí)現(xiàn)CORDIC旋轉(zhuǎn)，省去了反正切計(jì)算的過程。

令，由可以得到：

根據(jù)式（9）可知應(yīng)取初值為，，這樣根據(jù)此更新公式和初值就可以迭代求得符號(hào)集。

將式（6）～式（8）展開可以得到：

則當(dāng)符號(hào)集確定可以將符號(hào)集作為地址，通過查表的方式得到，從而實(shí)現(xiàn)對(duì)角元素式（13）的更新，而對(duì)于要掃描的非對(duì)角元素可以直接置零。式（14）中行（列）行（列）其他元素的更新則要用1次CORDIC旋轉(zhuǎn)來實(shí)現(xiàn)，初值為。這里符號(hào)計(jì)算和CORDIC旋轉(zhuǎn)可以并行進(jìn)行，每算出一個(gè)值就進(jìn)行一個(gè)基本角的旋轉(zhuǎn)。總體結(jié)構(gòu)如圖2所示。

5 設(shè)計(jì)實(shí)現(xiàn)及結(jié)果

本設(shè)計(jì)所選用的器件是ALTERA公司的EP3C120F780C7器件，使用Verilog語言完成電路描述之后，在Quartus Ⅱ軟件平臺(tái)上進(jìn)行編譯，在Modelsim中進(jìn)行仿真。以矩陣為例，采用21 b定點(diǎn)數(shù)，將如下矩陣作為輸入矩陣進(jìn)行系統(tǒng)仿真。

將所得定點(diǎn)數(shù)結(jié)果化為對(duì)應(yīng)的浮點(diǎn)數(shù)數(shù)據(jù)在表1中給出，在Matlab中直接調(diào)用eig函數(shù)計(jì)算出的特征值和特征向量在表2中給出。表1和表2分別表示設(shè)計(jì)結(jié)果與理論結(jié)果。進(jìn)行比較可以看出，本設(shè)計(jì)所得結(jié)果與理論結(jié)果基本一致；雖然第4列和第5列特征向量的符號(hào)相反，但由于各特征向量間的乘積為一個(gè)很小的接近于0的數(shù)，故其正負(fù)不影響其正交性。

假設(shè)有任意兩個(gè)信號(hào)入射到五元均勻圓陣上，入射方位角分別為120°和160°，圖3給出了用MUSIC算法進(jìn)行DOA估計(jì)的結(jié)果，并對(duì)局部進(jìn)行了放大。其中，圖3中（a）為理論估計(jì)結(jié)果，其EVD計(jì)算部分采用Matkab浮點(diǎn)計(jì)算，圖3中（b）為本文設(shè)計(jì)估計(jì)結(jié)果，其EVD計(jì)算部分采用本設(shè)計(jì)的定點(diǎn)操作，除了EVD計(jì)算，所有MUSIC算法中的其他計(jì)算都采用浮點(diǎn)操作。從圖中可以看出，本設(shè)計(jì)可以有效估計(jì)出入射信號(hào)的DOA信息，證明此硬件結(jié)構(gòu)可以應(yīng)用于MUSIC算法的FPGA實(shí)現(xiàn)中。

6 結(jié) 語

本文采用CORDIC算法和Jacobi算法，在FPGA上實(shí)現(xiàn)了數(shù)據(jù)協(xié)方差矩陣的特征值和特征向量計(jì)算，為MUSIC算法的FPGA實(shí)現(xiàn)提供了保證。實(shí)對(duì)稱矩陣特征分解應(yīng)用廣泛，其硬件實(shí)現(xiàn)對(duì)其他算法的工程應(yīng)用有重要作用。

參考文獻(xiàn)

[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] 徐德琛，劉志文，徐友根.某測向系統(tǒng)中MUSIC算法的FPGA實(shí)現(xiàn)[J].北京理工大學(xué)學(xué)報(bào)，2010，30（9）：1107?1111.

[10] 于正文，尹慶莉.求特征值的Jacobi方法[J].山東科學(xué)，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算法的實(shí)現(xiàn)[J].西安郵電學(xué)院學(xué)報(bào)，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.