黃家俊，王建新

（南京理工大學(xué)電子工程與光電技術(shù)學(xué)院，江蘇南京 210094）

陣列信號(hào)處理中的一類(lèi)典型問(wèn)題是如何從接收到的多路觀測(cè)數(shù)據(jù)中解析出被混疊著的原始信號(hào)[1]。常用的陣列信號(hào)分離算法有波達(dá)方向（Direction of Arrival，DOA）估計(jì)[2-4]、盲源分離[5]、快速獨(dú)立成分分析（Fast Independent Component Analysis，F(xiàn)astICA）[6-9]等。這些算法都需要涉及大量的矩陣運(yùn)算，而矩陣特征值分解就是其中的關(guān)鍵步驟之一。

常規(guī)的特征值分解算法主要針對(duì)實(shí)對(duì)稱(chēng)矩陣，如Jacobi 旋轉(zhuǎn)法、QR 迭代法[10]，其中Jacobi 法計(jì)算簡(jiǎn)單、精度高，應(yīng)用更為廣泛。傳感器陣列觀測(cè)到的信號(hào)經(jīng)數(shù)字正交下變頻后變成I/Q 兩路的復(fù)信號(hào)，其協(xié)方差矩陣是一個(gè)復(fù)Hermitian 矩陣，無(wú)法直接應(yīng)用傳統(tǒng)的算法。所以需要通過(guò)一定的運(yùn)算先將復(fù)Hermitian 矩陣轉(zhuǎn)化為實(shí)對(duì)稱(chēng)矩陣。

FPGA 硬件相較于DSP 硬件，運(yùn)行速度快，并行運(yùn)算能力強(qiáng)，更適合解決陣列信號(hào)處理中的矩陣運(yùn)算。因此，該文研究陣列信號(hào)處理中典型的復(fù)Hermitian 矩陣特征值分解問(wèn)題，在FPGA 中實(shí)現(xiàn)矩陣轉(zhuǎn)化和特征值分解的運(yùn)算。

1 矩陣轉(zhuǎn)化

假設(shè)由M個(gè)天線組成的天線陣列觀測(cè)到的信號(hào)經(jīng)正交下變頻后得到的復(fù)觀測(cè)數(shù)據(jù)為x(t)，如式（1）所示。

其中，N為x(t)的采樣長(zhǎng)度。顯然矩陣Rx是一個(gè)M×M階的復(fù)Hermitian 矩陣。

復(fù)Hermitian 矩陣向?qū)崒?duì)稱(chēng)矩陣的轉(zhuǎn)化有兩種處理方法。第一種轉(zhuǎn)化方法是將M×M階的復(fù)Hermitian 矩陣轉(zhuǎn)換為2M×2M階的實(shí)對(duì)稱(chēng)矩陣[11]，首先用實(shí)對(duì)稱(chēng)矩陣Rr和反對(duì)稱(chēng)實(shí)矩陣Ri作為實(shí)部和虛部來(lái)表示矩陣Rx。假設(shè)矩陣Rx的特征值σ對(duì)應(yīng)的特征向量實(shí)部為u、虛部為v，可以得到下式：

易證式（3）中的分塊矩陣是實(shí)對(duì)稱(chēng)的。綜上，M×M階的復(fù)Hermitian 矩陣Rx的特征值分解等價(jià)于求解由實(shí)矩陣Rr、Ri組成的2M×2M階分塊矩陣的特征值分解。這種轉(zhuǎn)化方法是嚴(yán)格等價(jià)的，沒(méi)有誤差。但是這種方法擴(kuò)大了待求解矩陣的階數(shù)，運(yùn)算量大。

為了減少運(yùn)算量，第二種轉(zhuǎn)化方法將M×M階復(fù)Hermitian 矩陣轉(zhuǎn)化為M×M階的實(shí)對(duì)稱(chēng)矩陣。定義一個(gè)M×M階分塊矩陣U,M為偶數(shù)時(shí)，結(jié)構(gòu)如下：

定義滿(mǎn)足如下等式的復(fù)Hermitian 矩陣P為中心對(duì)稱(chēng)的復(fù)Hermitian 矩陣：

不難證明可以通過(guò)下式將矩陣Rx轉(zhuǎn)化為一個(gè)中心對(duì)稱(chēng)的復(fù)Hermitian 矩陣R[2,4]：

定義一個(gè)由矩陣R變換而來(lái)的Hermitian矩陣K：

綜上，矩陣K是一個(gè)實(shí)對(duì)稱(chēng)矩陣。定義UK為矩陣K的特征向量矩陣，∑K為矩陣K的特征值矩陣，那么矩陣R的特征值分解可以表示為：

根據(jù)式（9）可以得到矩陣R的特征向量矩陣為UHUK，矩陣R的特征值矩陣就是∑K。

如果M為奇數(shù)，此時(shí)，矩陣U結(jié)構(gòu)類(lèi)似，矩陣I、J都是階的，上述證明依舊成立。

綜上，可以通過(guò)一定的運(yùn)算將原本的復(fù)Hermitian 矩陣Rx先變換為一個(gè)中心對(duì)稱(chēng)的復(fù)Hermitian 矩陣R，再將矩陣R等價(jià)變換為一個(gè)實(shí)對(duì)稱(chēng)矩陣K。在求解K的特征值分解后得到R的特征值分解，從而替代Rx的特征值分解。第二種轉(zhuǎn)化方法不擴(kuò)大原矩陣階數(shù)，計(jì)算量相較于第一種方法是大幅度減小的。當(dāng)然這種轉(zhuǎn)化并不是嚴(yán)格等價(jià)的，對(duì)于矩陣R和矩陣Rx之間的轉(zhuǎn)化誤差，可以證明在所有的復(fù)Hermitian中心對(duì)稱(chēng)矩陣中，矩陣R是對(duì)矩陣Rx在歐式距離（Euclidean Distance）下的最佳估計(jì)[12]。

2 Jacobi旋轉(zhuǎn)法

Jacobi 法的基本原理為如果一個(gè)方陣是三角矩陣或?qū)顷?，則該方陣主對(duì)角線上的元素就是該方陣的特征值。因此可以利用旋轉(zhuǎn)矩陣不斷地對(duì)實(shí)對(duì)稱(chēng)矩陣K進(jìn)行正交相似變換，通過(guò)變換將方陣K中所有非主對(duì)角線元素都消去，在精度足夠高的情況下，最后得到的近似的對(duì)角陣就是原矩陣K的特征值矩陣。矩陣K的特征向量矩陣也可以通過(guò)類(lèi)似的運(yùn)算來(lái)求得[13-14]。

首先基于單位矩陣I定義一個(gè)M×M階旋轉(zhuǎn)矩陣Qij(1 ≤i

易證矩陣Qij是一個(gè)正交矩陣，令K0=K，利用旋轉(zhuǎn)矩陣Qij對(duì)矩陣K0進(jìn)行一次正交相似變換，得到矩陣K1：

變換后的矩陣K1仍然是一個(gè)實(shí)對(duì)稱(chēng)矩陣，設(shè)矩陣K1中第p行第q列（1 ≤p

1）當(dāng)p≠i或j且q≠i或j時(shí)，變換前后的矩陣元素相等，不發(fā)生變化。

2）當(dāng)p=i,q=j時(shí)：

3）當(dāng)p=i或j，q≠i或j時(shí)：

為了提高效率，每次正交相似變化都應(yīng)該消去矩陣上三角部分中絕對(duì)值最大的元素。定義實(shí)對(duì)稱(chēng)矩陣G的特征值矩陣CV、特征向量矩陣EV，利用Jacobi 旋轉(zhuǎn)法求解矩陣G的特征值分解，一般步驟如下：

1）給矩陣CV和矩陣EV賦初值（矩陣CV初值為矩陣G，矩陣EV初值為單位矩陣I）；

2）尋找矩陣CV上三角部分中絕對(duì)值最大的元素cvij，將其與允許誤差ej進(jìn)行比較，小于ej進(jìn)行步驟6），否則進(jìn)行步驟3）；

3）根據(jù)元素cvij、cvii、cvjj的值計(jì)算旋轉(zhuǎn)角θ；

4）計(jì)算旋轉(zhuǎn)角θ的正余弦值，構(gòu)造旋轉(zhuǎn)矩陣Qij；

5）計(jì)算變換后的矩陣CV、矩陣EV，然后進(jìn)行步驟2）：

6）輸出矩陣CV和矩陣EV。

Jacobi 旋轉(zhuǎn)法的計(jì)算時(shí)間取決于兩點(diǎn)：一是ej的大小，ej越小表明計(jì)算出來(lái)的特征值矩陣越接近對(duì)角陣，算法誤差越小，需要的迭代次數(shù)就越多。二是待分解矩陣G的大小，G越大，待消去的元素?cái)?shù)量越多，迭代次數(shù)也就越多，并且單次迭代時(shí)矩陣乘法的運(yùn)算量也越大。

3 Jacobi旋轉(zhuǎn)法的FPGA實(shí)現(xiàn)

上一節(jié)的算法流程中，Jacobi 法的迭代次數(shù)不固定。在FPGA 中迭代次數(shù)不固定會(huì)導(dǎo)致時(shí)延的不確定，因此需要根據(jù)待分解矩陣的階數(shù)以及數(shù)據(jù)的量化精度來(lái)固定迭代次數(shù)。對(duì)于ej而言，在FPGA 中會(huì)受到量化精度和計(jì)算精度的影響，過(guò)小的ej是無(wú)法達(dá)到的。根據(jù)仿真，16 位量化精度時(shí)，4 階方陣需要15 次迭代；24 位量化精度時(shí)，6 階方陣需要25 次迭代。

利用FPGA 實(shí)現(xiàn)Jacobi 法求解特征值矩陣時(shí)，需要解決兩個(gè)問(wèn)題：一是三角函數(shù)的計(jì)算，二是矩陣乘法的實(shí)現(xiàn)。

Jacobi 迭代法中三角函數(shù)的計(jì)算主要包括兩個(gè)部分：一是利用反正切arctan 函數(shù)計(jì)算旋轉(zhuǎn)角θ的值，二是計(jì)算旋轉(zhuǎn)角θ的正、余弦值。這兩步計(jì)算都可以采用ISE 軟件中自帶的Cordic IP 核來(lái)完成，計(jì)算精度很高。

Jacobi 法中每次迭代需要進(jìn)行3 次矩陣乘法運(yùn)算。雖然可以用式(12)～(17)來(lái)直接計(jì)算每個(gè)元素值，但每次消去的元素位置不固定，導(dǎo)致每個(gè)元素的表達(dá)式也不固定，不如直接通過(guò)矩陣乘法來(lái)計(jì)算旋轉(zhuǎn)后的矩陣。在FPGA 中，可以采用脈動(dòng)陣列（Systolic Array）的結(jié)構(gòu)來(lái)構(gòu)造矩陣乘法器[16-17]。脈動(dòng)陣列由多個(gè)相同的處理單元（Processing Elements，PE）按照特定的規(guī)則互聯(lián)而成。

首先對(duì)矩陣乘法進(jìn)行拆分，設(shè)m×n階矩陣C是由m×l階矩陣A和l×n階矩陣B相乘得到的：

其中，i=1，2，…，m；j=1，2，…，l。

不難發(fā)現(xiàn)每一個(gè)PE 需要一個(gè)乘法器和一個(gè)累加器來(lái)完成對(duì)應(yīng)元素的計(jì)算。為了使整個(gè)運(yùn)算過(guò)程流水線化，輸入的數(shù)據(jù)需要在PE 間“流動(dòng)”起來(lái)，所以每個(gè)PE 都要對(duì)兩路輸入數(shù)據(jù)做一個(gè)單位的延時(shí)輸出來(lái)提供給下一級(jí)PE 進(jìn)行運(yùn)算[18]。綜上，矩陣乘法器中PE 的框圖如圖1 所示。

圖1 PE結(jié)構(gòu)圖

為了實(shí)現(xiàn)m×l階矩陣A與l×n階矩陣B相乘，需要m×n個(gè)這樣的PE 來(lái)組成脈動(dòng)陣列。每一個(gè)PE 的輸出對(duì)應(yīng)矩陣C中的一個(gè)元素，用脈動(dòng)陣列實(shí)現(xiàn)的矩陣乘法器結(jié)構(gòu)如圖2 所示。

圖2 脈動(dòng)陣列矩陣乘法器結(jié)構(gòu)圖

上圖輸入數(shù)據(jù)中的0 表示對(duì)輸入數(shù)據(jù)做一個(gè)時(shí)鐘的延遲以匹配上一級(jí)PE 的延時(shí)輸出。如果預(yù)先對(duì)輸入數(shù)據(jù)流進(jìn)行控制可以提高PE 的利用率，減少PE 的數(shù)量，從而節(jié)省乘法器和加法器資源。

在矩陣乘法運(yùn)算后，要對(duì)結(jié)果進(jìn)行截位，截位策略是保證未參與運(yùn)算的元素值不發(fā)生變化。根據(jù)仿真結(jié)果，特征值矩陣需要提前進(jìn)行位擴(kuò)展，以保證矩陣對(duì)角線上的特征值不溢出[19]。

4 仿真結(jié)果

以16 位定點(diǎn)數(shù)量化后的4 階復(fù)Hermitian 矩陣Rx為例，其實(shí)部為：

使用第二種矩陣轉(zhuǎn)化方法后用Jacobi 旋轉(zhuǎn)法求特征值分解，其中Jacobi 法的迭代次數(shù)為15 次，Matlab 軟件中計(jì)算出的特征值矩陣CV如下：

在ISE 軟件中編寫(xiě)好程序后，用Modelsim 仿真得到的特征值矩陣CVM和特征向量矩陣EVM，如圖3所示。模塊端口圖如圖4 所示。

圖3 Modelsim仿真結(jié)果圖

圖4 模塊端口圖

其中，clk 為輸入時(shí)鐘端口，en 為模塊輸入使能端口，rst_n 為模塊異步清零端口，Snr_in1～Snr_in4 是矩陣Rx實(shí)部輸入端口，Sni_in1～Sni_in4 是矩陣Rx虛部輸入端口，每一個(gè)端口輸入待求解矩陣的一行數(shù)據(jù)；valid 端口為輸出有效端口，計(jì)算結(jié)果輸出時(shí)同步置高，CV_out1～CV_out4 是矩陣CVM的輸出端口（特征值矩陣在迭代開(kāi)始時(shí)就已經(jīng)做了防溢出位擴(kuò)展），EVr_out1～EVr_out4 是矩陣EVM實(shí)部輸出端口，EVi_out1～EVi_out4 是矩陣EVM虛部輸出端口。整個(gè)模塊從數(shù)據(jù)輸入到結(jié)果輸出的延時(shí)為1 074 個(gè)時(shí)鐘周期。由于Cordic IP 核的計(jì)算存在細(xì)小的誤差，導(dǎo)致仿真結(jié)果與Matlab 定點(diǎn)仿真結(jié)果有一定的誤差，不過(guò)對(duì)整體精度影響非常小?！?‖?2表示對(duì)矩陣取二范數(shù)，Modelsim 仿真結(jié)果與Matlab 仿真結(jié)果的歸一化誤差為：

整個(gè)模塊的資源消耗情況如表1 所示。

表1 FPGA模塊資源消耗

5 結(jié)束語(yǔ)

該文對(duì)矩陣特征值分解算法進(jìn)行了研究并給出了求解復(fù)Hermitian 矩陣的特征值分解的一般方案。該方案利用FPGA 完成復(fù)Hermitian 矩陣到實(shí)對(duì)稱(chēng)矩陣的轉(zhuǎn)化后，再運(yùn)用Jacobi 旋轉(zhuǎn)法求解轉(zhuǎn)化后的實(shí)對(duì)稱(chēng)矩陣特征值分解。整個(gè)方案具有結(jié)構(gòu)簡(jiǎn)單、計(jì)算精度高、時(shí)延小、資源消耗少的特點(diǎn)。