摘 要:基于FPGA的協(xié)方差矩陣運(yùn)算的實(shí)現(xiàn)大多采用的是定點(diǎn)計(jì)算方式，在運(yùn)算過程中，存在數(shù)據(jù)處理動(dòng)態(tài)范圍小，容易溢出，截?cái)嗾`差等問題。加之以空間譜估計(jì)為研究背景的協(xié)方差矩陣運(yùn)算，大多得到的是針對(duì)特殊陣列模型的實(shí)對(duì)稱矩陣，不具備通用性。針對(duì)定點(diǎn)運(yùn)算的不足和該運(yùn)算的適用范圍，研究了浮點(diǎn)運(yùn)算和復(fù)數(shù)運(yùn)算的特點(diǎn)，提出了基于復(fù)數(shù)矢量的浮點(diǎn)協(xié)方差運(yùn)算的FPGA實(shí)現(xiàn)方案。在Altera\\\\ stratix \\\\EP1S20F780C7中的仿真和調(diào)試結(jié)果表明了該方案的有效性。

關(guān)鍵詞:協(xié)方差矩陣運(yùn)算; 浮點(diǎn)運(yùn)算; 復(fù)數(shù)運(yùn)算; FPGA

中圖分類號(hào):TN919; TP302文獻(xiàn)標(biāo)識(shí)碼:A

文章編號(hào):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 引 言

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

在所有運(yùn)算都是定點(diǎn)運(yùn)算的情況下，每次乘法之后數(shù)據(jù)位寬都要擴(kuò)大一倍。若相乘后的數(shù)據(jù)繼續(xù)做加減運(yùn)算，為了保證數(shù)據(jù)不溢出，還必須將數(shù)據(jù)位寬擴(kuò)展一位，而協(xié)方差矩陣的運(yùn)算核心就是乘累加單元，隨著采樣點(diǎn)數(shù)的增加，位寬擴(kuò)展呈線性增加。最終導(dǎo)致FPGA器件資源枯竭，無(wú)法實(shí)現(xiàn)設(shè)計(jì)。為了保證算法的實(shí)現(xiàn)，必須對(duì)中間運(yùn)算數(shù)據(jù)進(jìn)行截?cái)啵瑢⒚看卫奂拥慕Y(jié)果除2(可以通過移位運(yùn)算來(lái)實(shí)現(xiàn))，以避免溢出。

此外，在應(yīng)用MUSIC算法[5]時(shí)，各種計(jì)算都是復(fù)數(shù)運(yùn)算。為達(dá)到減少算法的計(jì)算量，提高M(jìn)USIC算法處理速度的目的，許多文獻(xiàn)致力于研究陣列的結(jié)構(gòu)特點(diǎn)，在保證測(cè)角精度的前提下，尋找一種簡(jiǎn)單而有效的數(shù)據(jù)預(yù)處理方法，將復(fù)數(shù)矩陣轉(zhuǎn)化為實(shí)數(shù)矩陣，把復(fù)矢量用一個(gè)實(shí)矢量來(lái)代替，從而將復(fù)數(shù)運(yùn)算轉(zhuǎn)化為實(shí)數(shù)運(yùn)算。

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

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

對(duì)于奇數(shù)陣元的均勻線陣，也有相關(guān)研究成果[9]表明，通過構(gòu)造一個(gè)酉矩陣，也可以達(dá)到數(shù)據(jù)預(yù)處理的目的。

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

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

因此，在充分應(yīng)用FPGA并行處理能力的同時(shí)，為了擴(kuò)展數(shù)據(jù)處理的動(dòng)態(tài)范圍，減少數(shù)據(jù)溢出機(jī)率，避免數(shù)據(jù)截?cái)嗨a(chǎn)生的誤差，提高協(xié)方差矩陣的運(yùn)算精度以及擴(kuò)展該運(yùn)算的通用性。本文以空間譜估計(jì)作為研究背景，研究了復(fù)數(shù)據(jù)運(yùn)算和浮點(diǎn)運(yùn)算的特點(diǎn)，提出了一種適用于任何陣列流型、任意陣元的基于復(fù)數(shù)浮點(diǎn)運(yùn)算的協(xié)方差矩陣的FPGA實(shí)現(xiàn)方案。

1 求解復(fù)數(shù)浮點(diǎn)協(xié)方差矩陣

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

圖1 復(fù)數(shù)浮點(diǎn)協(xié)方差矩陣實(shí)現(xiàn)原理圖

1.1 FIFO數(shù)據(jù)緩存器

在該設(shè)計(jì)方案中選擇FIFO作為數(shù)據(jù)存儲(chǔ)器，這是因?yàn)橐坏┒嗦方邮諜C(jī)有數(shù)據(jù)輸出，就會(huì)啟動(dòng)FIFO進(jìn)行存儲(chǔ)，進(jìn)而FIFO的不空信號(hào)有效(empty=0)，觸發(fā)后續(xù)的矩陣運(yùn)算;否則，運(yùn)算停止，一切狀態(tài)清零，F(xiàn)PGA恢復(fù)idle(空閑)狀態(tài)，等待新的快拍采樣數(shù)據(jù)的到來(lái)。這樣可以很方便地控制運(yùn)算的開始和結(jié)束。

矩陣運(yùn)算所需要的同步時(shí)鐘需要設(shè)計(jì)一個(gè)類似于單穩(wěn)態(tài)觸發(fā)器的模塊。當(dāng)檢測(cè)到empty=‘0’時(shí)，就觸發(fā)一個(gè)含有121個(gè)clk(對(duì)于串行方案而言)時(shí)鐘信號(hào)周期長(zhǎng)度的高電平。該高電平與主時(shí)鐘相與便可以得到運(yùn)算的同步時(shí)鐘。

1.2 數(shù)據(jù)共軛轉(zhuǎn)換

由于測(cè)向陣列的輸出矢量X(t)是一個(gè)復(fù)矢量，對(duì)其求協(xié)方差矩陣需用陣列輸出列矢量X(t)與其共軛轉(zhuǎn)置矢量XH(n)對(duì)應(yīng)相乘。如式(1)所示:

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

(1)

1.3 定點(diǎn)數(shù)到浮點(diǎn)數(shù)的轉(zhuǎn)換

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

圖2 IEEE 754標(biāo)準(zhǔn)常用的浮點(diǎn)數(shù)格式

1.4浮點(diǎn)復(fù)數(shù)乘累加器

1.4.1 復(fù)數(shù)乘法器

假設(shè)有兩個(gè)復(fù)數(shù)分別為a+jb和c+jd，這兩個(gè)數(shù)的乘積為:

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

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

1.4.2 浮點(diǎn)復(fù)數(shù)乘累加器

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

121階數(shù)據(jù)緩存器實(shí)際上就是121個(gè)數(shù)據(jù)鎖存器級(jí)聯(lián)形成的一個(gè)移位寄存器，初始狀態(tài)為零。當(dāng)浮點(diǎn)復(fù)數(shù)乘法器有輸出的時(shí)候，啟動(dòng)數(shù)據(jù)緩存器與之進(jìn)行加法操作，121個(gè)時(shí)鐘周期以后可以實(shí)現(xiàn)一次快拍采樣的矩陣?yán)奂?。累加清零信?hào)由時(shí)序控制器給出，當(dāng)所有的快拍采樣點(diǎn)運(yùn)算都結(jié)束之后，數(shù)據(jù)緩存器輸出累加結(jié)果(即協(xié)方差矩陣的運(yùn)算結(jié)果)，同時(shí)控制器送出一個(gè)清零信號(hào)，清零121階數(shù)據(jù)緩存器。

圖3 復(fù)數(shù)乘法器的原理框圖

圖4 浮點(diǎn)復(fù)數(shù)乘累加器原理框圖

2 仿真結(jié)果

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

2.1 串行處理方案仿真結(jié)果

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

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

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

圖5中，在350 ns時(shí)，讀使能有效(rd=‘0’)且有一個(gè)讀時(shí)鐘信號(hào)的上升沿到來(lái)，所以empty信號(hào)存在短暫的不空(empty=‘0’)狀態(tài)，捕獲到這個(gè)信息，便觸發(fā)單穩(wěn)態(tài)觸發(fā)器模塊，產(chǎn)生具有121個(gè)clk時(shí)鐘周期長(zhǎng)度，占空比為120∶1的q_clk信號(hào)，進(jìn)行FIFO的讀操作。

在350~535 ns時(shí)間段，因?yàn)閷憰r(shí)鐘信號(hào)沒有到來(lái)，所以FIFO為空(empty=‘1’)。從550 ns~24.75 μs時(shí)間段讀時(shí)鐘信號(hào)沒有上升沿到來(lái)，整個(gè)設(shè)計(jì)處于第一個(gè)矩陣的運(yùn)算過程中，即運(yùn)算一個(gè)矩陣所需要的時(shí)間為24.2 μs。與此同時(shí)，第二個(gè)數(shù)據(jù)寫入FIFO，empty一直處于不空狀態(tài)(empty=‘0’)。

在第一個(gè)矩陣運(yùn)算結(jié)束之后，即24.6 μs時(shí)，系統(tǒng)檢測(cè)到empty=‘0’，開始讀數(shù)據(jù)并觸發(fā)第二個(gè)矩陣運(yùn)算的時(shí)鐘控制信號(hào)。如圖6所示，在24.6 μs時(shí)，empty=‘1’。FIFO中的第二個(gè)數(shù)據(jù)被讀出，處于空狀態(tài)。從24.85~49.05 μs進(jìn)入第二個(gè)矩陣的運(yùn)算周期。

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

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

2.2 并行處理方案仿真結(jié)果

并行方案運(yùn)算原理與串行方案的一樣，只是在時(shí)鐘控制上有所區(qū)別，因?yàn)椴捎昧?1個(gè)浮點(diǎn)復(fù)數(shù)乘累加器，進(jìn)行一次矩陣運(yùn)算，只需要11個(gè)時(shí)鐘周期，如圖7，圖8所示。在仿真時(shí)，設(shè)置在寫使能信號(hào)有效(wr=‘0’)的同時(shí)，有3個(gè)寫時(shí)鐘信號(hào)(wr_clk)的上升沿到來(lái)，即分別向22個(gè)FIFO中存入3個(gè)數(shù)據(jù)，則輸出有3個(gè)矩陣。從圖7中還可以清楚地看出，運(yùn)算結(jié)果是矩陣的11行數(shù)據(jù)并行輸出，輸出結(jié)果是一個(gè)對(duì)稱矩陣。

圖7 并行設(shè)計(jì)方案仿真結(jié)果(一)

圖8 并行設(shè)計(jì)方案仿真結(jié)果(二)

3 結(jié) 語(yǔ)

在分析了目前應(yīng)用于空間譜估計(jì)的協(xié)方差矩陣運(yùn)算在硬件實(shí)現(xiàn)上的不足，如定點(diǎn)計(jì)算的數(shù)據(jù)動(dòng)態(tài)范圍小，運(yùn)算精度不高，且只適用于特定陣列模型和的陣元數(shù)，不具備通用性。在此基礎(chǔ)上提出了基于浮點(diǎn)運(yùn)算的通用型協(xié)方差矩陣的實(shí)現(xiàn)方案。仿真結(jié)果表明，本文所提出的實(shí)現(xiàn)方案采用的是復(fù)數(shù)乘法運(yùn)算，最終結(jié)果得到的是復(fù)共軛對(duì)稱矩陣，適合利用任意的陣列模型和陣元數(shù)得到與之相對(duì)應(yīng)的協(xié)方差矩陣。這就拓展了協(xié)方差矩陣運(yùn)算的應(yīng)用范圍，且整個(gè)運(yùn)算過程采用的是浮點(diǎn)運(yùn)算，提高了整個(gè)運(yùn)算的精度。

參考文獻(xiàn)

[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]吳仁彪.一種通用的高分辨率波達(dá)方向估計(jì)預(yù)處理新方法[J].電子科學(xué)學(xué)刊，1993，15(3):305-309.

[7]陳昊.空間譜估計(jì)算法的高速實(shí)現(xiàn)[D].成都:電子科技大學(xué)，2006.

[8]張樂.信號(hào)子空間分解的FPGA實(shí)現(xiàn)[D].成都:電子科技大學(xué)，2006.

[9]何勁，陳建文，劉中.利用空時(shí)相關(guān)矩陣酉變換估計(jì)波達(dá)方向[J].系統(tǒng)工程與技術(shù)，2005，27(5):811-813.

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