陳 靜，祝嬌嬌，劉盛典，劉 鵬，徐振興

(航天恒星科技有限公司，北京 102102)

0 引言

在電子信息快速發(fā)展的當(dāng)今社會(huì)，衛(wèi)星導(dǎo)航在軍民應(yīng)用領(lǐng)域占據(jù)越來越重要的地位，干擾和抗干擾技術(shù)的發(fā)展直接影響導(dǎo)航應(yīng)用的有效精度[1]。而在抗干擾的過程中，對(duì)空間信號(hào)波達(dá)方向(Direction of Arrival,DOA)估計(jì)、干擾個(gè)數(shù)檢測(cè)、最優(yōu)權(quán)矢量的實(shí)現(xiàn)精度影響著導(dǎo)航接收機(jī)的抗干擾性能，協(xié)方差矩陣的特征分解是這些算法實(shí)現(xiàn)過程中的核心部分[2-4]。

本文根據(jù)自適應(yīng)陣列天線獲得的協(xié)方差矩陣性能，將復(fù)數(shù)陣列協(xié)方差矩陣的特征分解進(jìn)行現(xiàn)場(chǎng)可編程門陣列(Field Programmable Gate Array，F(xiàn)PGA)實(shí)現(xiàn)，并通過將其應(yīng)用于空間信號(hào)DOA估計(jì)進(jìn)行驗(yàn)證。另外，在實(shí)現(xiàn)的過程中對(duì)直接調(diào)用CORDIC IP核的方式進(jìn)行了精度誤差分析，并用一種雙精度浮點(diǎn)的方式進(jìn)行修正，提高了FPGA矩陣特征分解實(shí)現(xiàn)精度。

1 陣列協(xié)方差矩陣特征分解

在抗干擾處理過程中，通常得到的數(shù)據(jù)是有限時(shí)間范圍內(nèi)的有限次快拍數(shù)。在這段時(shí)間內(nèi)假定接收到的導(dǎo)航信號(hào)的方向不發(fā)生變化，或者認(rèn)為導(dǎo)航信號(hào)的包絡(luò)雖然隨時(shí)間變化，但它是一個(gè)平穩(wěn)隨機(jī)的過程，其統(tǒng)計(jì)特性不隨時(shí)間變化。而且天線的陣元個(gè)數(shù)大于干擾信號(hào)個(gè)數(shù)，干擾信號(hào)的方向向量是線性獨(dú)立的，陣列中噪聲具有高分布特性，所以接收機(jī)接收到的信號(hào)向量的協(xié)方差矩陣是對(duì)角非奇異陣，也是復(fù)數(shù)正定Hermite陣。

1.1 復(fù)數(shù)協(xié)方差矩陣轉(zhuǎn)化為實(shí)矩陣

相對(duì)于實(shí)數(shù)計(jì)算，復(fù)數(shù)的處理增加了FPGA運(yùn)算的復(fù)雜度，所以為了降低計(jì)算復(fù)雜度，在設(shè)計(jì)的過程中需將復(fù)數(shù)協(xié)方差矩陣轉(zhuǎn)化為實(shí)數(shù)矩陣。

根據(jù)文獻(xiàn)[5]將n階協(xié)方差矩陣Rxx寫成實(shí)數(shù)矩陣和虛數(shù)矩陣兩部分，即

Rxx=A+iB

(1)

式中，A、B均為實(shí)矩陣。

由于協(xié)方差矩陣是復(fù)數(shù)正定Hermite陣，根據(jù)其性能可知

AT=A，BT=-B

(2)

設(shè)特征值為λ，對(duì)應(yīng)的特征向量為U+iV，其中U、V均為n維實(shí)列向量，則有

(A+iB)(U+iV)=λ(U+iV)

(3)

根據(jù)式(3)可寫為

(4)

1.2 實(shí)數(shù)矩陣特征分解

經(jīng)典的實(shí)數(shù)矩陣特征值計(jì)算算法中，雅克比(Jacobi)算法在同等計(jì)算精度要求下具有更快的速度[6-7]，且具有并行計(jì)算的優(yōu)勢(shì)，可將分解過程劃分為并行的子步驟同時(shí)進(jìn)行。另外根據(jù)應(yīng)用需求，需得到特征值對(duì)應(yīng)的特征向量，而雙邊Jacobi的計(jì)算結(jié)果可以直接給出特征值對(duì)應(yīng)的特征向量，所以采用雙邊Jacobi算法進(jìn)行特征值和特征向量的FPGA實(shí)現(xiàn)。

Jacobi算法通過對(duì)實(shí)數(shù)矩陣S進(jìn)行一系列平面旋轉(zhuǎn)變換來產(chǎn)生一系列對(duì)稱方陣Ak，每次旋轉(zhuǎn)變換使Ak中的2個(gè)元素變?yōu)?[8]。當(dāng)多次迭代后，Ak趨于一個(gè)對(duì)角陣D，D的各主對(duì)角元素就是S的特征值[9-10]。從Ak到Ak+1的旋轉(zhuǎn)變換公式為

(5)

(6)

T的選擇為

(7)

式中，p

(8)

進(jìn)行一次旋轉(zhuǎn)變換，矩陣Ak相對(duì)于矩陣Ak-1變化了的元素是

(9)

(10)

(11)

(12)

Jacobi算法在每次進(jìn)行旋轉(zhuǎn)變換后，根據(jù)式(5)、式(7)可以看出，實(shí)對(duì)稱矩陣的第p行、q列、q行、p列都發(fā)生了變化，變化結(jié)果如式(10)、式(11)所示，而矩陣的其他行、列的元素與變換前相應(yīng)的元素相同，并不會(huì)發(fā)生變化。

利用這一特性，對(duì)于n×n的矩陣，在FPGA設(shè)計(jì)時(shí)可采用n/2個(gè)處理器并行處理的方法同時(shí)消去2n個(gè)非對(duì)角元素，一次迭代只需n-1步完成，減少了FPGA的運(yùn)算時(shí)間。所以對(duì)于陣元數(shù)為M的陣列天線，協(xié)方差矩陣轉(zhuǎn)化為2M×2M的實(shí)矩陣后，可采用M個(gè)處理器并行處理。但每次旋轉(zhuǎn)后矩陣需按表1所示的規(guī)律進(jìn)行行列交換，以確保新矩陣的正確性，一次迭代需進(jìn)行2M-1次旋轉(zhuǎn)和2M-1次行列交換。

表1 一次迭代處理器循環(huán)次數(shù)及行列交換規(guī)律

Jacobi算法求解實(shí)數(shù)矩陣S特征值和特征向量的具體步驟為：

1)初始化特征向量為單位陣E，即主對(duì)角線元素為1，其他元素為0；

2)根據(jù)式(12)計(jì)算tan2θ，求sinθ、cosθ及矩陣Tpq；

3)利用式(5)～式(6)計(jì)算得到A1，用當(dāng)前特征向量矩陣E乘以矩陣Tpq得到特征向量V1；

5)若當(dāng)前迭代前的矩陣A的非對(duì)角線元素中最大值小于給定的閾值e時(shí)，停止計(jì)算；否則，令A(yù)=A′1，V=V′1，重復(fù)執(zhí)行步驟2)～5)。停止計(jì)算時(shí)，得到特征值Ii≈(A1)ij,i,j=1,2,…,n以及特征向量V。

6)根據(jù)特征值從大到小的順序重新排列矩陣的特征值和特征向量，對(duì)復(fù)數(shù)協(xié)方差矩陣的特征值/特征向量進(jìn)行提取。

2 特征求解FPGA實(shí)現(xiàn)設(shè)計(jì)

自適應(yīng)陣列天線接收的衛(wèi)星信號(hào)，經(jīng)累加后得到協(xié)方差矩陣進(jìn)行特征分解，F(xiàn)PGA實(shí)現(xiàn)流程圖如圖1所示。首先將協(xié)方差矩陣輸入change_rxx2r模塊，使得M×M的復(fù)數(shù)Hermite協(xié)方差矩陣轉(zhuǎn)化為2M×2M實(shí)矩陣，M為天線陣元數(shù)；對(duì)于得到的2M×2M的實(shí)矩陣，一方面根據(jù)第1節(jié)中特征求解的步驟2)，進(jìn)行旋轉(zhuǎn)角度的正余弦值求解，另一方面存入ram中實(shí)現(xiàn)循環(huán)迭代和矩陣元素的讀取。接下來，通過M個(gè)并行處理器compute實(shí)現(xiàn)步驟3)，對(duì)輸入矩陣進(jìn)行行、列2次角度旋轉(zhuǎn)，然后進(jìn)入exchange_value模塊實(shí)現(xiàn)步驟4)，完成行、列交換，其中陣元數(shù)不同，一次迭代(sweep)交換的次數(shù)和規(guī)律不同(參考表1)。每次sweep結(jié)束后，通過find_max模塊尋找新矩陣非對(duì)角元素的最大值ξ，進(jìn)行閾值判斷。ξ大于閾值時(shí)，將sweep后的新矩陣存入ram中進(jìn)行下一次sweep迭代運(yùn)算；ξ小于等于閾值時(shí)，矩陣的對(duì)角線上元素即為實(shí)矩陣特征值。最后，將計(jì)算的特征值由大到小排序，去掉重根，根據(jù)式(4)得到復(fù)數(shù)協(xié)方差矩陣的特征值,模塊sort_self完成該過程實(shí)現(xiàn)。

特征向量的實(shí)現(xiàn)，依賴于特征值實(shí)現(xiàn)過程中得到的正余弦值。當(dāng)協(xié)方差矩陣進(jìn)入特征分解模塊時(shí)eig_vector_start=1，特征向量求解開始，初始化特征向量矩陣為2M×2M的單位矩陣。根據(jù)第1節(jié)步驟3)，在eigvector模塊中對(duì)單位矩陣進(jìn)行1次角度旋轉(zhuǎn)，將旋轉(zhuǎn)矩陣存入ram中，根據(jù)讀地址的順序?qū)πD(zhuǎn)矩陣進(jìn)行行列交換(參考表1)。最后，根據(jù)特征值由大到小的排列順序，將對(duì)應(yīng)特征向量進(jìn)行排序，并根據(jù)式(4)將其轉(zhuǎn)化為復(fù)數(shù)協(xié)方差矩陣對(duì)應(yīng)的特征向量，模塊sort_self完成該過程實(shí)現(xiàn)。

圖1 矩陣特征求解的FPGA實(shí)現(xiàn)流程圖Fig.1 The flow chart of matrix eigendecomposition FPGA implementation

3 仿真分析及應(yīng)用

根據(jù)式(10)、式(11)可以看出，矩陣的每次旋轉(zhuǎn)需計(jì)算旋轉(zhuǎn)角度的正余弦值，在FPGA實(shí)現(xiàn)過程中分別用兩種方法進(jìn)行計(jì)算。

3.1 CORDIC算法實(shí)現(xiàn)仿真

CORDIC算法是一系列與運(yùn)算基數(shù)相關(guān)的角度不斷偏擺，從而逼近所需旋轉(zhuǎn)的角度，可通過該算法進(jìn)行向量旋轉(zhuǎn)、三角函數(shù)、乘、除等運(yùn)算[11-12]。CORDIC算法分為旋轉(zhuǎn)模式和矢量模式，兩種模式的輸入輸出方式如表2所示，其中X0、Y0、Z0為輸入初值。

表2 CORDIC算法兩種模式下輸出結(jié)果

根據(jù)式(10)～式(12)所示，在特征求解的過程中需要進(jìn)行2次旋轉(zhuǎn)模式運(yùn)算和1次矢量模式運(yùn)算，在FPGA實(shí)現(xiàn)時(shí)可通過調(diào)用CORDIC IP核來實(shí)現(xiàn)矩陣的旋轉(zhuǎn)[13]。

對(duì)CORDIC算法實(shí)現(xiàn)的特征分解FPGA工程進(jìn)行仿真分析。仿真設(shè)置陣元數(shù)為7的陣列天線接收信號(hào)，進(jìn)行單寬帶干擾，干擾和信號(hào)的功率比為80dB。經(jīng)計(jì)算協(xié)方差矩陣的特征分解情況如表3所示。根據(jù)對(duì)比可以看出，F(xiàn)PGA定點(diǎn)求出的特征值，在小特征值時(shí)會(huì)有200左右的誤差，從而使部分小特征值不能區(qū)分重根，影響協(xié)方差矩陣特征值和特征向量的提取。

表3 基于COEDIC算法FPGA實(shí)現(xiàn)和MATLAB自實(shí)現(xiàn)特征值數(shù)據(jù)比較

進(jìn)一步對(duì)FPGA工程分析可以看出，定點(diǎn)特征分解的誤差主要集中在CORDIC核產(chǎn)生的旋轉(zhuǎn)角度θ，由于每次生成θ都會(huì)產(chǎn)生誤差，隨著生成次數(shù)和θ的減小，誤差會(huì)逐漸增大，如圖2所示。

圖2 旋轉(zhuǎn)角度隨旋轉(zhuǎn)次數(shù)的變化Fig.2 Variation of rotation angle with rotation times

由于CORDIC核輸入/輸出數(shù)據(jù)位寬最大為48bit，精度位寬為48bit。在抗干擾過程中，為了提高抗干擾性能，模擬信號(hào)轉(zhuǎn)換為16bit位寬的數(shù)字信號(hào)，經(jīng)FPGA處理協(xié)方差矩陣數(shù)據(jù)位寬為18bit。為提高特征分解的精度，需擴(kuò)大特征分解模塊的輸入輸出位寬，但由于CORDIC核的計(jì)算精度為:

P=輸出位寬+輸入位寬+log2(輸出位寬)

(13)

則特征分解的輸入/輸出位寬最大可擴(kuò)到21bit，經(jīng)FPGA仿真實(shí)現(xiàn)的特征值精度較差，不能區(qū)分小特征值，如表3所示。若提高特征值分解精度需擴(kuò)大輸入位寬，直接調(diào)用CORDIC IP核的實(shí)現(xiàn)方法不能滿足，可以采用雙精度浮點(diǎn)型的方式進(jìn)行FPGA實(shí)現(xiàn)，以滿足大位寬、高精度特征分解需求。

3.2 雙精度浮點(diǎn)實(shí)現(xiàn)仿真

根據(jù)分析可知調(diào)用CORDIC IP核產(chǎn)生的誤差是不可避免的，但可以通過對(duì)式(12)進(jìn)一步的推導(dǎo)，實(shí)現(xiàn)雙精度浮點(diǎn)正余弦值的求解，以避免CORDIC核的使用，從而消除誤差。

根據(jù)式

(14)

(15)

在圖1所示的FPGA實(shí)現(xiàn)流程中，cos_sin模塊通過式(14)、式(15)進(jìn)行雙精度浮點(diǎn)型正余弦值計(jì)算，通過實(shí)際采集的衛(wèi)星信號(hào)累加計(jì)算得到的協(xié)方差矩陣進(jìn)行仿真分析。在實(shí)際采集的信號(hào)中，單寬帶干擾和信號(hào)的功率比為80dB，接收七陣元陣列天線，經(jīng)過雙精度浮點(diǎn)特征求解后得到的特征值如表4所示，與MATLAB自實(shí)現(xiàn)的特征值完全相同，滿足精度要求，表5所示為對(duì)應(yīng)求得的特征向量。

表4 基于雙精度浮點(diǎn)FPGA實(shí)現(xiàn)和MATLAB自實(shí)現(xiàn)特征值數(shù)據(jù)比較

表5 基于雙精度浮點(diǎn)FPGA實(shí)現(xiàn)的特征值對(duì)應(yīng)的特征向量

協(xié)方差矩陣特征分解一個(gè)重要的應(yīng)用是空間信號(hào)DOA估計(jì)，對(duì)FPGA實(shí)現(xiàn)的特征向量可以通過經(jīng)典DOA估計(jì)-多重信號(hào)分類(Multiple Signal Classification, MUSIC)進(jìn)行驗(yàn)證[14]。

MUSIC算法[15]主要是通過尋求陣列空間譜函數(shù)PMUSIC的峰值來得到DOA的估計(jì)。

(16)

其中，UN為協(xié)方差矩陣小特征值對(duì)應(yīng)的特征向量，即噪聲子空間；a(θ)為導(dǎo)向矢量。

仿真設(shè)置為七陣元陣列天線接收信號(hào)，單窄帶干擾，干擾與接收信號(hào)的功率比值為80dB，干擾的俯仰角設(shè)為20°，方位角設(shè)為160°。經(jīng)FPGA進(jìn)行協(xié)方差特征求解的特征值和特征向量進(jìn)行MUSIC估計(jì)，結(jié)果如圖3所示，估計(jì)的俯仰角與方位角與設(shè)定值完全相同，基于FPGA實(shí)現(xiàn)的特征分解精度符合DOA估計(jì)應(yīng)用。

圖3 FPGA實(shí)現(xiàn)特征分解DOA估計(jì)Fig.3 DOA estimation for FPGA implementation of eigendecomposition

4 結(jié)論

導(dǎo)航抗干擾接收機(jī)接收到的信號(hào)，其協(xié)方差矩陣是復(fù)數(shù)正定Hermite陣。在特征分解的FPGA實(shí)現(xiàn)的過程中，為了簡(jiǎn)化實(shí)現(xiàn)過程，需將復(fù)數(shù)矩陣轉(zhuǎn)化為實(shí)數(shù)矩陣。

通過雙邊Jacobi算法對(duì)實(shí)數(shù)矩陣進(jìn)行特征值和特征向量的求解，在FPGA設(shè)計(jì)時(shí)，每次旋轉(zhuǎn)角度正余弦的計(jì)算誤差直接影響實(shí)現(xiàn)精度。分別用CORDIC算法和雙精度浮點(diǎn)直接計(jì)算方法對(duì)旋轉(zhuǎn)角度正余弦進(jìn)行計(jì)算，通過對(duì)比仿真，CORDIC IP核實(shí)現(xiàn)方式具有一定的局限性，適用于處理數(shù)據(jù)較小矩陣，對(duì)于目前使用的導(dǎo)航接收機(jī)，信號(hào)的協(xié)方差矩陣數(shù)據(jù)為18bit，處理精度不能滿足要求。經(jīng)仿真分析，雙精度浮點(diǎn)直接計(jì)算方法的FPGA實(shí)現(xiàn)滿足特征分解精度要求，并通過DOA估計(jì)進(jìn)行應(yīng)用驗(yàn)證，DOA估計(jì)結(jié)果正確，滿足應(yīng)用需求?；贔PGA的特征分解實(shí)現(xiàn)的方法設(shè)計(jì)和仿真分析，為導(dǎo)航接收機(jī)抗干擾性能的提升提供了有效的工程基礎(chǔ)。