ESPRIT 算法廣義逆矩陣求解的快速FPGA 實(shí)現(xiàn)

王衛(wèi)江，張拓鋒，蔣榮堃,2,3，李澤英，王曉華，譚志昕，薛丞博,2

（1. 北京理工大學(xué) 集成電路與電子學(xué)院，北京 100081；2. 北京理工大學(xué) 重慶創(chuàng)新中心，重慶 401120；3. 北京理工大學(xué) 重慶微電子中心，重慶 401332；4. 北京遙測(cè)技術(shù)研究所，北京 100094）

相較于多信號(hào)分類（multiple signal classification，MUSIC）算法[1]等其他空間譜估計(jì)算法，基于旋轉(zhuǎn)不變子空間的信號(hào)參數(shù)估計(jì)（estimating signal parameters via rotational invariance techniques ，ESPRIT）算法由于省去了譜峰搜索這一步驟，減少了很大一部分計(jì)算量且仍舊能夠保持良好的精確性和穩(wěn)定性. 但由于ESPRIT 算法需要對(duì)復(fù)數(shù)矩陣進(jìn)行特征值分解和廣義特征值分解，依舊具有很大的計(jì)算量. 而此類高分辨測(cè)向系統(tǒng)一般都對(duì)實(shí)時(shí)性有較高要求，僅僅依靠軟件實(shí)現(xiàn)可能難以滿足現(xiàn)實(shí)應(yīng)用的需求，因此使用硬件實(shí)現(xiàn)來(lái)完成系統(tǒng)的加速是時(shí)下空間譜估計(jì)領(lǐng)域一個(gè)熱點(diǎn)[2-4].

ESPRIT 算法的基本原理是利用子空間的旋轉(zhuǎn)不變性來(lái)求解信號(hào)的入射角度[5]，其中比較復(fù)雜的部分是利用最小二乘法來(lái)求解信號(hào)子空間的廣義特征值，這一部分也可以轉(zhuǎn)化成求解信號(hào)子空間矩陣的廣義逆矩陣. 由于信號(hào)子空間矩陣為復(fù)數(shù)矩陣且階數(shù)大，其廣義逆矩陣的求解運(yùn)算量較大，使用MATLAB 等軟件很難達(dá)到實(shí)時(shí)性要求. 而FPGA 平臺(tái)以其功耗低、計(jì)算速度快等優(yōu)點(diǎn)使得其目前越來(lái)越多的被應(yīng)用于矩陣的各種運(yùn)算[6]，如定點(diǎn)、浮點(diǎn)運(yùn)算、特征值求解、矩陣求逆、矩陣相乘等.

目前基于FPGA 的矩陣求逆設(shè)計(jì)大部分都只應(yīng)用于非奇異矩陣，而針對(duì)一般矩陣的廣義逆求解的硬件實(shí)現(xiàn)很少有人研究. 在現(xiàn)有的硬件廣義逆求解研究中，SAHOO 等[7]提出了一種基于正交三角（QR）分解的廣義逆矩陣求解算法，但該算法注重提高精確度，在速度上沒有過(guò)高要求. 應(yīng)俊等[8]提出使用基于坐標(biāo)旋轉(zhuǎn)數(shù)字計(jì)算機(jī)（coordinate rotation digital computer，CORDIC）算法來(lái)對(duì)奇異值（SVD）分解進(jìn)行優(yōu)化，從而降低硬件資源的消耗，但同樣引入了其他計(jì)算單元，較為復(fù)雜. 陳曉東等[9]使用Cholesky 分解來(lái)計(jì)算矩陣求逆，該算法計(jì)算相較于SVD 分解和QR分解計(jì)算簡(jiǎn)單，但同樣涉及大量開方運(yùn)算，不利于硬件實(shí)現(xiàn). 周金強(qiáng)等[10]對(duì)Cholesky 分解進(jìn)行了改進(jìn)，通過(guò)引入輔助矩陣的方法規(guī)避掉了開方運(yùn)算，具有更好的實(shí)時(shí)性能，但只能處理共軛對(duì)稱矩陣，局限性較大.

本文針對(duì)ESPRIT 算法中信號(hào)子空間矩陣的特點(diǎn)，提出了使用廣義逆公式來(lái)構(gòu)建一個(gè)一般復(fù)數(shù)矩陣的廣義逆求解系統(tǒng). 在FPGA 上進(jìn)行了硬件實(shí)現(xiàn)，最后結(jié)合MATLAB 仿真結(jié)果對(duì)該系統(tǒng)的性能進(jìn)行了分析. 經(jīng)驗(yàn)證，該系統(tǒng)計(jì)算簡(jiǎn)單，相較于其他方法減少了硬件消耗，并且在保證結(jié)果精確度的同時(shí)有效減少了運(yùn)算時(shí)間.

1 矩陣的廣義逆求解

廣義逆矩陣是逆矩陣的推廣形式. 逆矩陣只針對(duì)非奇異矩陣，而奇異矩陣或者非方陣并不存在逆矩陣，因此引入廣義逆矩陣的概念,在MATLAB 中可以使用函數(shù)B= pinv（A）來(lái)進(jìn)行求解，結(jié)果返回一個(gè)與矩陣A轉(zhuǎn)置同型的矩陣B且滿足ABA=A，BAB=B，稱矩陣B為矩陣A的廣義逆.

1.1 廣義逆矩陣的定義

1.2 廣義逆矩陣的解法

對(duì)于矩陣的廣義逆求解，一般有公式法、奇異值（SVD）分解法和正交三角（QR）分解法.

1.2.1 公式法

稱為右逆；

如果A不滿秩，則只能進(jìn)行奇異值分解，使用SVD分解法.

1.2.2 SVD 分解法

因?yàn)橛暇仃嚨哪婢仃嚨扔谄涔曹椶D(zhuǎn)置矩陣，對(duì)角矩陣的逆為對(duì)角線元素取倒數(shù)，則矩陣A的廣義逆為

1.2.4 3 種方法的比較與選擇

3 種方法中，公式法可以看成是對(duì)求解方程組Ax=b的推廣，當(dāng)矩陣A滿秩（行滿秩或列滿秩）時(shí)都可以使用公式法計(jì)算廣義逆矩陣；QR 分解也只應(yīng)用在矩陣滿秩的情況下，一般常用的QR 分解算法如Householder 變換、Givens 旋轉(zhuǎn)和Schmidt 正交化等計(jì)算都較為復(fù)雜；當(dāng)矩陣A不滿秩時(shí)只能選用SVD 分解法通過(guò)奇異值求解廣義逆矩陣，這一方法對(duì)于任何A∈Cm×n都適用. 但在計(jì)算矩陣奇異值分解時(shí)，隨著矩陣的增大計(jì)算量會(huì)以指數(shù)型增長(zhǎng)，使得系統(tǒng)的計(jì)算時(shí)間大大增加. MARKKANDAN 等[12]中對(duì)SVD分解和QR 分解算法的計(jì)算復(fù)雜度進(jìn)行了詳細(xì)的分析，QR 分解相較于SVD 分解計(jì)算更加簡(jiǎn)單.

在ESPRIT 算法中需要使用最小二乘法，最小二乘的本質(zhì)就是一個(gè)超定方程組的求解問題，超定方程組的求解方式之一就是通過(guò)求解廣義逆的形式，即求解信號(hào)子空間矩陣的廣義逆. 由于信號(hào)子空間矩陣由若干個(gè)線性無(wú)關(guān)的信號(hào)子空間向量組成，因此必定為列滿秩矩陣，因此可以使用式（4）直接進(jìn)行求解. 該方法相較于SVD 分解和QR 分解復(fù)雜計(jì)算更少，更加適合用于硬件實(shí)現(xiàn)，因此本設(shè)計(jì)選擇使用公式法求解廣義逆矩陣.

2 公式法求解過(guò)程

公式法求解信號(hào)子空間矩陣廣義逆步驟分為：

（1）輸入信號(hào)子空間矩陣S∈Cm×n；

（3）計(jì)算 (SHS)-1， 即求解SHS的逆矩陣；

（4）矩陣 (SHS)-1與 矩陣SH相乘得到信號(hào)子空間矩陣的廣義逆矩陣S+.

矩陣SH為矩陣S的共軛轉(zhuǎn)置，因此相乘后的矩陣SHS為厄米特矩陣，其主對(duì)角線上的元素都是實(shí)數(shù)，矩陣中每一個(gè)第i行第j列的元素都與第j行第i列的元素的共軛相等. 針對(duì)方陣求逆的問題一般常用的方法為Cholesky 分解法、QR 分解法和LU 分解法.其中Cholesky 分解法雖然計(jì)算量較小，但只能對(duì)正定矩陣進(jìn)行分解，適用范圍過(guò)??； QR 分解雖然適用于任意型的矩陣，應(yīng)用范圍廣泛，但 QR 分解的過(guò)程過(guò)于復(fù)雜，不易于實(shí)現(xiàn)硬件化；LU 分解雖然也對(duì)矩陣有一定要求，但處理的矩陣SHS必定為非奇異矩陣，滿足LU 分解的充分條件，且計(jì)算格式簡(jiǎn)單，計(jì)算量小，易于硬件實(shí)現(xiàn). MARKKANDAN 等[12]中同樣對(duì)LU 分解的計(jì)算復(fù)雜度進(jìn)行了分析，隨著矩陣維度的增大，其所消耗的資源相較于QR 分解算法越來(lái)越小，因此本文使用LU 分解來(lái)求解矩陣SHS的逆矩陣.

2.1 LU 分解

設(shè)矩陣A∈Cn×n， 將矩陣A分解成一個(gè)單位下三角矩陣L和一個(gè)上三角矩陣U的乘積，即A=LU，該種分解被稱為Doolittle 分解或LU 分解[13].

2.2 下三角矩陣求逆

使用LU 分解的矩陣求逆運(yùn)算表達(dá)式為

求解出L和U之后，需要對(duì)單位下三角矩陣L和上三角矩陣U分別求逆，并求U-1和L-1的乘積. 采用初等變換法可以得到三角矩陣的求逆公式. 對(duì)于n階下三角矩陣L，與n階單位矩陣E組成增廣矩陣（L|E）：

由于L、U矩陣都為可逆矩陣，根據(jù)矩陣轉(zhuǎn)置的逆矩陣等于矩陣逆的轉(zhuǎn)置矩陣這一性質(zhì)：

可以先將上三角陣U轉(zhuǎn)置成與L同型的下三角矩陣再求其逆矩陣，再將結(jié)果取轉(zhuǎn)置就能得到U的逆矩陣，這樣硬件實(shí)現(xiàn)時(shí)更加快速簡(jiǎn)潔.

2.3 矩陣乘法

3 FPGA 系統(tǒng)設(shè)計(jì)

基于上述公式法求解信號(hào)子空間矩陣廣義逆的算法，將FPGA 硬件設(shè)計(jì)分為了3 個(gè)主要模塊：矩陣乘法模塊Ⅰ、LU 分解求逆模塊和矩陣乘法模塊Ⅱ.FPGA 系統(tǒng)的結(jié)構(gòu)框圖如圖1 所示，第一個(gè)矩陣乘法模塊接收初始信號(hào)子空間矩陣S的信息，計(jì)算共軛轉(zhuǎn)置矩陣SH與原矩陣S的乘法并將結(jié)果發(fā)送到LU分解求逆模塊；LU 分解求逆子模塊將接收到的矩陣分解成L和U矩陣，對(duì)L和U矩陣分別求逆再相乘得到輸入矩陣的逆矩陣，并將逆矩陣發(fā)送給第二個(gè)矩陣乘法模塊；矩陣乘法模塊Ⅱ計(jì)算逆矩陣(SHS)-1與矩陣SH相 乘，最終得到的結(jié)果就是廣義逆矩陣S+.各模塊中矩陣的存儲(chǔ)使用的都是FPGA 中BRAM IP核的雙口RAM，根據(jù)矩陣元素的行列信息來(lái)獲得對(duì)應(yīng)數(shù)據(jù)在RAM 中的地址，以此來(lái)完成數(shù)據(jù)的讀寫和運(yùn)算.

廣義逆求解系統(tǒng)設(shè)計(jì)的模塊示意圖如圖1 所示.

圖1 廣義逆求解系統(tǒng)結(jié)構(gòu)框圖Fig. 1 Structure block diagram of generalized inverse solution system

在FPGA 中數(shù)據(jù)有定點(diǎn)數(shù)和浮點(diǎn)數(shù)兩種表示方式. 浮點(diǎn)數(shù)能夠表達(dá)的數(shù)值范圍更廣，但是與此同時(shí)計(jì)算的復(fù)雜度也更高并且誤差會(huì)隨著數(shù)據(jù)的增大而增大；定點(diǎn)數(shù)能夠表示的數(shù)值范圍相對(duì)較小但誤差恒定不會(huì)隨著數(shù)據(jù)的變化而變化. 尤其涉及到硬件實(shí)現(xiàn)時(shí)，浮點(diǎn)運(yùn)算十分復(fù)雜、會(huì)占用大量的資源. 因此在本設(shè)計(jì)中的數(shù)據(jù)都使用定點(diǎn)數(shù)表示，數(shù)據(jù)使用二進(jìn)制表示，包括符號(hào)位、整數(shù)位和小數(shù)位3 部分，涉及到的復(fù)數(shù)的加減乘除運(yùn)算也都為定點(diǎn)運(yùn)算.

3.1 矩陣乘法模塊

矩陣乘法模塊的主要運(yùn)算單元由一系列乘累加組成，計(jì)算量較大、耗時(shí)也比較長(zhǎng). 在本次FPGA 硬件設(shè)計(jì)中，使用狀態(tài)機(jī)來(lái)進(jìn)行乘累加計(jì)算的控制，當(dāng)狀態(tài)機(jī)跳轉(zhuǎn)到工作狀態(tài)時(shí)，求解乘積矩陣AB的第i行j列元素，需要將矩陣A的第i行所有元素與矩陣B第j列的所有元素一一對(duì)應(yīng)相乘并求和. 計(jì)算完畢后狀態(tài)機(jī)將轉(zhuǎn)入空閑狀態(tài)直到下一個(gè)使能信號(hào)到來(lái)，開始計(jì)算乘積矩陣AB的下一個(gè)元素. 圖2 為矩陣乘模塊硬件結(jié)構(gòu).

圖2 矩陣乘模塊硬件結(jié)構(gòu)Fig. 2 Hardware structure of matrix multiplication module

矩陣乘法模塊Ⅰ的功能為計(jì)算信號(hào)子空間矩陣S的共軛轉(zhuǎn)置矩陣SH與其自身的乘積. 由于矩陣SH是矩陣S的共軛轉(zhuǎn)置，相乘后得到的矩陣SHS為厄米特矩陣，因此在計(jì)算時(shí)可以直接省去對(duì)應(yīng)主對(duì)角線元素中虛部的計(jì)算過(guò)程，以此減少部分計(jì)算量，節(jié)約了相對(duì)應(yīng)的資源. 該模塊首先讀取存儲(chǔ)在RAM_S 和RAM_S_i 的信號(hào)子空間矩陣，其中RAM_S 存儲(chǔ)信號(hào)子空間矩陣的實(shí)部，RAM_S_i 存儲(chǔ)信號(hào)子空間矩陣的虛部. 經(jīng)過(guò)乘累加計(jì)算后將得到的矩陣SHS的元素存儲(chǔ)到RAM_SHS中.

3.2 LU 分解求逆模塊

根據(jù)上文得到的式（12）來(lái)設(shè)計(jì)LU 分解求逆硬件實(shí)現(xiàn)的算法. 首先需要生成L、U矩陣的行列信息，從第一行第一列開始（i=1，j=1）計(jì)算矩陣L和U的第一個(gè)元素l11，u11. 計(jì)算完畢后列數(shù)加一繼續(xù)計(jì)算第一行第二列的元素l12，u12，當(dāng)列數(shù)增加到大于矩陣維數(shù)時(shí)行數(shù)加一，列數(shù)置1. 以此類推，直到行數(shù)也大于矩陣維數(shù)時(shí)便可以將矩陣三角分解后的L和U矩陣完全求解. 由于L和U分別是下三角矩陣和上三角矩陣，這兩個(gè)矩陣對(duì)應(yīng)位置上的元素一般總有一個(gè)是0，所以在計(jì)算元素前首先對(duì)當(dāng)前的行列信息進(jìn)行判斷以此來(lái)減少不必要的計(jì)算：

圖3 LU 分解求逆模塊結(jié)構(gòu)框圖Fig. 3 Structure block diagram of LU decomposition inversion module

（1）如果行列相等（i=j），此時(shí)矩陣L元素li j=1，矩陣U元素uij不為0；

（2）如果行數(shù)大于列數(shù)（i＞j），此時(shí)矩陣L元素li j不為0，而矩陣U元素ui j=0；

（3）如果行數(shù)小于列數(shù)（i＜j），此時(shí)矩陣L元素li j=0， 而矩陣U元素uij不為0.

此外根據(jù)遞推公式可以看出在計(jì)算矩陣L和U當(dāng)前行列的元素時(shí)需要用到對(duì)應(yīng)位置的原矩陣的元素以及此行列前求解的L、U矩陣的元素. 因此可以根據(jù)行列信息來(lái)產(chǎn)生地址值，如使用行數(shù)i乘以當(dāng)前計(jì)算的元素位數(shù)s再加上列數(shù)j，以此來(lái)讀取存儲(chǔ)在BRAM 中的L、U矩陣和原矩陣SHS中所對(duì)應(yīng)的元素.

4 仿真結(jié)果與分析

4.1 廣義逆求解模塊功能仿真

本次設(shè)計(jì)中軟件程序的編寫和仿真都使用MATLAB R2020a 在Intel(R) Core(TM) i5-3470 處 理 器 上完成，內(nèi)存為8 GB、主頻為3.20 GHz. 而系統(tǒng)的硬件編譯和仿真都在Vivado 平臺(tái)上進(jìn)行，F(xiàn)PGA 選用的是Virtex-7 的xc7vx485tffg1157-1 開發(fā)板，工作頻率為200 MHz. 使用MATLAB 隨機(jī)生成用于測(cè)試的信號(hào)子空間矩陣，該矩陣為行列數(shù)分別為59 和15 的復(fù)矩陣.

4.2 FPGA 邏輯資源分析

廣義逆求解模塊的LUT（Look-Up-Table）、Block RAM、Registers 和DSP（digital signal processor）的資源占用情況如表1 所示. 根據(jù)表格可以看出LU 分解求逆模塊需要的邏輯資源是最大的，而矩陣乘法模塊需要的資源較少，僅為前者的1/4. 這是由于LU 分解求逆涉及到了大量的除法運(yùn)算，相較于僅是大量乘累加運(yùn)算的矩陣乘運(yùn)算占用的資源更多.

表1 各模塊所占用的邏輯資源Tab. 1 Logical resources of each module

4.3 運(yùn)算時(shí)間分析

通過(guò)硬件仿真可以得出使用FPGA 的設(shè)計(jì)在求解矩陣廣義逆時(shí)用時(shí)約為2.18 ms. 各模塊具體的工作時(shí)間如表2 所示, 可以看出矩陣乘模塊運(yùn)行的時(shí)間較長(zhǎng)，約占總體時(shí)間的80.7%.

表2 各模塊運(yùn)行時(shí)間對(duì)比Tab. 2 Operation time of each module

這是由于矩陣乘法模塊需要計(jì)算大量的乘累加運(yùn)算，需要較長(zhǎng)的運(yùn)算時(shí)間且矩陣越大需要的時(shí)間也就越長(zhǎng). 而在MATLAB 中使用pinv 函數(shù)對(duì)同樣的矩陣進(jìn)行廣義逆求解，計(jì)算5 次取平均時(shí)間約為15.7 ms.可以得出當(dāng)FPGA 工作在200 MHz 的時(shí)鐘下時(shí)，求解廣義逆矩陣的時(shí)間比軟件端快了大約7.2 倍，可以節(jié)省大約86%的計(jì)算時(shí)間.

4.4 ESPRIT 算法結(jié)果分析

將FPGA 求出的廣義逆矩陣重新輸入到MATLAB 中完成ESPRIT 算法的后續(xù)步驟并且與在MATLAB 中使用pinv 函數(shù)求得信號(hào)子空間矩陣的廣義逆的ESPRIT 算法最終的結(jié)果進(jìn)行誤差分析如表3，可以看出最終所得角度的平均誤差大概為0.04°，滿足工程的需求.

表3 ESPRIT 算法結(jié)果及誤差分析Tab. 3 ESPRIT algorithm results and error analysis

4.5 與其他算法對(duì)比

應(yīng)俊等[8]中使用了CORDIC 算法對(duì)矩陣進(jìn)行SVD 分解，系統(tǒng)運(yùn)行的總時(shí)間為547.6 μs. 使用行列數(shù)都為8 的輸入矩陣在與此文獻(xiàn)相同的工作頻率250 MHz 下對(duì)本設(shè)計(jì)進(jìn)行仿真，最終得到本系統(tǒng)運(yùn)行的時(shí)間為139.4 μs. 可以看出本設(shè)計(jì)使用的算法相較于傳統(tǒng)的SVD 分解法在計(jì)算廣義逆矩陣上具有更好的時(shí)間性能.

周金強(qiáng)等[10]中使用了改進(jìn)Cholesky 分解來(lái)求解廣義逆矩陣. 使用行列數(shù)都為4 的信號(hào)子空間矩陣，對(duì)本設(shè)計(jì)進(jìn)行仿真，在矩陣維度和工作頻率都相同的條件下，兩種算法的邏輯資源和系統(tǒng)運(yùn)算時(shí)間如表4.

表4 邏輯資源與計(jì)算時(shí)間Tab. 4 Logical resources and computing time

由表4 可以看出本設(shè)計(jì)相較于改進(jìn)的Cholesky算法消耗的邏輯資源更少但是計(jì)算時(shí)間更長(zhǎng). 但文獻(xiàn)[10]中的算法僅能計(jì)算共軛對(duì)稱矩陣，局限性較大，而本設(shè)計(jì)的算法能夠計(jì)算任意矩陣的廣義逆矩陣，具有更多的應(yīng)用場(chǎng)景.

5 結(jié) 論

本文針對(duì)ESPRIT 算法中信號(hào)子空間矩陣的特點(diǎn)，利用廣義逆矩陣的性質(zhì)提出了一種新的矩陣廣義逆求解算法. 這種算法比一般的QR 分解和SVD分解法計(jì)算更加簡(jiǎn)單，用時(shí)更短. 利用FPGA 技術(shù)對(duì)該廣義逆矩陣求解系統(tǒng)進(jìn)行了硬件實(shí)現(xiàn)，將整個(gè)系統(tǒng)分為了矩陣乘法模塊Ⅰ、LU分解求逆模塊和矩陣乘法模塊Ⅱ共3 個(gè)模塊，對(duì)每個(gè)模塊分別進(jìn)行了設(shè)計(jì)和測(cè)試，給出了總體設(shè)計(jì)框圖. 并且使用Vivado對(duì)廣義逆矩陣求解系統(tǒng)進(jìn)行了硬件的仿真得到時(shí)序圖和計(jì)算結(jié)果. 將其與MATLAB 仿真的計(jì)算結(jié)果和平均時(shí)間進(jìn)行對(duì)比，發(fā)現(xiàn)使用FPGA 的設(shè)計(jì)在求解矩陣廣義逆時(shí)可以節(jié)省約86%的計(jì)算時(shí)間，且最終所得角度的平均誤差大概為0.04°. 此結(jié)果可以驗(yàn)證該硬件系統(tǒng)在滿足一般工程精度需求的同時(shí)可以一定程度上減少計(jì)算所需時(shí)間，具有更好的實(shí)時(shí)性.