基于FPGA的自適應陣列協方差矩陣特征分解的實現

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

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

0 引言

在電子信息快速發展的當今社會，衛星導航在軍民應用領域占據越來越重要的地位，干擾和抗干擾技術的發展直接影響導航應用的有效精度[1]。而在抗干擾的過程中，對空間信號波達方向(Direction of Arrival,DOA)估計、干擾個數檢測、最優權矢量的實現精度影響著導航接收機的抗干擾性能，協方差矩陣的特征分解是這些算法實現過程中的核心部分[2-4]。

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

1 陣列協方差矩陣特征分解

在抗干擾處理過程中，通常得到的數據是有限時間范圍內的有限次快拍數。在這段時間內假定接收到的導航信號的方向不發生變化，或者認為導航信號的包絡雖然隨時間變化，但它是一個平穩隨機的過程，其統計特性不隨時間變化。而且天線的陣元個數大于干擾信號個數，干擾信號的方向向量是線性獨立的，陣列中噪聲具有高分布特性，所以接收機接收到的信號向量的協方差矩陣是對角非奇異陣，也是復數正定Hermite陣。

1.1 復數協方差矩陣轉化為實矩陣

相對于實數計算，復數的處理增加了FPGA運算的復雜度，所以為了降低計算復雜度，在設計的過程中需將復數協方差矩陣轉化為實數矩陣。

根據文獻[5]將n階協方差矩陣Rxx寫成實數矩陣和虛數矩陣兩部分，即

Rxx=A+iB

(1)

式中，A、B均為實矩陣。

由于協方差矩陣是復數正定Hermite陣，根據其性能可知

AT=A，BT=-B

(2)

設特征值為λ，對應的特征向量為U+iV，其中U、V均為n維實列向量，則有

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

(3)

根據式(3)可寫為

(4)

1.2 實數矩陣特征分解

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

Jacobi算法通過對實數矩陣S進行一系列平面旋轉變換來產生一系列對稱方陣Ak，每次旋轉變換使Ak中的2個元素變為0[8]。當多次迭代后，Ak趨于一個對角陣D，D的各主對角元素就是S的特征值[9-10]。從Ak到Ak+1的旋轉變換公式為

(5)

(6)

T的選擇為

(7)

式中，p

(8)

進行一次旋轉變換，矩陣Ak相對于矩陣Ak-1變化了的元素是

(9)

(10)

(11)

(12)

Jacobi算法在每次進行旋轉變換后，根據式(5)、式(7)可以看出，實對稱矩陣的第p行、q列、q行、p列都發生了變化，變化結果如式(10)、式(11)所示，而矩陣的其他行、列的元素與變換前相應的元素相同，并不會發生變化。

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

表1 一次迭代處理器循環次數及行列交換規律

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

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

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

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

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

6)根據特征值從大到小的順序重新排列矩陣的特征值和特征向量，對復數協方差矩陣的特征值/特征向量進行提取。

2 特征求解FPGA實現設計

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

特征向量的實現，依賴于特征值實現過程中得到的正余弦值。當協方差矩陣進入特征分解模塊時eig_vector_start=1，特征向量求解開始，初始化特征向量矩陣為2M×2M的單位矩陣。根據第1節步驟3)，在eigvector模塊中對單位矩陣進行1次角度旋轉，將旋轉矩陣存入ram中，根據讀地址的順序對旋轉矩陣進行行列交換(參考表1)。最后，根據特征值由大到小的排列順序，將對應特征向量進行排序，并根據式(4)將其轉化為復數協方差矩陣對應的特征向量，模塊sort_self完成該過程實現。

圖1 矩陣特征求解的FPGA實現流程圖Fig.1 The flow chart of matrix eigendecomposition FPGA implementation

3 仿真分析及應用

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

3.1 CORDIC算法實現仿真

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

表2 CORDIC算法兩種模式下輸出結果

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

對CORDIC算法實現的特征分解FPGA工程進行仿真分析。仿真設置陣元數為7的陣列天線接收信號，進行單寬帶干擾，干擾和信號的功率比為80dB。經計算協方差矩陣的特征分解情況如表3所示。根據對比可以看出，FPGA定點求出的特征值，在小特征值時會有200左右的誤差，從而使部分小特征值不能區分重根，影響協方差矩陣特征值和特征向量的提取。

表3 基于COEDIC算法FPGA實現和MATLAB自實現特征值數據比較

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

圖2 旋轉角度隨旋轉次數的變化Fig.2 Variation of rotation angle with rotation times

由于CORDIC核輸入/輸出數據位寬最大為48bit，精度位寬為48bit。在抗干擾過程中，為了提高抗干擾性能，模擬信號轉換為16bit位寬的數字信號，經FPGA處理協方差矩陣數據位寬為18bit。為提高特征分解的精度，需擴大特征分解模塊的輸入輸出位寬，但由于CORDIC核的計算精度為:

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

(13)

則特征分解的輸入/輸出位寬最大可擴到21bit，經FPGA仿真實現的特征值精度較差，不能區分小特征值，如表3所示。若提高特征值分解精度需擴大輸入位寬，直接調用CORDIC IP核的實現方法不能滿足，可以采用雙精度浮點型的方式進行FPGA實現，以滿足大位寬、高精度特征分解需求。

3.2 雙精度浮點實現仿真

根據分析可知調用CORDIC IP核產生的誤差是不可避免的，但可以通過對式(12)進一步的推導，實現雙精度浮點正余弦值的求解，以避免CORDIC核的使用，從而消除誤差。

根據式

(14)

(15)

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

表4 基于雙精度浮點FPGA實現和MATLAB自實現特征值數據比較

表5 基于雙精度浮點FPGA實現的特征值對應的特征向量

協方差矩陣特征分解一個重要的應用是空間信號DOA估計，對FPGA實現的特征向量可以通過經典DOA估計-多重信號分類(Multiple Signal Classification, MUSIC)進行驗證[14]。

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

(16)

其中，UN為協方差矩陣小特征值對應的特征向量，即噪聲子空間；a(θ)為導向矢量。

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

圖3 FPGA實現特征分解DOA估計Fig.3 DOA estimation for FPGA implementation of eigendecomposition

4 結論

導航抗干擾接收機接收到的信號，其協方差矩陣是復數正定Hermite陣。在特征分解的FPGA實現的過程中，為了簡化實現過程，需將復數矩陣轉化為實數矩陣。

通過雙邊Jacobi算法對實數矩陣進行特征值和特征向量的求解，在FPGA設計時，每次旋轉角度正余弦的計算誤差直接影響實現精度。分別用CORDIC算法和雙精度浮點直接計算方法對旋轉角度正余弦進行計算，通過對比仿真，CORDIC IP核實現方式具有一定的局限性，適用于處理數據較小矩陣，對于目前使用的導航接收機，信號的協方差矩陣數據為18bit，處理精度不能滿足要求。經仿真分析，雙精度浮點直接計算方法的FPGA實現滿足特征分解精度要求，并通過DOA估計進行應用驗證，DOA估計結果正確，滿足應用需求。基于FPGA的特征分解實現的方法設計和仿真分析，為導航接收機抗干擾性能的提升提供了有效的工程基礎。