胡鐵喬，張毛毛，李陽波

（中國民航大學(xué)天津市智能信號與圖像處理重點實驗室，天津 300300）

可配置Cholesky分解矩陣求逆的FPGA實現(xiàn)

胡鐵喬，張毛毛，李陽波

（中國民航大學(xué)天津市智能信號與圖像處理重點實驗室，天津 300300）

在陣列信號抗干擾處理算法中，功率倒置法作為簡單有效的干擾抑制手段被廣泛應(yīng)用，本算法的關(guān)鍵是求出對應(yīng)協(xié)方差矩陣的逆矩陣。鑒于FPGA的速度優(yōu)勢，目前陣列信號抗干擾處理算法多采用FPGA實現(xiàn)。本研究針對常用陣列信號抗干擾系統(tǒng)的實現(xiàn)，設(shè)計出一種基于FPGA的Cholesky矩陣分解模塊，可適用于3～8階正定Hermitian矩陣的求逆，配置簡單方便。測試結(jié)果表明，在不影響計算精度的情況下，該方法節(jié)省了FPGA中的硬件資源。

Cholesky分解；矩陣求逆；陣列信號抗干擾；FPGA實現(xiàn)

目前，陣列天線抗干擾比較常用的算法包括功率倒置法、自適應(yīng)波束形成法和盲自適應(yīng)波束形成法等。由于功率倒置法不需知道信號的波達(dá)方向，只在干擾來向上形成零陷，因此該方法得到了廣泛應(yīng)用[1]。功率倒置法的關(guān)鍵是計算協(xié)方差矩陣的逆矩陣。

通過矩陣分解實現(xiàn)矩陣求逆的運算量小且易于實現(xiàn)，本文采用矩陣分解方法得到協(xié)方差矩陣的逆矩陣。常用的矩陣分解方法有：Cholesky分解、LU分解和QR分解。其中，Cholesky分解適用于共軛對稱正定矩陣；LU分解適用于順序主子式不為0的任意矩陣；QR分解適用于任意矩陣。在功率倒置法的算法中，協(xié)方差矩陣為共軛對稱特性，所以采用Cholesky分解方法求解。

文獻(xiàn)[2-3]采用基于Cholesky分解的矩陣求逆算法，其算法實現(xiàn)中使用較多的乘法器和除法器，這樣大大提高了運算速度。文獻(xiàn)[4]同樣采用基于Cholesky分解的矩陣求逆算法，雖然也是可配置，但配置的矩陣維度較少，只能配置2、3、4三種維度的矩陣，且其算法在計算不同維度矩陣時消耗相同的資源。

目前，常用的陣列信號抗干擾處理系統(tǒng)使用的陣元數(shù)目多為4、5和7陣元，本文針對常用的陣列信號處理算法，對協(xié)方差矩陣求逆提出了基于FPGA的Cholesky分解，其適用于3～8階的協(xié)方差矩陣求逆算法，鑒于協(xié)方差矩陣計算耗時較多，算法對時間的敏感度不強，于是采用分時復(fù)用計算，節(jié)省了硬件資源。

1 基于Cholesky分解求逆方法

根據(jù)Cholesky分解定理可得[4-5]：對角線元素為實數(shù)的對角矩陣D和對角線元素為1的下三角矩陣L及其共軛轉(zhuǎn)置矩陣 LH，A=LDLH，即

其中：A為n×n正定Hermitian矩陣；W為A的逆矩陣；L為上三角矩陣；D為對角矩陣；C和Y均為下三角矩陣。

第1步由矩陣A分解得到矩陣L和矩陣C，令C=LD，即可得

由式（2）根據(jù)矩陣乘法公式可得，矩陣C的第一列與矩陣A的第一列元素相同，即得矩陣C的第一列。然后由矩陣C第一列的元素計算矩陣L的第一列元素，即

然后通過交叉求解按列由A求出矩陣C和矩陣L的其他列元素，具體計算公式為

第2步求矩陣C的逆矩陣Y。由于W為A的逆矩陣，AW=CLHW=E，令 LHW=Y，得 CY=E，且 Y為下三角矩陣，即

由式（6）根據(jù)矩陣乘法公式，按照對角線元素優(yōu)先求取，可得矩陣Y的對角線元素為

然后計算與對角線元素平行元素依次求取，具體計算公式為

第3步由矩陣L和矩陣Y求矩陣W，即

由矩陣性質(zhì)可知，共軛對稱矩陣的逆矩陣也為共軛對稱矩陣，即矩陣W為共軛對稱矩陣，只需求解其一半元素即可獲得矩陣W，從下至上依次計算，具體計算公式為

由以上3步即可完成一個正定Hermitian矩陣的求逆過程[6-7]。

2 可配置矩陣求逆的硬件實現(xiàn)

2.1 總體設(shè)計結(jié)構(gòu)

本文提出的針對3～8階正定Hermitian矩陣求逆的可配置Cholesky分解模塊電路，包括矩陣輸入、輸出模塊，矩陣運算模塊和控制模塊，其中矩陣運算模塊包括矩陣A分解（矩陣C、L交叉計算）、矩陣C的求逆電路（矩陣Y的計算）、逆矩陣W的計算。根據(jù)式（4）、式（7）和式（10）的運算可知，矩陣的求逆過程需進(jìn)行大量的加法、乘法和除法運算，而乘法器和除法器會占用大量的計算資源，所以在實現(xiàn)過程中需在“資源-速度”上做出權(quán)衡，系統(tǒng)采用了分時復(fù)用技術(shù)。硬件結(jié)構(gòu)示意圖如圖1所示（圖中實線箭頭表示數(shù)據(jù)的計算流程，虛線箭頭表示控制信號的流程）。

圖1 硬件實現(xiàn)框圖Fig.1 Hardware implementation block diagram

矩陣A的實部和虛部與使能信號進(jìn)入Cholesky求逆系統(tǒng)，使能信號進(jìn)入控制單元完成計數(shù)器清0和控制信號復(fù)位，矩陣求逆計算由控制單元控制矩陣A分解，交叉求解得到矩陣L和矩陣C；然后將矩陣L緩存，同時計算矩陣C的逆矩陣Y；矩陣Y計算完成后由控制單元控制矩陣Y和矩陣L計算得到逆矩陣W，最后輸出逆矩陣和計算完成信號。

2.2 控制模塊

該模塊的功能是完成對Cholesky分解求逆模塊的運算控制。本模塊信號包括使能信號、輸入矩陣維度、可配置Count_Chol計數(shù)器和控制字。圖2是5×5矩陣的一個除法器的分時復(fù)用控制圖。

圖2 除法器分時復(fù)用控制圖Fig.2 Divider time-sharing control chart

當(dāng)使能信號輸入，計數(shù)器Count_Chol從0開始計數(shù)。當(dāng)計數(shù)器計數(shù)到某一數(shù)值時，產(chǎn)生信號對除法器進(jìn)行賦值。控制字分為3種，分別為矩陣C和矩陣L交叉計算的控制字Scl_1、Scl_2…；矩陣C的逆矩陣Y計算的控制字Sy_1、Sy_2…；以及矩陣W計算的控制字Sw_1、Sw_2…?？刂谱址柕臄?shù)字則代表計算的某一行或某一列，如Scl_1代表開始計算式（3）（即矩陣C的第二列），Sy_1則代表開始計算式（6）（即矩陣Y的對角線元素）。與此相次，復(fù)乘法器也通過相同計數(shù)器進(jìn)行賦值。

以計算5階矩陣的逆矩陣中除法器的賦值為例，如圖2所示，根據(jù)式（3）和式（4）計算矩陣C的各列元素，當(dāng)計數(shù)器 Count=Scl_1+divP（i-2）（divP 表示兩次除法之間的間隔）時，即 Count=[0、5、10、15]時生成使能信號給除法器賦值c21/a11、c31/a11…；當(dāng)Count=Scl_2+divP（i-3）時，即 Count=[20、25、30]時生成使能信號給除法器賦值c32/a22、c42/a22…，依此類推。根據(jù)式（7）計算矩陣Y的對角線元素，當(dāng)Count=Sy_1+divP（i-2）時，即 Count=[75、80、85、90、95]時生成使能信號給除法器賦值1/c11、1/c22…，乘法器的賦值與此類似。

對于不同維度的矩陣，其求逆的計算過程相似，但計算的元素數(shù)目不同，通過改變控制字來控制計算過程。根據(jù)矩陣的階數(shù)來計算不同階矩陣的控制字，且每兩個控制字之間相互聯(lián)系。以5階矩陣為例，介紹其中兩個控制字之間的計算關(guān)系。Scl_3=Scl_2+divP（ChNO-1），其中 ChNO 為矩陣階數(shù)，divP（ChNO-1）則表示計算矩陣C的第三列和矩陣L的第二列所需要的時延。Sy_1=Scl_3+divP（ChNO-2），其中ChNO為矩陣階數(shù)。

不同維度矩陣的控制字?jǐn)?shù)目也不同，隨著矩陣的維度遞減，需要的控制信號數(shù)目也遞減。以8×8矩陣為例，整個計算過程需要21個控制字（Scl_1～Scl_6、Sy_1～Sy_8和 Sw_1～Sw_7），而對于 5× 5矩陣，整個計算過程則只需要12個控制信號（Scl_1～Scl_3、Sy_1～Sy_5和Sw_1～Sw_4），以此類推。這樣做不僅簡單靈活，而且節(jié)省了硬件資源。

2.3 矩陣運算模塊

矩陣計算模塊中采用分時復(fù)用的方法，由于系統(tǒng)處理的信號均為復(fù)信號，所以本系統(tǒng)占用了2個除法器（實部和虛部）和3個復(fù)乘模塊，數(shù)據(jù)通過控制模塊發(fā)出控制信號作為各個乘法器和除法器的控制輸入使能。其中矩陣C、矩陣Y和矩陣W在計算過程中分別復(fù)用了一個復(fù)乘模塊，矩陣L、矩陣C和矩陣Y在計算過程中共同復(fù)用了除法器模塊。

由于FPGA計算浮點數(shù)時資源消耗較大，且計算速度較慢，所以在算法實現(xiàn)的過程中并保證精度的情況下采用移位為定點數(shù)計算。除法器和乘法器的運算結(jié)果均移位為整數(shù)。對于不同的系統(tǒng)要求，可通過修改移位位寬來滿足不同的精度要求。

在算法實現(xiàn)過程中的復(fù)乘模塊，要實現(xiàn)兩個復(fù)數(shù)的共軛相乘，若用常規(guī)方法計算，根據(jù)

可知，該模塊需占用4個乘法器和2個加法器。本文所用計算法為

采用修改后的方法，雖然加法器的數(shù)目增加了4個，但是乘法器只用了3個，這樣雖增加了邏輯資源占用，但是節(jié)省了更多的運算資源。

對于不同階矩陣，采用相同的計算模塊。通過不同的控制字控制運算模塊的賦值，這樣可做到矩陣階數(shù)越低，用時越短。若計算8×8矩陣的逆矩陣，需620個時鐘周期，而如果計算3×3矩陣，則只需230個時鐘周期。

3 資源占比及實驗結(jié)果分析

完成算法設(shè)計后，驗證算法的可行性，本文所使用的硬件編程語言為Verilog HDL，硬件編程環(huán)境為Xilinx ISE 13.4軟件。并在Chip Scope Pro Analyzer中驗證運算結(jié)果。

測試數(shù)據(jù)為7×7協(xié)方差矩陣，輸入矩陣A數(shù)據(jù)位寬為32位，輸出矩陣W數(shù)據(jù)位寬為20位。

3.1 資源占比

由于本方法采用的FPGA芯片為Xilinx Virtex-5 XC5VSX95T，此芯片含有58 880個片寄存單元（Slice Registers）和 58 880 個可編程邏輯單元（Slice LUTs），以及640個DSP核。表1為不同維度矩陣求逆所占用的資源對比。

表1 硬件資源占用情況對比（占用數(shù)目）Tab.1 Comparison of coccupied hardware resource

從表1中可以看出，計算的矩陣維度越少，系統(tǒng)占用的片寄存單元和可編程邏輯單元越少。

3.2 結(jié)果驗證

結(jié)果驗證方法采用FPGA硬件運算的結(jié)果與Matlab調(diào)用函數(shù)求逆的結(jié)果相比較。為了測試對比，inv_w為直接調(diào)用Matlab inv函數(shù)的結(jié)果，并擴(kuò)大230倍的取整值，即

硬件計算結(jié)果采用按行輸出，結(jié)果如圖3所示，圖中wr（1～7）_mux表示逆矩陣每一行的實部，wi（1～7）_mux表示逆矩陣每一行的虛部，Count_Chol即為控制計數(shù)器。對比inv_w和圖3可以看出，硬件實現(xiàn)結(jié)果與Matlab計算結(jié)果相吻合。

圖3 Cholesky矩陣求逆計算結(jié)果Fig.3 Calculation result of Cholesky matrix inversion

4 結(jié)語

本文基于陣列天線抗干擾功率倒置法的實現(xiàn)，綜合考慮常用陣元數(shù)目，提出了一種基于Cholesky的協(xié)方差矩陣求逆算法，可靈活配置3～8階協(xié)方差矩陣，降低了矩陣求逆的實現(xiàn)難度。在FPGA實現(xiàn)過程中，采用分時復(fù)用乘法器和除法器，節(jié)省了FPGA的計算資源。FPGA計算結(jié)果與Matlab直接調(diào)用函數(shù)計算結(jié)果相一致。另外，更多維度的矩陣求逆也可參照本算法實現(xiàn)。

[1]吳仁彪,盧 丹,李 嬋.用于高動態(tài)GPS接收機(jī)的微分約束最小功率抗干擾方法[J].中國科學(xué),2011,41(8):968-977.

[2]CHU Xuezheng,KHALED B,JOHN T.Rapid Prototyping of an Improved Cholesky Decomposition Based MIMO Detector on FPGAs[R].AHS,NASA/ESA Conference on,2009.

[3]魏嬋娟,張春水,劉 健.一種基于Cholesky分解的快速矩陣求逆方法設(shè)計[J].電子設(shè)計工程,2014,22(1):159-164.

[4]潘 曉,許有云,甘小鶯.基于Cholesky分解的可配置矩陣求逆FPGA實現(xiàn)[J].信息技術(shù),2009,5(11):141-143.

[5]張賢達(dá).矩陣分析與應(yīng)用[M].北京:清華大學(xué)出版社,2004:225-227.

[6]李 濤,張忠培.矩陣求逆的FPGA實現(xiàn)[J].通信技術(shù),2010,43(11):147-149.

[7]郭 磊.矩陣運算的硬件加速技術(shù)研究[D].長沙:國防科學(xué)技術(shù)大學(xué),2010.

（責(zé)任編輯：楊媛媛）

FPGA implementation of matrix inversion based on configurable Cholesky decomposition

HU Tieqiao,ZHANG Maomao,LI Yangbo
（Intelligent Signal and Image Processing Key Lab of Tianjin,CAUC,Tianjin 300300,China）

As a simple and effective interference suppression algorithm,power inversion algorithm is widely used in various arraysignalprocessingalgorithms.Thekeyofthealgorithmistofindtheinversematrixofcorrespondingcovariance matrix.Given the high speed of FPGA,array signal anti-jamming processing algorithm is mostly accomplished by FPGA at present.In order to realize the commonly used system of array signal anti-jamming system,a Cholesky matrix decomposition based on FPGA is designed,which can be applied to the inversion of 3～8 positive definite Hermitian matrixes.Test results show that this algorithm can save hardware resource without affecting calculation precision.

Cholesky decomposition;matrix inversion;array signal anti-jamming;FPGA implementation

TN919；V243.1

：A

：1674－5590（2017）04－0007－04

2017-02-22；

：2017-04-09

：國家重點研發(fā)計劃（2016YFB0502402）

胡鐵喬（1970—），男，河南洛陽人，副教授，工學(xué)碩士，研究方向為基于FPGA、DSP高速處理平臺的陣列信號處理系統(tǒng).