魏嬋娟，張春水，劉 健

（中國空間技術(shù)研究院航天恒星科技有限公司 北京 100086）

抗干擾接收機由于增加了抗干擾處理的部分將會對導航信號的實時解算造成一定的影響，求取抗干擾濾波權(quán)值用時越長，對解算的延遲越大，造成的定位誤差也將越大。在權(quán)值更新以及抗干擾濾波的整個處理過程中，矩陣求逆的用時占用了90%以上的時間，矩陣求逆用時過長是導致接收機權(quán)值更新過慢的最主要因素，因此若采用直接矩陣求逆對抗干擾算法進行實現(xiàn)，必須保證矩陣求逆的用時極短。本文針對大維數(shù)赫米特矩陣求逆實現(xiàn)方法進行FPGA設計及實現(xiàn)。

1 Cholesky分解求逆方法及其流水設計

A為正定的赫米特矩陣，采用Cholesky分解矩陣求逆方法的實現(xiàn)過程為[1]:

1）由Cholesky分解得到下三角矩陣L及其共軛轉(zhuǎn)置矩陣LH，A =L×LH；

2） 對下三角矩陣L進行矩陣求逆計算，得到其逆矩陣IL=L-1；

3） 由IL矩 陣 求 得A的 逆 矩 陣IA，IA =（LH）-1×L-1=LH×IL。

在FPGA實現(xiàn)赫米特矩陣求逆的關(guān)鍵在于實現(xiàn)方法的流水設計，合理高效的流水方法可以最大程度地節(jié)約資源提高實現(xiàn)速度，以下分3部分實現(xiàn)赫米特矩陣Cholesky分解求逆的實現(xiàn)方法設計[2-4]。

1.1 矩陣的Cholesky分解

對A矩陣進行Cholesky分解，得到下三角矩陣L，A=LLH。采用分塊的方法計算L矩陣[5]。

將矩陣L及A進行分塊表示如下：

由A=LLH推導可得：

根據(jù)如上公式表示，L的第一列的值可由A的第一列的值獲得，其計算方法為：

在計算完L的第一列的之后，之后需要計算L1，根據(jù)式2可得:

在將A矩陣的部分更新為A1-llH后，可采用相同的方法計算出L1矩陣的第一列的值，即L矩陣的第2列的值，依次類推可計算出L的所有值。

對以上的Cholesky矩陣分解方法進行FPGA實現(xiàn)的流水設計，以4×4矩陣為例，可分為四步分析各個數(shù)據(jù)之間的依賴關(guān)系并計算L各列的值。

圖1 矩陣Cholisky分解步驟圖Fig.1 Four steps of Cholisky decomposition

步驟1：計算L矩陣的第1列；之后更新A矩陣（圖中，左側(cè)的矩陣為A矩陣，右側(cè)為L矩陣，箭頭所指向的方向為計算當前值所依賴的數(shù)據(jù)，例如將頭由a指向b，表示計算a時需要用到b，需在b計算完之后才可進行a的計算）。

步驟2：計算L矩陣的第2列；之后更新A矩陣。

步驟3：計算L矩陣的第3列；之后更新A矩陣。

步驟4：計算L矩陣的第4列；A矩陣無需再更新。

對于4×4矩陣，在經(jīng)過4步之后可獲得L矩陣。

1.2 下三角矩陣求逆

對下三角矩陣L進行求逆計算得到的逆矩陣IL[6]，根據(jù)L*IL=E，分別求得IL各列的值。

L矩陣求逆數(shù)據(jù)依賴關(guān)系及流水設計如圖2所示。

步驟1：計算IL對角線上的值（圖中，左側(cè)的矩陣為L矩陣，右側(cè)的矩陣為IL矩陣）。步驟2：計算IL第二對角線上的值。步驟3：計算IL第三對角線上的值。步驟4：計算IL第四對角線上的值。

通過對對角線上的值進行同步計算，在經(jīng)過4個步驟之后可計算出IL矩陣的所有值。

1.3 矩陣相乘

上三角矩陣LH與下三角矩陣L相乘得到A的逆矩陣IA，IA=ILHIL。

矩陣相乘流水設計數(shù)據(jù)間的依賴關(guān)系（圖中左側(cè)矩陣為L矩陣，右側(cè)矩陣為IA矩陣）。

矩陣相乘較為簡單，其各個數(shù)據(jù)間的計算無需等待，僅僅依賴于已知的下三角矩陣IL。

圖2 下三角矩陣求逆數(shù)據(jù)依賴關(guān)系圖Fig.2 Data dependency relationship of lower triangular matrix inversion

圖3 矩陣相乘的數(shù)據(jù)依賴關(guān)系圖Fig.3 Data dependency relationship of matrix vector multiplication

矩陣相乘較為簡單，其各個數(shù)據(jù)間的計算無需等待，僅僅依賴于已知的下三角矩陣IL。

2 FPGA實現(xiàn)

2.1 狀態(tài)控制及功能描述

與1節(jié)Cholesky分解求逆流水方法對應，F(xiàn)PGA硬件實現(xiàn)電路包括3個主要模塊：Cholesky分解、下三角矩陣求逆和矩陣相乘。其在FPGA中實現(xiàn)的狀態(tài)控制及存儲如下圖所示：

圖4 FPGA實現(xiàn)Cholesky矩陣分解求逆狀態(tài)控制圖Fig.4 FPGA implementation of Cholesky decomposition and inversion

功能描述：

1）Start_inv為矩陣分解求逆開始標志，save狀態(tài)下，將生成的協(xié)方差矩陣存入RAM_A；

2）Chol狀態(tài)下完成協(xié)方差矩陣的Cholesky分解，得到的下三角矩陣L存入RAM_B；

3）Inv狀態(tài)下完成下三角矩陣的求逆計算，將得到的逆矩陣IL存回RAM_A；

4）Inv_mult狀態(tài)下完成ILH矩陣與IL矩陣的相乘，即上三角陣與下三角陣的相乘，其結(jié)果存入RAM_B。

5）在Inv_mult狀態(tài)結(jié)束后進入Output狀態(tài)將RAM_B中的逆矩陣輸出。

6 ）在Ouptput狀態(tài)結(jié)束后重新返回start狀態(tài)，開始對下一組協(xié)方差矩陣的分解求逆。

2.2 用時測試

跟據(jù)以上的流水設計進行FPGA代碼的編寫實現(xiàn)，矩陣求逆的整體用時如下表所示：

表1 FPGA實現(xiàn)Cholesky矩陣求逆用時Tab.1 FPGA implementation time table

由上表看出，28×28維矩陣的求逆用時小于0.4 ms，35×35維矩陣的求逆用時僅需不到0.8 ms。

3 結(jié) 論

文中針對大維數(shù)赫米特矩陣的求逆的問題展開，介紹了適用于正定赫米特矩陣的Cholesky矩陣分解求逆方法，并對該方法進行了高效的流水設計，分3步實現(xiàn)矩陣的FPGA分解求逆，經(jīng)過試驗測試，用時極短。

[1]張賢達.矩陣分析與應用[M].北京:清華大學出版社, 2004.

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

[3]Thompson, J. Benkrid, K. Xuezheng Chu. Rapid Prototyping of an Improved Cholesky Decomposition Based MIMO Detector on FPGAs [J].Adaptive Hardware and Systems, 2009(4):369-375.

[4]Mach D M,Koshak W J.General matrix inversion for the calibration of electric field sensor arrays on aircraft platforms [J].Marshall Space Flight Center, 2006(6):1576-1587.

[5]魯翠仙.分塊矩陣在矩陣求逆中的應用[D].昆明:云南大學,2009.

[6]彭玲,彭大芹. 一種下三角復矩陣求逆方法的IP設計與實現(xiàn)[J].電子測試, 2011(10):9-12.

PENG Ling, PENG Da-qin.IP design and implementarion of lower trlangular matrx invertion[J].Electronic Test,2011(10):9-12.