多子陣波束域Root-MUSIC算法并行處理實(shí)現(xiàn)方法

劉千里

（海軍駐武漢四三八廠軍事代表室，武漢，430060）

多子陣波束域Root-MUSIC算法并行處理實(shí)現(xiàn)方法

劉千里

（海軍駐武漢四三八廠軍事代表室，武漢，430060）

在對MSB-RMU算法計(jì)算量分析的基礎(chǔ)上，研究了基于并行Jacobi方法的特征值分解算法，并設(shè)計(jì)了一種基于FPGA的采用雙邊CORDIC結(jié)構(gòu)的并行Jacobi實(shí)現(xiàn)方法，為特征分解類算法提供了一種高效的單芯片解決方案。最后討論了MSB-RMU算法在一塊由6片DSP和2片F(xiàn)PGA組成的并行信號處理平臺的實(shí)現(xiàn)過程。試驗(yàn)結(jié)果表明，設(shè)計(jì)的并行處理方法能夠滿足多子陣波束域Root-MUSIC算法的實(shí)時(shí)實(shí)現(xiàn)。

多子陣波束域；Root-MUSIC；并行Jacobi；FPGA；DSP

對于海底地形地貌聲吶探測系統(tǒng)來說，關(guān)鍵的要素之一就是算法的研究與實(shí)現(xiàn)。因此，除了算法本身的高性能之外，考察海底探測算法另一個(gè)重要的標(biāo)準(zhǔn)就是其是否具有可實(shí)現(xiàn)性，即算法是否能實(shí)時(shí)實(shí)現(xiàn)。

多子陣波束域CAATI算法（MSB-CAATI）和多子陣波束域Root-MUSIC算法（MSB-RMU），兩者都能夠同時(shí)估計(jì)回波的方向和強(qiáng)度，非常適用于海底地形地貌探測系統(tǒng)［1］。但后者具有更好的估計(jì)性能，且實(shí)現(xiàn)過程較前者簡單［2］。本文在對MSB-RMU算法計(jì)算量分析的基礎(chǔ)上，研究了基于并行Jacobi方法的特征值分解算法，并設(shè)計(jì)了一種基于FPGA的采用雙邊CORDIC結(jié)構(gòu)的并行Jacobi實(shí)現(xiàn)方法，為特征分解類算法提供了一種高效的單芯片解決方案。最后討論了MSB-RMU算法在一塊由6片DSP和2片F(xiàn)PGA組成的并行信號處理平臺的實(shí)現(xiàn)過程。

1 算法計(jì)算量分析

MSB-RMU算法中計(jì)算量主要集中在多個(gè)子陣的波束形成、子陣波束域輸出的相關(guān)矩陣運(yùn)算、相關(guān)矩陣的特征分解以及多項(xiàng)式求根得到DOA，再根據(jù)陣元域協(xié)方差矩陣估計(jì)回波強(qiáng)度就可以進(jìn)行海底成像。MSB-RMU算法較傳統(tǒng)Root-MUSIC算法增加了多個(gè)子陣的波束形成，但相關(guān)矩陣階數(shù)減少，從而特征分解的運(yùn)算量大大減小。由此可見，對于大陣元數(shù)系統(tǒng)而言，總的運(yùn)算量要低于傳統(tǒng)算法。即便如此，算法還是具有相當(dāng)?shù)倪\(yùn)算量，采用單獨(dú)的處理器仍難以實(shí)時(shí)實(shí)現(xiàn)，因此還必須采用大規(guī)模并行處理平臺進(jìn)行計(jì)算。

MSB-RMU算法將MUSIC算法的高分辨性能和FT波束形成算法的魯棒性有機(jī)結(jié)合，同時(shí)算法將陣元域處理轉(zhuǎn)換到多個(gè)子陣的波束域處理，從而大大減小了對平均協(xié)方差矩陣RM的特征值分解運(yùn)算量。盡管如此，對于大陣列情況，當(dāng)子陣個(gè)數(shù)較多時(shí)，平均協(xié)方差矩陣RM的階數(shù)仍然較大，也就成為了整個(gè)算法中運(yùn)算量相對集中的部分。因此高效地實(shí)現(xiàn)對平均協(xié)方差矩陣的特征值分解是MSB-RMU算法實(shí)時(shí)實(shí)現(xiàn)的關(guān)鍵。

2 并行Jacobi方法實(shí)現(xiàn)特征值分解（EVD）

2.1 Jacobi方法特征值分解

Jacobi方法簡單緊湊、精度高、收斂較快，是計(jì)算實(shí)對稱矩陣全部特征值和相應(yīng)特征向量的有效方法［3］。Jacobi算法采用平面旋轉(zhuǎn)矩陣進(jìn)行正交相似變換，將實(shí)對稱矩陣化為對角矩陣，用于求解實(shí)對稱矩陣的全部特征值和特征向量。下面給出了基于Jacobi方法的特征分解基本形式。

其中，R為N×N實(shí)對稱矩陣，Λ為對角矩陣，其對角線上的元素即為矩陣R的特征值，I為單位矩陣，U矩陣的列向量為所有特征值對應(yīng)的特征向量。Wpq為平面旋轉(zhuǎn)矩陣：

2.2 并行Jacobi方法特征值分解

并行Jacobi算法［4］的基本思路是每次旋轉(zhuǎn)消去多于一對的非對角元素。若矩陣R為N×N實(shí)對稱矩陣，令m=N/2，顯然R有N（N-1）/2對非對線元素需要消去，每次正交變換可同時(shí)消去的最大數(shù)為N/2。把每2m-1變換記為一次掃描，并行Jacobi算法若能滿足下面兩個(gè)條件，并行計(jì)算的效率將達(dá)到最高。

（1）構(gòu)造正交變換矩陣Qj，使每一個(gè)變換都能消去N/2對非對角元素。

（2）每次掃描，任一對非對角元素只被消去一次，即在一次掃描中，2m-1個(gè)正交變換的每一個(gè)變換消去的元素，都不同于其他變換消去的元素。

對于Qj的選取有很多種方法，文獻(xiàn)［5］給出一種選取方法。對每一次掃描，其2m-1個(gè)正交變換中的每一個(gè)變換Qj的非零元素按以下方法選取。

其中（p，q）的按照下面定義：

并行Jacobi算法每次正交相似變換可同時(shí)消去N/2對非對角元素，因此與經(jīng)典Jacobi算法相比，在Jacobi掃描次數(shù)相同的情況下，所用時(shí)間約為經(jīng)典Jacobi算法的2/N。

3 基于并行Jacobi方法的EVD在FPGA上的實(shí)現(xiàn)

利用CORDIC算法來實(shí)現(xiàn)Jacobi算法已被廣泛的研究［6］，而利用雙邊CORDIC在FPGA上實(shí)現(xiàn)并行Jacobi算法可以得到更高的效率?；?.2的分析，可以把每次正交相似變換記為一次消去過程，每一個(gè)變換消去的元素，都不同于其他變換消去的元素，因此消去過程可同時(shí)進(jìn)行。每個(gè)消去過程可由一個(gè)雙邊CORDIC來實(shí)現(xiàn)，因此需要N/2個(gè)雙邊CORDIC組成一個(gè)并行處理網(wǎng)絡(luò)。雙邊CORDIC解決如式（12）所示的雙邊旋轉(zhuǎn)問題。

單個(gè)雙邊CORDIC處理器原理框圖如圖1所示。

圖1 雙邊CORDIC原理框圖

圖中nσ為加法器的符號控制位，根據(jù)具體情況來選擇，針對Jacobi正交相似變換，符號的選取原則是使待消去的元素逐漸變?yōu)榱恪?/p>

特征向量的計(jì)算也采用CORDIC來實(shí)現(xiàn)，所不同的是每次旋轉(zhuǎn)過程中的符號由消去過程所確定。由于并行Jacobi算法是用于計(jì)算實(shí)對稱矩陣特征分解的方法，因此對于Hermite矩陣而言，還需經(jīng)過下面的變換［7］。

將N階Hermite矩陣H寫成實(shí)數(shù)矩陣和復(fù)數(shù)矩陣兩部分，即

可見，H矩陣的特征分解問題轉(zhuǎn)化為實(shí)對稱矩陣的特征分解問題。

由于CORDIC算法僅需移位和加法運(yùn)算即可完成，使得硬件實(shí)現(xiàn)起來非常簡單。XILINX公司的IP核生成軟件Core Generator中自帶CORDIC算法的IP核，該處理核的結(jié)構(gòu)如圖2所示。

圖2 CORDIC核結(jié)構(gòu)

該CORDIC核具有兩種可選的結(jié)構(gòu)配置選項(xiàng)：一種是完全并行配置模式，該模式具有單周期數(shù)據(jù)輸入輸出，但布局布線所占用的面積較大；另一種是字串行實(shí)現(xiàn)模式，該模式的數(shù)據(jù)輸入輸出需要多個(gè)時(shí)鐘周期，但布局布線所占用的面積較小。算法需要產(chǎn)生一個(gè)尺度因子來對輸出結(jié)果的幅度進(jìn)行校正，CORDIC核提供自動(dòng)補(bǔ)償?shù)倪x項(xiàng)來進(jìn)行幅度校正。CORDIC算法過程中，每個(gè)單位精度需要近似一次移位-加減操作。若選擇字串行結(jié)構(gòu)配置，N的輸出寬度需要N個(gè)時(shí)鐘周期的延遲，并在每N個(gè)時(shí)鐘周期后產(chǎn)生一個(gè)新的輸出；若選擇并行結(jié)構(gòu)配置，輸出寬度需要N個(gè)時(shí)鐘周期的延遲，并在每一個(gè)時(shí)鐘周期后產(chǎn)生一個(gè)新的輸出。因此在設(shè)計(jì)時(shí)應(yīng)注意時(shí)鐘延遲的影響，并根據(jù)具體工程需要在速度和實(shí)現(xiàn)面積上進(jìn)行選擇。

4 多子陣波束域Root-MUSIC算法的并行處理系統(tǒng)設(shè)計(jì)

4.1 并行處理系統(tǒng)組成

為了滿足系統(tǒng)運(yùn)算量需求，設(shè)計(jì)了一塊由6片DSP和2片F(xiàn)PGA組成的CPCI結(jié)構(gòu)并行處理系統(tǒng)，系統(tǒng)框圖如圖3所示。

圖3 并行處理平臺結(jié)構(gòu)框圖

系統(tǒng)主要完成對多個(gè)通道原始數(shù)據(jù)的正交變換和基于MSB-RMU算法的DOA估計(jì)，其中正交變換、多個(gè)子陣波束形成以及特征分解都由FPGA來完成，自相關(guān)矩陣的計(jì)算和最后的多項(xiàng)式求根由DSP陣列完成。下面按照信號流的順序分別討論各主要信號處理部分的實(shí)現(xiàn)方法。

4.2 基于CORDIC的數(shù)字下變頻器

數(shù)字下變頻器（DDC）的作用是將模數(shù)變換器采集到的原始數(shù)據(jù)載波頻率搬移到基帶進(jìn)行處理。傳統(tǒng)的DDC由正余弦查找表、乘法器和低通濾波器組成，如圖4（a）所示。也可以采用CORDIC算法加以實(shí)現(xiàn)，從而無需乘法器和存儲正余弦表。CORDIC算法由Volder［8］最早提出用于計(jì)算平面直角坐標(biāo)系和極坐標(biāo)系之間的變換。其基本思想是用一系列固定角度的不斷偏擺來逼近所需旋轉(zhuǎn)的角度，適當(dāng)選取一些固定的角度值，可以使CORDIC運(yùn)算只需進(jìn)行移位和加減操作，適合硬件實(shí)現(xiàn)。所以在本系統(tǒng)中采用FPGA來實(shí)現(xiàn)DDC。利用CORDIC算法實(shí)現(xiàn)的DDC原理框圖如圖4（b）所示。

圖4 DDC實(shí)現(xiàn)框圖

4.3 基于DSP陣列的協(xié)方差矩陣計(jì)算

完成多子陣波束形成后，將所有子陣大于平均波束能量的輸出波束送給4個(gè)DSP組成的DSP陣列，每個(gè)DSP從子陣中挑選出符合要求的輸出波束，由13個(gè)子陣的對應(yīng)方向輸出組成新的信號矢量，其協(xié)方差矩陣表示為Ri，i=1，2…，M，再將M個(gè)協(xié)方差矩陣平均，得到平均協(xié)方差矩陣RM，對RM進(jìn)行特征分解即可得到多子陣波束域噪聲子空間和信號子空間。

由于矩陣運(yùn)算的相對獨(dú)立性，所以協(xié)方差矩陣可以采用并行計(jì)算的方式來實(shí)現(xiàn)。對于本系統(tǒng)而言，協(xié)方差矩陣Ri，為13×l3的Hermite陣，因此只需計(jì)算矩陣的上三角陣的元素即可。將所有上三角陣中的元素按照圖5中所示分成4組，分別由4片DSP進(jìn)行計(jì)算。

圖5 協(xié)方差矩陣計(jì)算分配

4.4 平均協(xié)方差矩陣特征分解及方位估計(jì)

對平均協(xié)方差矩陣進(jìn)行特征分解，采用前面介紹的并行Jacobi方法在FPGA2上實(shí)現(xiàn)。方位估計(jì)可以采用經(jīng)典的MUSIC算法空間搜索方法來實(shí)現(xiàn)，也可以采用Root-MUSIC算法求根實(shí)現(xiàn)。對于空間掃描方法而言，由于前端的波束形成已經(jīng)初步確定了回波的方位，因此不需要對整個(gè)空間進(jìn)行搜索，而只需對特定角度范圍進(jìn)行搜索即可，用兩片DSP實(shí)現(xiàn)。對于求根方法，采用兩片DSP輪流實(shí)現(xiàn)求根過程。

5 總結(jié)

從MSB-RMU算法的并行實(shí)現(xiàn)的角度入手，針對小規(guī)模系統(tǒng)設(shè)計(jì)的基于FPGA的雙邊CORDIC結(jié)構(gòu)的并行Jacobi實(shí)現(xiàn)方法，為特征分解類算法提供了一種高效的單芯片解決方案，對于類似算法的并行實(shí)現(xiàn)提供了一種思路。設(shè)計(jì)的由6片DSP和2片F(xiàn)PGA組成的并行信號處理平臺上較好的實(shí)現(xiàn)了MSB-RMU算法的實(shí)時(shí)處理，但對于大規(guī)模的并行處理平臺及算法的實(shí)時(shí)實(shí)現(xiàn)，如算法所用資源，工作頻率等因素，如何評估需要什么樣的硬件平臺，在FPGA上完成哪些功能等方面還需要進(jìn)行更深入的研究。

［1］周天.超寬覆蓋海底地形地貌高分辨探測技術(shù)研究［D］.哈爾濱工程大學(xué)，2005:43-50.

［2］周天，李海森，朱志德，等.多波束測深系統(tǒng)中多子陣檢測法的改進(jìn)及其性能分析［J］.聲學(xué)技術(shù)，2005，24（3）:138-143.

［3］任春麗，徐甲同，王俊平.實(shí)對稱三對角矩陣特征值的一種并行算法及實(shí)現(xiàn)［J］.西安電子科技大學(xué)學(xué)報(bào)，1999，26（2）:217-221.

［4］張麗君，金綏更.向量算法與并行算法［M］.北京:國防工業(yè)出版社，1993:123-130.

［5］袁生光.對稱矩陣特征值分解的硬件實(shí)現(xiàn)研究［D］.浙江大學(xué)，2008:30-37.

［6］MINSEOK KIM，KOICHI ICHIGE，HIROYUKI ARAI.Design of Jacobi EVD Processor based on CORDIC for DOA estimation with MUSIC algorithm ［C］.Personal，Indoor and Mobile Radio Communications，2002.The 13th IEEE International Symposium on，Lisboa，Portugal，2002:120-124

［7］李小波，薛王偉，孫志勇.一種求解復(fù)Hermite矩陣特征值的方法［J］.數(shù)據(jù)采集與處理，2005，20（4）: 403-406.

［8］謝倫國.并行計(jì)算機(jī)互聯(lián)網(wǎng)絡(luò)技術(shù)［M］.北京: 電子工業(yè)出版社，2004:154-169.