基于現(xiàn)場可編程門陣列的矩陣求逆算法設(shè)計

安國臣, 劉若凡, 趙滿, 袁玉鑫, 王曉君

(1.河北科技大學(xué)信息科學(xué)與工程學(xué)院, 石家莊 050018; 2.河北晶禾電子技術(shù)股份有限公司, 石家莊 050200)

近年來,北斗衛(wèi)星導(dǎo)航在人們的日常生活和軍事應(yīng)用領(lǐng)域中的地位越來越重要[1]。但是由于北斗導(dǎo)航衛(wèi)星在軌高度大約為21 500 km,導(dǎo)致衛(wèi)星信號抵達地面時的功率非常低,所以北斗導(dǎo)航信號很容易被外界干擾,進而失去定位能力[2-4]。因此,北斗導(dǎo)航抗干擾技術(shù)的研究與實現(xiàn)已成為相關(guān)領(lǐng)域的重大研究方向。

采用陣列信號處理方法來進行衛(wèi)星導(dǎo)航系統(tǒng)的抗干擾算法有很多,其中采樣矩陣求逆(sample matrix invert,SMI)算法具有很好的抗干擾性能,得到了廣泛的應(yīng)用。然而SMI算法需要進行矩陣求逆運算導(dǎo)致運算量過大,尤其是采用空時聯(lián)合抗干擾時,矩陣維數(shù)很大,導(dǎo)致運算量激增,很難滿足動態(tài)性能要求[5-6]。文獻[7-8]中均對SMI算法的性能進行分析,說明了該算法相比較于其他的自適應(yīng)處理算法有較好抗干擾性能,但是由于運算量較大導(dǎo)致硬件難以實現(xiàn)。

為了解決上面所說的運算量大的問題,常采用矩陣分解的方法進行求逆,其中矩陣分解常用的方法有QR分解、Cholesky分解等方法,QR分解需要將矩陣分解成一個正交矩陣和一個上三角矩陣,文獻[9]通過采用QR分解基于Vivado 擁有的高層次綜合(high level synthesis,HLS)開發(fā)工具,完成了64階矩陣求逆設(shè)計,且時鐘周期為4 ns時,矩陣求逆在2.3～9.2 ms完成,但矩陣求逆運算量較大,計算時間相對較長。而基于Cholesky分解的矩陣求逆算法,是將矩陣的求逆分解為三角矩陣的求逆,大大降低了矩陣求逆的運算量[10-13]。文獻[14]提出了一種基于現(xiàn)場可編程門陣列(field programmable gate array,FPGA)的Cholesky分解矩陣求逆實現(xiàn)方法,其求逆采用分時復(fù)用的方法進行實現(xiàn),節(jié)省了硬件資源,但該方法僅適用于3～8階的協(xié)方差矩陣求逆;文獻[15]基于FPGA首先利用Cholesky分解求解出下三角矩陣,然后使用Gauss消元法求解權(quán)值的方法,該文獻中提到的對4階協(xié)方差矩陣進行求解權(quán)值,所消耗時間為16 μs。文獻[16]利用Vivado HLS開發(fā)工具,基于C/C++語言實現(xiàn)了Cholesky分解矩陣求逆算法,該文獻對48階協(xié)方差矩陣進行求逆,僅求逆模塊所消耗時間為0.9 ms,并未對協(xié)方差矩陣計算模塊進行統(tǒng)計說明。

由于SMI抗干擾算法權(quán)值更新速率主要取決于逆矩陣的求解速度,而權(quán)值更新速率又對動態(tài)抗干擾性能有著重大影響。為了提高權(quán)值更新率,縮短逆矩陣求解延遲,現(xiàn)提出一種基于FPGA的快速采樣矩陣求逆算法實現(xiàn)架構(gòu)。該架構(gòu)采用流水線加狀態(tài)機的實現(xiàn)結(jié)構(gòu),該結(jié)構(gòu)在有限的FPGA資源下,實現(xiàn)計算速度最大化,實現(xiàn)運算速度與資源占用的平衡。根據(jù)項目需求,實現(xiàn)7陣元8抽頭即56階矩陣求逆算法的設(shè)計與仿真驗證。

1 Cholesky矩陣求逆算法原理

Cholesky矩陣求逆算法可分為三個步驟:計算Cholesky分解矩陣L、計算下三角矩陣L的逆矩陣S和三角矩陣相乘。下面對每個步驟進行具體說明。

1.1 計算Cholesky分解矩陣

設(shè)協(xié)方差矩陣R由一個下三角矩陣L與其共軛轉(zhuǎn)置LH的乘積組成,即

R=LLH

(1)

式(1)中:

(2)

(3)

(4)

式(4)中:*表示對數(shù)據(jù)取共軛;M為采樣矩陣階數(shù)。

通過R=LLH式關(guān)系,可以得到R中元素與LLH中元素的對應(yīng)關(guān)系為

(5)

(6)

(7)

式(5)表示第一列元素;式(6)表示除第一個元素外的對角線元素(i=j≠1);式(7)表示除第一列和對角線元素外的下三角元素(i>j,i≠1,j≠1)。

根據(jù)上面的對應(yīng)關(guān)系,可以反推L中的元素,將其分為4塊:第1塊為第1個元素l11;第2塊為第1列除第1個元素外的元素li1;第3塊為上三角陣中的0;第4塊為除第1列元素和0外的元素LM-1;如式(8)所示。

(8)

那么可以得到各個分塊的遞推式,其中第一列元素為

(9)

(10)

LM-1要分成兩個部分求取,第一部分為對角線上元素為

(11)

第二部分為非對角線上元素,即

(12)

至此,就完成了對下三角矩陣L的求取,即完成了對R矩陣的Cholesky分解。

1.2 計算下三角矩陣L的逆矩陣

現(xiàn)在L已知,接下來需要計算其逆矩陣L-1,假設(shè)有矩陣S=L-1,則滿足

(13)

式(13)中:IM為M階單位矩陣,分析式(13)可知,S也為下三角矩陣,即

(14)

重新整理式(13)可得等式兩側(cè)矩陣中元素的對應(yīng)關(guān)系,將其分為3塊求取,第1塊為對角線元素(i=j);第2塊為除對角線外的下三角矩陣中元素(i>j,i≠1);第3塊為上三角0元素。

那么可得各個分塊的遞推式,其中第一塊中元素為

(15)

第二塊中的元素為

(16)

至此,就完成了對下三角矩陣L的求逆,即得到了L矩陣的逆矩陣S。

1.3 三角矩陣相乘

計算出S后,根據(jù)R-1=(L-1)HL-1,可以得到R-1的表達式為

R-1=SHS

(17)

式(17)中:S為M×M維下三角矩陣。

1.4 計算權(quán)值

基于線性約束最小方差(linear constrained minimum variance, LCMV)準則的功率倒置(powerinversion,PI)算法將第一陣元的增益約束為1,這樣可避免加權(quán)矢量出現(xiàn)零解的情況[17-19],其算法下的權(quán)值為

Wopt=R-1C(CHR-1C)-1

(18)

式(18)中:C=[1,0,…,0]T為M×1維矢量。

2 FPGA實現(xiàn)架構(gòu)

針對上述的運算步驟,研究各環(huán)節(jié)算法特點,綜合考慮計算速度與資源占用情況的權(quán)衡,分別設(shè)計了5個運算模塊:協(xié)方差矩陣計算模塊、Cho-lesky分解模塊、下三角矩陣求逆模塊、三角矩陣相乘模塊和計算權(quán)值模塊,系統(tǒng)總體實現(xiàn)架構(gòu)如圖1所示。

圖1 系統(tǒng)總體實現(xiàn)架構(gòu)Fig.1 Overall implementation architecture of the system

系統(tǒng)實現(xiàn)架構(gòu)采用流水線加狀態(tài)機結(jié)構(gòu),流水線結(jié)構(gòu)可減少運算時間,提高計算效率,狀態(tài)機結(jié)構(gòu)用來實現(xiàn)系統(tǒng)中的算法流程控制。

2.1 協(xié)方差矩陣計算模塊設(shè)計

協(xié)方差矩陣計算模塊的實現(xiàn)框圖如圖2所示。從圖2可知,通道數(shù)據(jù)存儲模塊負責(zé)對通道數(shù)據(jù)進行存儲;通道數(shù)據(jù)讀取模塊負責(zé)對通道數(shù)據(jù)進行讀取并延遲形成空時矩陣;動態(tài)連接矩陣模塊中的切換模塊負責(zé)將空時矩陣中的數(shù)據(jù)進行切換后連接至57路復(fù)數(shù)乘累加模塊,為了減少計算時間,57路復(fù)數(shù)乘累加模塊采用了并行執(zhí)行結(jié)構(gòu),即一個時鐘可完成57次復(fù)數(shù)乘累加運算;存儲模塊負責(zé)對計算出來的協(xié)方差矩陣數(shù)據(jù)進行存儲。

AD1、AD2、AD7為通道數(shù)據(jù)信號;AD1_0、AD1_7分別為第1通道第一、第八抽頭數(shù)據(jù)信號;AD7_0、AD7_7分別為第7通道第一、第八抽頭數(shù)據(jù)信號;Rxx1、Rxx57為計算出的協(xié)方差矩陣數(shù)據(jù)信號圖2 協(xié)方差矩陣計算模塊結(jié)構(gòu)框圖Fig.2 Block diagram of covariance matrix calculation module structure

計算協(xié)方差矩陣狀態(tài)機控制模塊流程圖如圖3所示。從圖3可知,首先當(dāng)狀態(tài)機接收到S1_Start信號時,狀態(tài)機由S0狀態(tài)轉(zhuǎn)換到S1狀態(tài),完成通道數(shù)據(jù)的存儲;存儲完成之后,S1_Done信號拉高,狀態(tài)機轉(zhuǎn)換到S2狀態(tài),讀取通道數(shù)據(jù)并形成空時矩陣;之后S2_Done信號拉高,狀態(tài)機轉(zhuǎn)換到S3狀態(tài),此時ctr_line信號從0開始增加1,通過該信號的數(shù)值切換進入57路復(fù)數(shù)乘累加模塊的數(shù)據(jù)用來計算出一組協(xié)方差矩陣結(jié)果,當(dāng)ctr_line信號大于28時轉(zhuǎn)換到空閑狀態(tài),否則轉(zhuǎn)換回S2狀態(tài),進行下一次迭代過程。

圖3 計算協(xié)方差矩陣狀態(tài)機控制模塊流程圖Fig.3 Flow chart of state machine control module for calculating covariance matrix

計算協(xié)方差矩陣中數(shù)據(jù)的規(guī)律圖如圖4所示。其中紅線內(nèi)代表第一次計算的一組對角線數(shù)據(jù),綠線內(nèi)代表第二次計算的一組對角線數(shù)據(jù),按照此規(guī)律,一共切換28次后,協(xié)方差矩陣中的下三角矩陣計算并存儲完成,其存儲方式使用RAM進行存儲,首先定義RAM讀寫地址為12位,將RAM地址的高六位代表行地址,低六位代表列地址,然后存儲時依據(jù)計算出來的數(shù)據(jù)在協(xié)方差矩陣中實際位置來確定行地址和列地址,方便后續(xù)讀取,例如,在圖4中,r11為協(xié)方差矩陣中第1行第1列的數(shù)據(jù),所以在進行存儲時,定義它的行地址為1,列地址為1,即RAM地址為65(000001_000001)。本文研究中后續(xù)數(shù)據(jù)存儲模塊都是采用此規(guī)律進行數(shù)據(jù)存儲。

圖4 計算協(xié)方差矩陣數(shù)據(jù)規(guī)律框圖Fig.4 Block diagram of data law for calculating covariance matrix

2.2 Cholesky分解矩陣求逆模塊

2.2.1 Cholesky分解矩陣

Cholesky分解模塊的結(jié)構(gòu)框圖如圖5所示。從圖5可知,存儲模塊負責(zé)存儲更新后的協(xié)方差矩陣數(shù)據(jù),即主要為更新協(xié)方差矩陣模塊輸出數(shù)據(jù);計算L矩陣是通過浮點開平方模塊和浮點除法模塊,按列依次計算出L矩陣數(shù)據(jù),然后將數(shù)據(jù)傳輸給L矩陣存儲模塊;更新協(xié)方差矩陣模塊中通過浮點復(fù)數(shù)乘累加模塊和浮點減法模塊,按列依次計算出式(11)中根號里以及式(12)中的分子數(shù)據(jù)。

W_Rxx_i_1、W_Rxx_i_2為協(xié)方差矩陣實部數(shù)據(jù)寫信號;W_Rxx_q_1、W_Rxx_q_2為協(xié)方差矩陣虛部數(shù)據(jù)寫信號;R_Rxx_i、R_Rxx_q分別為協(xié)方差矩陣實部、虛部數(shù)據(jù)信號;L_i、L_q分別為計算出的L矩陣實部、虛部數(shù)據(jù)信號;R_L_i、R_L_q分別為L矩陣實部、虛部數(shù)據(jù)讀信號圖5 Cholesky分解矩陣模塊結(jié)構(gòu)框圖Fig.5 Block diagram of Choleskydecomposition matrix module structure

Cholesky分解矩陣狀態(tài)機控制模塊流程圖如圖6所示。從圖6可知,首先當(dāng)狀態(tài)機接收到L_S1_Start信號時,狀態(tài)機由S0狀態(tài)轉(zhuǎn)換到S1狀態(tài),完成協(xié)方差矩陣數(shù)據(jù)的存儲;存儲完成之后L_S1_Done信號拉高,狀態(tài)機轉(zhuǎn)換到S2狀態(tài),計算出一列L矩陣數(shù)據(jù);之后L_S2_Done信號拉高,狀態(tài)機轉(zhuǎn)換到S3狀態(tài),從第二列開始,計算出一列新的協(xié)方差矩陣數(shù)據(jù);之后L_S3_Done信號拉高,狀態(tài)機轉(zhuǎn)換到S4狀態(tài),此時Order信號從0開始增加1,當(dāng)Order>55時,轉(zhuǎn)換到空閑狀態(tài),否則轉(zhuǎn)換到S2狀態(tài)。

圖6 Cholesky分解矩陣狀態(tài)機控制模塊流程圖Fig.6 Flow chart of state machine control module for Cholesky decomposition matrix

2.2.2 計算下三角矩陣L的逆矩陣

計算下三角矩陣L的逆矩陣模塊框圖如圖7所示。從圖7可知,L矩陣存儲模塊負責(zé)存儲Cho-lesky分解矩陣數(shù)據(jù);計算S矩陣模塊中通過復(fù)數(shù)乘累加模塊和浮點除法模塊計算出一組S矩陣數(shù)據(jù)后連接至S矩陣存儲模塊1;S矩陣存儲模塊2負責(zé)存儲S矩陣數(shù)據(jù)。其中S矩陣存儲模塊1存儲的數(shù)據(jù)用于S矩陣的計算,S矩陣存儲模塊2存儲更新后的S矩陣數(shù)據(jù)。

R_S_i_1、R_S_i_2為S矩陣實部數(shù)據(jù)讀信號;R_S_q_1、R_S_q_2為S矩陣虛部數(shù)據(jù)讀信號;R_Lii_i為L矩陣對角線實部數(shù)據(jù)讀信號;S_i、S_q分別為計算出的S矩陣實部、虛部數(shù)據(jù)信號圖7 計算下三角矩陣L的逆矩陣模塊結(jié)構(gòu)框圖Fig.7 Block diagram of calculating the inverse matrix module of the lower triangular matrix L

下三角矩陣L的逆矩陣狀態(tài)機控制模塊流程圖如圖8所示。從圖8可知,首先當(dāng)狀態(tài)機接收到S_S1_Start信號時,狀態(tài)機由S0狀態(tài)轉(zhuǎn)換到S1狀態(tài),完成 Cholesky分解矩陣數(shù)據(jù)的存儲;存儲完成后S_S1_Done信號拉高,狀態(tài)機轉(zhuǎn)換到S2狀態(tài),計算出一組S矩陣數(shù)據(jù);之后狀態(tài)機轉(zhuǎn)換到S3狀態(tài)并返回給狀態(tài)機一個結(jié)束信號S_S2_Done,此時cycle信號從0開始增加1,當(dāng)cycle=55并且S_S2_D-one為高電平時,轉(zhuǎn)換到空閑狀態(tài),否則S_S3_Done信號拉高,轉(zhuǎn)換到S2狀態(tài)。

圖8 計算下三角矩陣L的逆矩陣狀態(tài)機控制模塊流程圖Fig.8 Flow chart of the state machine control module for calculating the inverse matrix of the lower triangular matrix L

計算下三角矩陣的逆矩陣中數(shù)據(jù)的規(guī)律圖如圖9所示。其中紅線內(nèi)代表第一次計算的一組對角線數(shù)據(jù),綠線內(nèi)代表第二次計算的一組對角線數(shù)據(jù),按照此規(guī)律,一共循環(huán)55次后,下三角矩陣L的逆矩陣就全部計算完成。

圖9 計算下三角矩陣的逆矩陣數(shù)據(jù)規(guī)律框圖Fig.9 Block diagram of the inverse matrix data law of the calculation of the lower triangular matrix

2.2.3 三角矩陣相乘及權(quán)值計算

三角矩陣相乘的實現(xiàn)框圖如圖10所示。從圖10可知,三角矩陣相乘模塊中讀取存儲S矩陣中第一列數(shù)據(jù)模塊和S矩陣中的所有數(shù)據(jù)后,連接至復(fù)數(shù)乘累加模塊。

Rxx_inv_i、Rxx_inv_q分別為三角矩陣相乘計算出的實部、虛部數(shù)據(jù)信號圖10 三角矩陣相乘模塊結(jié)構(gòu)框圖Fig.10 Block diagram of triangular matrix multiplication module structure

在計算權(quán)值模塊中,根據(jù)式(18)推導(dǎo)出的計算權(quán)值方法,調(diào)用浮點數(shù)除法核,用來計算出56個權(quán)值數(shù)據(jù)。

3 性能仿真

3.1 數(shù)據(jù)對比

使用同一組通道數(shù)據(jù),分別在MATLAB和FPGA上進行權(quán)值計算,計算出56個權(quán)值數(shù)據(jù)。如表1和表2所示,從表1和表2可以看出,本文所設(shè)計的實現(xiàn)架構(gòu)計算出的權(quán)值數(shù)據(jù)和MATLAB理論仿真的相位誤差很小,滿足進行加權(quán)時對權(quán)值精度的要求。

表1 權(quán)值實部數(shù)據(jù)比較Table1 Comparison of weights real part data

表2 權(quán)值虛部數(shù)據(jù)比較Table 2 Comparison of weights imaginary part data

3.2 權(quán)值更新時間

在100 MHz時鐘下,計算協(xié)方差矩陣模塊消耗約30 000個時鐘,Cholesky分解矩陣求逆模塊消耗約90 000個時鐘,所以計算出一次權(quán)值所消耗約120 000個時鐘,即所消耗時間約為1.2 ms。

在FPGA開發(fā)平臺上進行了系統(tǒng)的整體仿真,以驗證設(shè)計的正確性和時效性,權(quán)值更新時間仿真結(jié)果如圖11所示,可以看出每次權(quán)值更新時間都約為1.2 ms。

圖11 權(quán)值更新時間仿真圖Fig.11 Simulation graph of weight update time

3.3 外場測試

采用以上矩陣求逆算法實現(xiàn)架構(gòu),進行了衛(wèi)星導(dǎo)航信號空時抗干擾系統(tǒng)的整體設(shè)計與實現(xiàn),并進行了外場實測。當(dāng)存在干信比為95 dB的單個白噪聲干擾時,接收機收星效果如圖12所示。由圖12可知,在該強度干擾信號下,系統(tǒng)完成了干擾信號的濾波,并實現(xiàn)了很好的收星效果,定位結(jié)果準確穩(wěn)定,同時也驗證了本文提出的矩陣求逆算法實現(xiàn)架構(gòu)的有效性。

4 結(jié)論

針對自適應(yīng)抗干擾算法在更新最優(yōu)權(quán)值時存在時間延時問題,提出了基于采樣矩陣求逆的自適應(yīng)抗干擾算法的實現(xiàn)架構(gòu),可實現(xiàn)高階協(xié)方差矩陣的求逆運算,能夠有效提高權(quán)值更新速率。該方法已經(jīng)在實際抗干擾系統(tǒng)中得到應(yīng)用,在該硬件平臺上達到了抗95 dB單白噪聲干擾的效果,滿足動態(tài)環(huán)境下權(quán)值更新率的要求。