基于FPGA的復數浮點協方差矩陣實現

摘 要:基于FPGA的協方差矩陣運算的實現大多采用的是定點計算方式，在運算過程中，存在數據處理動態范圍小，容易溢出，截斷誤差等問題。加之以空間譜估計為研究背景的協方差矩陣運算，大多得到的是針對特殊陣列模型的實對稱矩陣，不具備通用性。針對定點運算的不足和該運算的適用范圍，研究了浮點運算和復數運算的特點，提出了基于復數矢量的浮點協方差運算的FPGA實現方案。在Altera\\\\ stratix \\\\EP1S20F780C7中的仿真和調試結果表明了該方案的有效性。

關鍵詞:協方差矩陣運算; 浮點運算; 復數運算; FPGA

中圖分類號:TN919; TP302文獻標識碼:A

文章編號:1004-373X(2010)15-0121-04

Realization of Complex Floating-point Covariance Matrix Operation Based on FPGA

WEI Yu-mei， LIU Shuai， JIN Ming

(Harbin Institute of Technology， Weihai 264209， China)

Abstract: The covariance matrix operation based on FPGA is realized with the fix-point calculation mode mostly， but there are some problems in the computing process， such as small dynamic range of data processing， easy overflowing， truncation error and so on. Additionally， covariance matrix operation that takes the spatial spectrum estimation as the research background mostly gets the real symmetric matrix which is in connection with special array model， and does not have universality. To direct at the shortness and application scope of the fixed-point operation， and implementing scheme of complex floating-point covariance matrix operation based on FPGA is put forward after the study on the characteristics of the floating point operation and complex operation. The results of simulation and debugging on EP1S20F780C7 indicate that the scheme is effective.

Keywords: covariance matrix operation; floating-point arithmetic; complex arithmetic; FPGA

0 引 言

協方差矩陣的計算是信號處理領域的典型運算，是實現多級嵌套維納濾波器、空間譜估計、相干源個數估計以及仿射不變量模式識別的關鍵部分，廣泛應用于雷達、聲吶、數字圖像處理等領域。采用FPGA(Field Programmable Gate Array)可以提高該類數字信號處理運算的實時性，是算法工程化的重要環節。但是FPGA不適宜對浮點數的處理，對復雜的不規則計算開發起來也比較困難。故目前國內外[1-4]協方差運算的FPGA實現都是采用定點運算方式。

在所有運算都是定點運算的情況下，每次乘法之后數據位寬都要擴大一倍。若相乘后的數據繼續做加減運算，為了保證數據不溢出，還必須將數據位寬擴展一位，而協方差矩陣的運算核心就是乘累加單元，隨著采樣點數的增加，位寬擴展呈線性增加。最終導致FPGA器件資源枯竭，無法實現設計。為了保證算法的實現，必須對中間運算數據進行截斷，將每次累加的結果除2(可以通過移位運算來實現)，以避免溢出。

此外，在應用MUSIC算法[5]時，各種計算都是復數運算。為達到減少算法的計算量，提高MUSIC算法處理速度的目的，許多文獻致力于研究陣列的結構特點，在保證測角精度的前提下，尋找一種簡單而有效的數據預處理方法，將復數矩陣轉化為實數矩陣，把復矢量用一個實矢量來代替，從而將復數運算轉化為實數運算。

接收陣元模型可分為任意離散陣、均勻圓弧陣、均勻圓陣和均勻線陣。在實際應用中，比較常見的是均勻線陣和均勻圓陣。每種陣列模型都有各自的特點，加之陣元數目的取值不同，也會導致陣列流型的對稱性變化。針對不同的陣元模型和陣元數，數據預處理的方法也會有所不同。

對于數據預處理的研究，目前已經有了一些比較成熟的算法。對于一個偶數陣元的對稱陣列(包括均勻線陣和均勻圓陣)，相關研究表明[6-8]，可利用其對稱性，分成兩個完全對稱的子陣，選擇合適的參考點，構造互為共軛對稱的方向矩陣，進而構造一個線性變換矩陣，即可達到將復數矩陣轉化為實數矩陣的目的。

對于奇數陣元的均勻線陣，也有相關研究成果[9]表明，通過構造一個酉矩陣，也可以達到數據預處理的目的。

由于均勻圓陣的陣列流型矩陣不是Vandermonde矩陣，即不具備旋轉不變性，因此適用于奇數陣元的均勻線陣的預處理理論不能直接用于奇數陣元的均勻圓陣，需要將圓陣先轉換到模式空間——虛擬線陣[10]，而轉換需要第一類Bessel函數，不適宜用硬件實現。

以上研究表明，目前除了奇數陣元的均勻圓陣外，其他常用陣列模型都可以通過預處理的方法將復數運算轉換為實數運算。若在某些特定的情況下，必須采用奇數陣元的均勻圓陣。此時，基于復數運算的協方差矩陣的實現就成為一種必然。

因此，在充分應用FPGA并行處理能力的同時，為了擴展數據處理的動態范圍，減少數據溢出機率，避免數據截斷所產生的誤差，提高協方差矩陣的運算精度以及擴展該運算的通用性。本文以空間譜估計作為研究背景，研究了復數據運算和浮點運算的特點，提出了一種適用于任何陣列流型、任意陣元的基于復數浮點運算的協方差矩陣的FPGA實現方案。

1 求解復數浮點協方差矩陣

以11陣元的均勻圓陣為例，其協方差矩陣的求解方案原理框圖如圖1所示。

圖1 復數浮點協方差矩陣實現原理圖

1.1 FIFO數據緩存器

在該設計方案中選擇FIFO作為數據存儲器，這是因為一旦多路接收機有數據輸出，就會啟動FIFO進行存儲，進而FIFO的不空信號有效(empty=0)，觸發后續的矩陣運算;否則，運算停止，一切狀態清零，FPGA恢復idle(空閑)狀態，等待新的快拍采樣數據的到來。這樣可以很方便地控制運算的開始和結束。

矩陣運算所需要的同步時鐘需要設計一個類似于單穩態觸發器的模塊。當檢測到empty=‘0’時，就觸發一個含有121個clk(對于串行方案而言)時鐘信號周期長度的高電平。該高電平與主時鐘相與便可以得到運算的同步時鐘。

1.2 數據共軛轉換

由于測向陣列的輸出矢量X(t)是一個復矢量，對其求協方差矩陣需用陣列輸出列矢量X(t)與其共軛轉置矢量XH(n)對應相乘。如式(1)所示:

Rx(i，j)=1N∑Nk=1[xi(k)xj(k)H]

(1)

1.3 定點數到浮點數的轉換

定點計算在硬件上實現簡單，計算速度比浮點計算要快，但是表示操作數的動態范圍受到限制，浮點數計算硬件實現比較困難，一次計算花費的時間也遠大于定點計算的花費，但是其表示的操作數動態范圍大，精度高。在本設計中，考慮到系統的數據動態范圍和運算精度，選擇浮點計算。由于運算數據是直接從接收機I，Q兩路通道的A/D變換器的輸出獲得，為定點數，因此必須要有一個將A/D采樣的定點數據轉換為浮點數的過程。設計中將16位定點數轉換為IEEE 754標準的單精度格式。32位單精度格式如圖2所示，最高位為符號位，其后8位為指數e(用移碼表示，基數f=2，偏移量為127)，余下的23位為尾數m。

圖2 IEEE 754標準常用的浮點數格式

1.4浮點復數乘累加器

1.4.1 復數乘法器

假設有兩個復數分別為a+jb和c+jd，這兩個數的乘積為:

(a+jb)×(c+jd)=(ac-bd)+j(ad+bc)

復數乘法器的工作原理如圖3所示，其中所用到的加法、減法和乘法器都是基于浮點的運算。值得一提的是，在實現浮點加減法的時候，可以將尾數連同符號位轉化為變形補碼形式后再進行加減運算。這樣做的目的是方便判斷數據是否溢出(變形補碼判斷溢出的規則是:當兩位符號位不同時表示溢出，否則無溢出。無論數據是否溢出，第一位符號位永遠代表真正的符號)，若溢出，則將尾數右歸，指數部分加1，若沒有溢出，則將尾數左歸(規格化)。浮點乘法相對較簡單，對應階碼相加，尾數相乘可以采用定點小數的任何一種乘法運算來完成，只是在限定只取一倍字長時，乘積的若干低位將會丟失，引入誤差。

1.4.2 浮點復數乘累加器

以11個陣元的圓陣為例，實現串行處理方案的浮點復數乘累加器的原理如圖4所示，實部和虛部(雙通道)的乘累加器模塊工作原理一樣。

121階數據緩存器實際上就是121個數據鎖存器級聯形成的一個移位寄存器，初始狀態為零。當浮點復數乘法器有輸出的時候，啟動數據緩存器與之進行加法操作，121個時鐘周期以后可以實現一次快拍采樣的矩陣累加。累加清零信號由時序控制器給出，當所有的快拍采樣點運算都結束之后，數據緩存器輸出累加結果(即協方差矩陣的運算結果)，同時控制器送出一個清零信號，清零121階數據緩存器。

圖3 復數乘法器的原理框圖

圖4 浮點復數乘累加器原理框圖

2 仿真結果

可編程邏輯設計有許多內在規律可循，其中一項就是面積和速度的平衡與互換原則。面積和速度是一對對立統一的矛盾體，要求一個設計同時具備設計面積最小，運行頻率最高，這是不現實的。于是基于面積優先原則和速度優先原則，本文分別設計了協方差矩陣的串行處理方案和并行處理方案，并用Altera\\\\stratix\\\\EP1S20F780C7進行板上調試。其調試結果表明，串行處理方案占用的資源是并行處理方案的1/4，但其運算速度卻是后者的11倍。

2.1 串行處理方案仿真結果

如圖5所示，clk為運算的總控制時鐘;reset為復位控制信號，高電平有效;rd為讀使能信號，低電平有效;wr為寫使能信號，低電平有效;wr_clk為寫時鐘信號，上升沿觸發;q_clk為讀時鐘信號，上升沿觸發;ab_re(31:0)和ab_im(31:0)為乘法器輸出的實部和虛部。q_t2為矩陣乘累加模塊的同步時鐘信號; clk11，state(3:0)，clk1和state(3:0)是狀態機的控制信號，控制矩陣運算規則。

如圖5所示，在100 ns時reset信號有效(即reset=‘1’)，所有狀態清零。從335~635 ns間，寫使能信號有效(wr=‘0’)且有兩個寫時鐘信號的上升沿到來，即向任意一個通道的FIFO中存入兩個快拍采樣數據，最后輸出結果應該有兩個矩陣，如圖6所示。當FIFO為空時，運算停止，所有狀態清零。等待新采樣數據的到來。

圖5 串行方案仿真結果(一)

圖5中，在350 ns時，讀使能有效(rd=‘0’)且有一個讀時鐘信號的上升沿到來，所以empty信號存在短暫的不空(empty=‘0’)狀態，捕獲到這個信息，便觸發單穩態觸發器模塊，產生具有121個clk時鐘周期長度，占空比為120∶1的q_clk信號，進行FIFO的讀操作。

在350~535 ns時間段，因為寫時鐘信號沒有到來，所以FIFO為空(empty=‘1’)。從550 ns~24.75 μs時間段讀時鐘信號沒有上升沿到來，整個設計處于第一個矩陣的運算過程中，即運算一個矩陣所需要的時間為24.2 μs。與此同時，第二個數據寫入FIFO，empty一直處于不空狀態(empty=‘0’)。

在第一個矩陣運算結束之后，即24.6 μs時，系統檢測到empty=‘0’，開始讀數據并觸發第二個矩陣運算的時鐘控制信號。如圖6所示，在24.6 μs時，empty=‘1’。FIFO中的第二個數據被讀出，處于空狀態。從24.85~49.05 μs進入第二個矩陣的運算周期。

圖6 串行方案仿真結果(二)

在仿真時，輸入數據為16位的定點數(1+j1;0+j0;2+j2;3+j3;4+j4;5+j5;6+j6;7+j7;8+j8;9+j9;A+jA)，輸出結果為32位的單精度浮點數。選擇的主時鐘周期為200 ns。在實際調試過程中，整個系統可以在50 MHz主時鐘頻率下正常工作。

2.2 并行處理方案仿真結果

并行方案運算原理與串行方案的一樣，只是在時鐘控制上有所區別，因為采用了11個浮點復數乘累加器，進行一次矩陣運算，只需要11個時鐘周期，如圖7，圖8所示。在仿真時，設置在寫使能信號有效(wr=‘0’)的同時，有3個寫時鐘信號(wr_clk)的上升沿到來，即分別向22個FIFO中存入3個數據，則輸出有3個矩陣。從圖7中還可以清楚地看出，運算結果是矩陣的11行數據并行輸出，輸出結果是一個對稱矩陣。

圖7 并行設計方案仿真結果(一)

圖8 并行設計方案仿真結果(二)

3 結 語

在分析了目前應用于空間譜估計的協方差矩陣運算在硬件實現上的不足，如定點計算的數據動態范圍小，運算精度不高，且只適用于特定陣列模型和的陣元數，不具備通用性。在此基礎上提出了基于浮點運算的通用型協方差矩陣的實現方案。仿真結果表明，本文所提出的實現方案采用的是復數乘法運算，最終結果得到的是復共軛對稱矩陣，適合利用任意的陣列模型和陣元數得到與之相對應的協方差矩陣。這就拓展了協方差矩陣運算的應用范圍，且整個運算過程采用的是浮點運算，提高了整個運算的精度。

參考文獻

[1]KIM Minseok Kim， ICHIGE Koichi. Implementation of FPGA based fast DOA estimator using unitary MUSIC algorithm[J]. IEEE Trans.AP， 2003，1(58): 213-217.

[2]SONG Wen-miao， ZHANG Jing-ying. Research and FPGA implementation of MUSIC DOA of smart antenna array[J]. Intermational Conference on Networks Security， 2009， 2: 189-192.

[3]CHEN Hao， WEI ping. The high-speed implementation of direction-of-arrival estimation algorithm[J]. IEEE International Conference， 2002(6): 922-925.

[4]ZHOU Zou， WANG Hong-yuan， YU Guo-wen. An improved MUSIC algorithm implemented with hign-speed parallel optimization for FPGA[J]. Antennas， Propagation EM Therory， 2006，9: 1-4.

[5]SCHMIDT R O. Multiple emitter location and signal parameter estimation[J]. IEEE Trans.on AP， 1986:3(3): 276-280.

[6]吳仁彪.一種通用的高分辨率波達方向估計預處理新方法[J].電子科學學刊，1993，15(3):305-309.

[7]陳昊.空間譜估計算法的高速實現[D].成都:電子科技大學，2006.

[8]張樂.信號子空間分解的FPGA實現[D].成都:電子科技大學，2006.

[9]何勁，陳建文，劉中.利用空時相關矩陣酉變換估計波達方向[J].系統工程與技術，2005，27(5):811-813.

[10]鐘誠，何培宇.基于空間平滑的均勻圓陣對相干信源的DOA估計[J].成都信息工程學院學報，2007，22(6):723-726.