超分辨測(cè)向中奇異值分解的FPGA實(shí)現(xiàn)*

談文韜,朱業(yè)騰,黎仁剛,林 明

(1.江蘇科技大學(xué) 電子信息學(xué)院,江蘇 鎮(zhèn)江 212003;2.中船重工集團(tuán)723研究所,江蘇 揚(yáng)州 225101)

?

超分辨測(cè)向中奇異值分解的FPGA實(shí)現(xiàn)*

談文韜1,朱業(yè)騰2,黎仁剛2,林 明**1

(1.江蘇科技大學(xué) 電子信息學(xué)院,江蘇 鎮(zhèn)江 212003;2.中船重工集團(tuán)723研究所,江蘇 揚(yáng)州 225101)

奇異值分解是超分辨測(cè)向技術(shù)的核心組成部分，現(xiàn)有的并行實(shí)現(xiàn)方案適用范圍窄，運(yùn)算量大，迭代時(shí)間長(zhǎng)。為了滿足測(cè)向接收機(jī)系統(tǒng)的高實(shí)時(shí)性需求，結(jié)合雙邊Jacobi算法的交換策略和單邊Jacobi算法的求角結(jié)構(gòu)，提出了一種改進(jìn)的實(shí)現(xiàn)方法。該實(shí)現(xiàn)方法修正了脈動(dòng)陣列的收斂性問(wèn)題，提高了復(fù)數(shù)矩陣的收斂速度。同時(shí)，給出了算法的現(xiàn)場(chǎng)可編程門陣列(FPGA)實(shí)現(xiàn)結(jié)構(gòu)。仿真結(jié)果證明該方案耗時(shí)在百微秒以內(nèi)，能夠應(yīng)用于電子偵察設(shè)備。

電子偵察；超分辨測(cè)向；奇異值分解；脈動(dòng)陣列

1 引 言

由于現(xiàn)代戰(zhàn)場(chǎng)電磁環(huán)境的日益復(fù)雜，超分辨測(cè)向技術(shù)以其優(yōu)越的算法性能得到了越來(lái)越多的關(guān)注[1-2]，奇異值分解(Singular Value Decomposition,SVD)成為了超分辨測(cè)向技術(shù)實(shí)現(xiàn)與應(yīng)用的掣肘。與在其他方面的應(yīng)用不同[3-4]，應(yīng)用于電子對(duì)抗的奇異值分解[5]對(duì)精度敏感，具有高實(shí)時(shí)性、高可靠性的嚴(yán)格要求。因?yàn)榉治鰡?wèn)題角度不同，一般應(yīng)用于該領(lǐng)域的奇異值分解更注重實(shí)時(shí)性，運(yùn)行時(shí)間通常約束在100 μs以內(nèi)；而常規(guī)實(shí)現(xiàn)方法更注重資源消耗，運(yùn)行時(shí)間通常為1～100 ms級(jí)[6-7]，與電子偵察系統(tǒng)的支援應(yīng)用需求相距甚遠(yuǎn)。并行化脈動(dòng)(Brent-Lu-Van，BLV)陣列在20世紀(jì)80年代被應(yīng)用于SVD的計(jì)算[8]，大大縮短了超分辨測(cè)向算法的整體運(yùn)行時(shí)間。Bravo等人[9]于2006年對(duì)算法提出了改進(jìn)，但BLV陣列存在非對(duì)稱矩陣收斂慢甚至不收斂的情況，且僅適用偶數(shù)維度實(shí)數(shù)矩陣。在工程應(yīng)用中，有限位寬數(shù)據(jù)的舍入操作，導(dǎo)致計(jì)算過(guò)程中的矩陣數(shù)值并非完全對(duì)稱，傳統(tǒng)BLV結(jié)構(gòu)會(huì)放大該誤差的影響，影響算法的計(jì)算精度。

本文以多重分類(Multiple Signal Classification,MUSIC)算法中奇異值分解的應(yīng)用為出發(fā)點(diǎn)，在雙邊Jacobi算法BLV并行化時(shí)間結(jié)構(gòu)的基礎(chǔ)上，綜合考慮系統(tǒng)實(shí)時(shí)性、可靠性以及資源等需求，給出了一種適用超分辨測(cè)向算法的SVD并行化實(shí)現(xiàn)結(jié)構(gòu)。

2 MUSIC算法

假設(shè)M元任意天線陣列的所有陣元均位于坐標(biāo)系XOY平面內(nèi)，第k個(gè)陣元坐標(biāo)為(xk,yk,0)，第i個(gè)窄帶信號(hào)波長(zhǎng)為λi，來(lái)波方向?yàn)?θi,φi)，如圖1，則第k個(gè)陣元到圓心(即原點(diǎn))的相位差Δrik為

(1)

圖1 天線陣列Fig.1 Antenna array

天線陣列的接收公式為

(2)

x(t)=As(t)+n(t)。

(3)

式中：P表示信號(hào)源數(shù)目，s(t)表示P個(gè)入射信號(hào)的集合，x(t)表示M個(gè)陣元處接收信號(hào)的組合，A表示天線陣的陣列流形。

MUSIC算法基本原理為

Rxx=E[x(t)x(t)H]=ARssAH+σ2In=UsΣsVs+UnΣnVn，

(4)

(5)

(6)

式中：Rxx表示接收數(shù)據(jù)的協(xié)方差矩陣;Rss表示入射信號(hào)協(xié)方差矩陣;σ2表示噪聲功率;In表示單位矩陣;Us、Un分別為信號(hào)與噪聲協(xié)方差矩陣的左奇異矢量，也即信號(hào)子空間與噪聲子空間;Vs、Vn為右奇異矢量，對(duì)于方陣，同樣可作為子空間，通常只計(jì)算左奇異矢量;上標(biāo)“^”表示極大似然估計(jì);Pmusic為空間譜。

3 SVD并行實(shí)現(xiàn)算法

奇異值分解算法具有較為明確的物理意義，設(shè)Rxx為M陣元天線陣列接收信號(hào)的自相關(guān)矩陣，則Rxx為一個(gè)M×M的Hermite矩陣。因?yàn)镠ermite矩陣是正規(guī)矩陣，所以Rxx的特征向量中可找出一組正交基，且特征值均為實(shí)數(shù)；U與V均為由Rxx特征向量構(gòu)成的M×M酉矩陣，包含對(duì)應(yīng)信號(hào)在不同陣元間的相位差信息。同時(shí)，由于噪聲與計(jì)算誤差的存在，Rxx的特征值非零，為非奇異矩陣。

目前應(yīng)用最為廣泛的是雙邊Jacobi算法[10]與單邊Jacobi算法[11]。

雙邊Jacobi算法是一種通過(guò)平面旋轉(zhuǎn)完成消元計(jì)算，使矩陣收斂于由奇異值構(gòu)成的對(duì)角矩陣的可并行化算法。

以Rxx為例，雙邊Jacobi算法可以表示為

(7)

式中：S為矩陣Rxx對(duì)應(yīng)的奇異值對(duì)角矩陣，Gk為對(duì)應(yīng)m、n行列交點(diǎn)數(shù)值的Givens變換矩陣的一種變形。

應(yīng)用于FPGA的雙邊Jacobi算法并行化實(shí)現(xiàn)方法，是一種基于BLV脈動(dòng)陣列的實(shí)現(xiàn)結(jié)構(gòu)。該方法由兩種2×2的運(yùn)算單元構(gòu)成，分別為主對(duì)角線上角求取主運(yùn)算單元與其他位置的旋轉(zhuǎn)乘法單元。

在BLV陣列結(jié)構(gòu)中，每次旋轉(zhuǎn)后進(jìn)行相鄰行列的陣列交換，該交換策略的元素遍歷性較差，如圖3所示。圖中，U12與U21單元的b、c元素?zé)o法傳遞到主對(duì)角線上的旋轉(zhuǎn)單元內(nèi)，僅能通過(guò)傳遞的旋轉(zhuǎn)角進(jìn)行迭代消元，且計(jì)算誤差會(huì)隨著迭代的增加在類似位置的值內(nèi)累積，導(dǎo)致算法收斂慢，計(jì)算結(jié)果誤差較大。

圖2 傳統(tǒng)BLV陣列結(jié)構(gòu)Fig.2 Structure of traditional BLV array

單邊算法旋轉(zhuǎn)角度基于行列正交性計(jì)算得出，這與雙邊算法的旋轉(zhuǎn)消元方法有一定差異。本文引入了單邊Jacobi算法的旋轉(zhuǎn)角計(jì)算方法,其第m列與第n列元素間旋轉(zhuǎn)角的計(jì)算公式為

(8)

(9)

式中:Rcol(*)表示矩陣對(duì)應(yīng)列，sign(*)表示取符號(hào)。式(8)與式(9)對(duì)復(fù)數(shù)矩陣同樣適用。

4 實(shí)現(xiàn)方案選擇

針對(duì)雙邊Jacobi算法中BLV陣列只適用奇數(shù)維度實(shí)數(shù)矩陣、元素遍歷性差,以及單邊Jacobi算法收斂速度慢的問(wèn)題，結(jié)合兩種算法的計(jì)算特點(diǎn)，下面給出符合應(yīng)用背景的實(shí)現(xiàn)方案。

4.1 雙邊算法交換策略

傳統(tǒng)BLV陣列具有元素遍歷性差的特點(diǎn)。在運(yùn)算過(guò)程中，兩個(gè)主運(yùn)算單元副對(duì)角線上所有乘法單元內(nèi)的元素收斂速度慢，甚至出現(xiàn)發(fā)散的情況;而雙邊算法的行列交換策略具有最優(yōu)的元素遍歷性，算法收斂速度快，結(jié)果可靠。

兩兩組合主對(duì)角線上元素，共有M(M-1)/2種組合；每次向BLV中輸入?M/2」個(gè)不含重復(fù)列的組合，并按組合位置交換對(duì)應(yīng)行列；按此方式進(jìn)行一次計(jì)算，稱作一次向量正交化，完成所有正交化組合稱作“一掃”，完成一掃需要2?(M-1)/2」+1次并行向量正交化計(jì)算。將生成的交換策略轉(zhuǎn)變?yōu)榫仃嚨牡刂?，每個(gè)數(shù)據(jù)有行地址與列地址兩個(gè)坐標(biāo)信息，所以完成一掃所需的地址數(shù)據(jù)量L可由式(10)計(jì)算：

L=2M2(2?(M-1)/2」+1)。

(10)

每個(gè)地址的位寬為「lb(M+1)?，「*?表示向上取整。

組合的行列排列信息可以離線計(jì)算，并按矩陣元素地址格式量化，以常數(shù)的形式存儲(chǔ)，不會(huì)增加計(jì)算量；或者，通過(guò)設(shè)計(jì)一個(gè)地址發(fā)生器，進(jìn)行行列的遍歷。

4.2 復(fù)數(shù)計(jì)算方案

常規(guī)的復(fù)數(shù)矩陣分解方法是將維度為M的復(fù)數(shù)矩陣化為維度2M的實(shí)對(duì)稱矩陣，采用實(shí)數(shù)矩陣方法進(jìn)行分解;而奇異值分解的計(jì)算復(fù)雜度為O(M3)，算法的復(fù)雜度變?yōu)樵瓉?lái)的8倍以上，尤其是在需要計(jì)算特征向量時(shí)，大大增加了計(jì)算成本與實(shí)現(xiàn)難度。同時(shí)使用實(shí)數(shù)化處理時(shí)，增加了每一掃需要的向量正交化次數(shù)，也增加了算法的迭代時(shí)間。

此處采用單邊Jacobi算法的旋轉(zhuǎn)角計(jì)算方式，根據(jù)復(fù)數(shù)計(jì)算與現(xiàn)場(chǎng)可編程邏輯門陣列(Field Programmable Gate Array,FPGA)結(jié)構(gòu)對(duì)公式進(jìn)行變形優(yōu)化。將主運(yùn)算單元中的2×2矩陣進(jìn)行旋轉(zhuǎn)消元，需要滿足

(11)

因?yàn)镽xx為Hermite復(fù)數(shù)矩陣，所以a、d為實(shí)數(shù)，b、c相互共軛。根據(jù)式(8)和式(9)，復(fù)數(shù)旋轉(zhuǎn)角求取公式化簡(jiǎn)為

(12)

式中:“*′”表示復(fù)數(shù)的共軛。式(12)結(jié)構(gòu)更適于FPGA實(shí)現(xiàn)。同時(shí)不再向同行同列傳遞旋轉(zhuǎn)角，而是經(jīng)過(guò)復(fù)數(shù)修正的旋轉(zhuǎn)矩陣，即先乘以實(shí)數(shù)化矩陣的共軛轉(zhuǎn)置將矩陣還原：

(13)

4.3 奇數(shù)行(列)乘法單元

在實(shí)際運(yùn)算過(guò)程中，由于奇數(shù)末位行位沒(méi)有配對(duì)組合，不必計(jì)算旋轉(zhuǎn)角，但在旋轉(zhuǎn)相乘時(shí)，必須連接旋轉(zhuǎn)矩陣傳遞相乘，來(lái)保證行列變換的一致性。

本文通過(guò)增加單邊末行(列)單邊乘法運(yùn)算模塊，如圖3虛線框中部分，對(duì)奇數(shù)余行(列)的運(yùn)算進(jìn)行補(bǔ)充:

(14)

圖3 修正的BLV陣列Fig.3 Modified BLV array

最終的奇異值分解并行實(shí)現(xiàn)方案如圖2所示。該算法結(jié)構(gòu)中，誤差累積因素只有三角函數(shù)正、余弦匹配誤差，其他過(guò)程誤差僅影響收斂速度，與計(jì)算精度無(wú)關(guān)，保證了數(shù)值的穩(wěn)定性。

4.4 實(shí)現(xiàn)方案性能分析

此處分別對(duì)傳統(tǒng)BLV實(shí)現(xiàn)方案與本文實(shí)現(xiàn)方案進(jìn)行分析對(duì)比。

(15)

以9元天線陣接收數(shù)據(jù)實(shí)部的協(xié)方差矩陣作為BLV陣列的輸入數(shù)據(jù)，則ek隨正交化迭代次數(shù)的變化曲線如圖4和圖5所示。

圖4 傳統(tǒng)BLV方案的曲線Fig.4 Curve of traditional BLV solution

圖5 本文方案的曲線Fig.5 Curve of the proposed solution

從圖4中可以看出，即便算法的輸入數(shù)據(jù)維度M相對(duì)較小(M<10)，BLV陣列的部分副對(duì)角線上的值無(wú)法消元，導(dǎo)致奇異值的均方誤差不能收斂于0值，與第3節(jié)中的描述一致,因此,該實(shí)現(xiàn)方法無(wú)法滿足超分辨測(cè)向算法的應(yīng)用需求。

以9元天線陣接收復(fù)數(shù)矩陣作為本文算法的輸入數(shù)據(jù)，在圖5中，可以發(fā)現(xiàn)，均方誤差值的下降曲線較為平滑，且在放大框圖中可以看出，到第3掃(第27次正交變換)完成時(shí)，均方誤差已經(jīng)接近10-5。

試驗(yàn)結(jié)果說(shuō)明了本文的方案選擇，收斂速度更快，穩(wěn)健性更好。該性能分析的過(guò)程中部分?jǐn)?shù)據(jù)為雙精度浮點(diǎn)數(shù)，比起FPGA內(nèi)單精度塊浮點(diǎn)的實(shí)現(xiàn)方法存在可控誤差。

5 硬件實(shí)現(xiàn)方案

傳統(tǒng)BLV陣列的主運(yùn)算模塊至少需要3個(gè)坐標(biāo)旋轉(zhuǎn)計(jì)算法(Coordinate Rotation Digital Computer,CORDIC)模塊，CORDIC算法模塊周期長(zhǎng)，大大增加了算法輸出的延時(shí)。此處使用組合邏輯與查找表方式分別設(shè)計(jì)了2個(gè)CORDIC模塊，大大縮減了旋轉(zhuǎn)角求取時(shí)間。旋轉(zhuǎn)角與sp矩陣計(jì)算結(jié)構(gòu)相似，求取順序不同，通過(guò)一次復(fù)用減少了算法的資源消耗。旋轉(zhuǎn)矩陣求取模塊實(shí)現(xiàn)結(jié)構(gòu)如圖6所示。

圖6 旋轉(zhuǎn)角求取模塊結(jié)構(gòu)Fig.6 Structure of rotation angle calculation module

基于雙邊算法交換策略的奇異值分解迭代流程如圖7所示。

圖7 奇異值分解迭代流程Fig.7 Iterative process of SVD

本文以200 MHz時(shí)鐘、塊浮點(diǎn)數(shù)據(jù)結(jié)構(gòu)(符號(hào)位1,指數(shù)位寬8,溢出保護(hù)位1，數(shù)據(jù)位寬24)進(jìn)行算法硬件實(shí)現(xiàn)。每個(gè)旋轉(zhuǎn)角計(jì)算單元需要 7 421個(gè)查找表(Look-up Table,LUT)、790個(gè)觸發(fā)器(Flip-flop,FF)、16個(gè)DSP48E1(一種FPGA乘法器)，每個(gè)乘法單元需要858個(gè)LUT、438個(gè)FF、40個(gè)DSP48E1。

以ε=10-6為向量正交精度，完成一次9行9列的奇異值分解計(jì)算僅需3 172個(gè)周期，約15.86 μs。驗(yàn)證結(jié)果證明該實(shí)現(xiàn)方案優(yōu)于文獻(xiàn)[4-5]中的毫秒級(jí)實(shí)現(xiàn)方案。

6 結(jié) 論

本文以雙邊Jacobi算法的BLV實(shí)現(xiàn)方案為基礎(chǔ)，結(jié)合單邊算法與奇數(shù)維度矩陣實(shí)現(xiàn)方法的優(yōu)點(diǎn)，完成了一種遍歷性、穩(wěn)健性、收斂速度均較優(yōu)的實(shí)現(xiàn)方案。該方案耗時(shí)被約束在百μs以內(nèi)，遠(yuǎn)低于原有的毫秒級(jí)實(shí)現(xiàn)方案，使超分辨測(cè)向算法在高實(shí)時(shí)性電子偵察系統(tǒng)中的應(yīng)用成為可能。本文的實(shí)現(xiàn)方法消耗的片上資源較多，仍有進(jìn)一步優(yōu)化的空間。

[1] 胡子揚(yáng),任淵.一種最小冗余線陣的目標(biāo)DOA估計(jì)方法[J].電訊技術(shù),2014,54(11):1493-1498. HU Ziyang,REN Yuan. A DOA estimation method based on sub-minimum redundancy linear array[J].Telecommunication Engineering,2014,54(11):1493-1498.(in Chinese)

[2] 付淑娟,景小榮,張祖凡,等.基于虛擬陣列改進(jìn)MUSIC算法的相干信源DoA估計(jì)[J].電訊技術(shù),2011,51(11):63-67. FU Shujuan,JING Xiaorong,ZHANG Zufan,et al.DoA estimation of coherent sources by using virtual array-based improved MUSIC algorithm[J].Telecommunication Engineering,2011,51(11):63-67.(in Chinese)

[3] MOHANTY R,ANIRUDH G,PRADHAN T. Design and performance analysis of fixed-point Jacobi SVD algorithm on reconfigurable system[J].IERI Procedia,2014(7):21-27.

[4] QIAO H L. New SVD based initialization strategy for non-negative matrix factorization[J].Pattern Recognition Letters,2015,63(C):71-77.

[5] LIU Y,CUI H Y. Antenna array signal direction of arrival estimation on digital signal processor(DSP)[J].Procedia Computer Science,2015,55(7):782-791.

[6] MILFORD D,SANDELL M. Singular value decomposition using an array of CORDIC processors[J].Signal Processing,2014,102(9):163-170.

[7] TAI Y G,DAN C T,PSARRIS K. Scalable matrix decompositions with multiple cores on FPGAs[J].Microprocessors and Microsystems,2013,37(8):887-898.

[8] LUK F T,BRENT R P. The solution of singular-value and symmetric eigenvalue problems on multiprocessor arrays[J].SIAM Journal on Scientific and Statistical Computing,1985,6(1)：69-84．

[9] BRAVO I,JIMNEZ P,MAZO M,et al.Implementation in FPGAs of Jacobi method to solve the eigenvalue and eigenvector problem[C]//Proceedings of 2006 International Conference on Field Programmable Logic & Applications.Madrid,Spain:IEEE,2006:1-4.

[10] JAMES D,KRESIMIR V. Jacobi's method is more accurate than QR[J].SIAM Journal on Matrix Analysis and Applications,1992,13(4):1204-1245.

[11] 郭強(qiáng),趙雷. 一種基于動(dòng)態(tài)序列的單邊Jacobi方法[J].蘇州大學(xué)學(xué)報(bào)(工科版),2011,31(4):16-22. GUO Qiang,ZHAO Lei. A new one-side Jacobi basedon dynamic ordering[J].Journal of Soochow University(Engineering Science Edition),2011,31(4):16-22.(in Chinese)

FPGA Implementation of Singular Value Decomposition Algorithm in Super-resolution Direction-finding

TAN Wentao1,ZHU Yeteng2,LI Rengang2,LIN Ming1

(1.School of Electronics Information,Jiangsu University of Science and Technology,Zhenjiang 212003,China;2.No.723 Research Institute of China Shipbuilding Industry Corporation,Yangzhou 225101,China)

Singular value decomposition(SVD) is the core of super-resolution direction-finding(DF) technique.The existing parallel implementation methods have narrow applicability,intensive computation and lengthy iteration. In order to meet the application requirements of DF receiver system,a modified implementation method is proposed according to the exchange policy of bilateral Jacobi algorithm and the angle calculation structure of unilateral Jacobi algorithm. In this implementation method,the ergodic problem in Brent-Lu-Van(BLV) array is modified,and the convergence rate in complex matrix is improved.An algorithm structure on field programmable gate array(FPGA) is presented.Simulation results show that the method consumes less than 100 ms and can be applied to electronic reconnaissance equipment.

electronic reconnaissance;super-resolution direction-finding;sigular value decomposition(SVD);Brent-Lu-Van(BLV) array

10.3969/j.issn.1001-893x.2017.07.012

談文韜,朱業(yè)騰,黎仁剛,等.超分辨測(cè)向中奇異值分解的FPGA實(shí)現(xiàn)[J].電訊技術(shù)，2017,57(7):801-805.[TAN Wentao,ZHU Yeteng,LI Rengang,et al.FPGA implementation of singular value decomposition algorithm in super-resolution direction-finding[J].Telecommunication Engineering,2017,57(7):801-805.]

2016-10-26；

2017-02-22 Received date:2016-10-26；Revised date：2017-02-22

國(guó)家自然科學(xué)基金青年科學(xué)基金項(xiàng)目(61401179)

TN971

A

1001-893X(2017)07-0801-05

談文韜(1991—),男,江蘇連云港人,碩士研究生,主要研究方向?yàn)槔走_(dá)信號(hào)與信息處理理論與技術(shù);

Email:DorusTan@163.com

朱業(yè)騰(1988—),男,江蘇揚(yáng)州人,工程師,主要研究方向?yàn)閿?shù)字信號(hào)處理;

黎仁剛(1978—),男,江蘇揚(yáng)州人,博士,研究員,主要研究方向?yàn)殛嚵行盘?hào)處理;

林 明(1960—),男,遼寧大連人,教授、碩士生導(dǎo)師,主要研究方向?yàn)槔走_(dá)信號(hào)處理。

Email:jskdlm@qq.com

**通信作者：jskdlm@qq.com Corresponding author：jskdlm@qq.com