基于迭代優化算法的MVDR空間譜估計的SoC實現

廖　斌，李思敏，2,唐智靈,2

(1.桂林電子科技大學 信息與通信學院,廣西 桂林　541004；2.桂林電子科技大學 認知無線電與信息處理教育部重點實驗室，廣西 桂林　541004)

基于迭代優化算法的MVDR空間譜估計的SoC實現

廖斌1，李思敏1，2,唐智靈1,2

(1.桂林電子科技大學 信息與通信學院,廣西 桂林541004；2.桂林電子科技大學 認知無線電與信息處理教育部重點實驗室，廣西 桂林541004)

摘要：針對MVDR算法中矩陣求逆運算量大、耗時長，提出一種基于迭代優化的坐標下降法。該算法將逆矩陣求逆轉化為線性方程組的最優解，以迭代的方式代替傳統的矩陣運算。同時采用FPGA的SoC設計方式，在處理器上搜索譜峰，并通過專門的硬件實現迭代優化算法，以加速求逆運算。仿真結果和硬件測試表明，計算結果與實際結果的誤差很小，基于SoC的設計能有效地搜索出譜峰。

關鍵詞：陣列信號處理；MVDR；DOA；逆矩陣；SoC

來波方向(direction of arrival,簡稱DOA)估計主要用于測量入射信號的來波方位。在智能天線應用中，一旦有了信號的方位信息，就可以根據需求判斷是否將最大波束對準該方向。而在電子偵查中，則可以通過來波方位信息對輻射源定位，以確定輻射源的相對位置[1]。

典型的測向算法有MUSIC(multiple signal classification)算法、ESPRIT(estimate signal parameter via rotational invariance techniques)算法以及MVDR(minimum variance distortionless response)算法。其中，MUSIC算法利用信號子空間和噪聲子空間的正交性，構建空間譜函數[2]；ESPRIT算法利用自相關矩陣信號子空間的旋轉不變特性構建空間譜函數[3]；而MVDR算法是在MVDR波束成形基礎上得到的空間譜函數，該波束成形是在保持期望信號增益不變的情況下，使其他觀測方向的輸出平均功率最小[4]。

通常基于MVDR的DOA估計需要完成協方差的求逆運算，而典型的求逆方式是采用一種基于systolic結構的QR分解算法[5-7]，通過并行運算減少處理時間。但是，由于占用了大量的矩陣運算，使其消耗更多的硬件資源。為此，提出一種基于迭代優化的坐標下降法，通過求解線性方程的迭代最優解代替矩陣的分解運算，減少硬件資源，縮短運算時間。同時，為了便于完成譜峰搜索，采用FPGA的SoC(system on chip)設計方案，利用基于迭代優化算法的硬件IP核完成求逆加速運算，并通過軟件的方式進行譜峰搜索算法的處理，簡化了硬件設計。

1相關工作

1.1陣列信號模型

假設由N個陣元組成均勻線陣，線陣間距為d，若有D個窄帶信號sk(t),k=1,2,…,D入射,其對應來波方向θk=(θ1,θ2,…,θD)，則第i個陣元接收到的信號:

(1)

其中：ni(t)為第i個陣元的加性噪聲；λ為信號波長。

整個陣列接收的信號表示為:

(2)

陣列的加權輸出：

(3)

其中w=[w1,w2,…,wN]T。

1.2基于MVDR的DOA估計

MVDR是一種最小方差法，其在信源方向θk上信號功率不變，而非信源方向的任何干擾所貢獻的功率為最小，即

(4)

根據Lagrange乘數法，可得

(5)

其中Rxx為陣列的協方差矩陣，有

(6)

實際處理中，無法得到理想的協方差矩陣，通常采用最大似然估計：

(7)

其中：x(n)=[x1(n),x2(n),…,xN(n)]T為第n時刻陣元的輸出，稱為一次快拍；K為快拍數。

根據式(4)，可得空間功率譜函數：

(8)

DOA估計就是對所有θ進行式(8)的計算，得到每個角度對應的功率，并搜索其譜峰，譜峰對應的角度就是其波達角。

1.3坐標下降法

根據式(7)，設x=R-1a(θ)，則有Rx=a(θ)，可得

(9)

因此，可將MVDR的求逆變換為線性方程組Rx=a(θ)的求解。

求解線性方程組有直接法和迭代法，其中，坐標下降法是求解線性方程組的一類迭代算法，其基于準確線性搜尋方式來求解線性方程組Ax=b，A∈RNN的通用算法步驟如下[8]：

1)r(0)←b-Ax(0)，p(k)←f(k)，一般取x(0)=0。

2)Fork=1,2,…

其中：a(k)為搜索步徑；p(k)為搜索方向矢量。

令p(k)=p(k+N)，且滿足p(n-1)=en，n=1,2,…,N，en∈RN，其中en表示第n個元素為1、其余元素為0的單位列向量，可得坐標下降法，其有如下關系：

(10)

(11)

(12)

1)r(0)=b-Ax(0)，p(0)=e1，k=0，x(0)=0。

2)Forn=1∶N

k=k+1;

if( k>kmax)stop

End

2坐標下降法的迭代優化及仿真

2.1基于迭代優化的坐標下降法

(13)

(14)

故可得到基于迭代優化的坐標下降法的偽代碼形式：

1)初始化變量：a(k)=2p，r=b，x=0，m=1，k=1，M、H為常數。

2)For m=1∶M

a(k)=a(k)/4;

For k=1∶H

For n=1∶N

xn=xn+sign(rn)a(k);

r=r-sign(rn)A:,n;

End

End

End

偽代碼中的常數M為對搜索步徑進行倍數縮放的次數，常數H為按確定步徑進行搜索的次數。改進后的迭代優化算法，由于a(k)為2的倍數次冪，可用移位方式代替乘法器，故算法共用了N×(N+1)×H×M個加法，而之前的坐標下降法共用了(N+1)×k個加法、N×k個乘法、k個除法。

(15)

其中：

2.2迭代優化算法的仿真

仿真條件：由四陣元構成均勻線陣，載波頻率為15 MHz，間距d=10 m，快拍數K=256，有2個信號源，入射角度分別為-10°、25°。改進的坐標下降法的初始步徑設為a(k)=2。

圖1、圖2為不同迭代參數下迭代優化算法和基于Matlab直接求逆函數的空間譜。由圖2可知，當參數M=8，H=8，基于迭代優化的坐標下降法的計算曲線幾乎與理論的數值曲線重合，即迭代優化算法的計算結果與理論計算值誤差很小。比較圖1和圖2可得，迭代算法的參數M、H影響計算結果的準確性，M和H的值越大，計算結果越接近理論值。

圖1　MVDR空間譜(H=4)Fig.1　MVDR spatial spectrum(H=4)

圖2　MVDR空間譜(H=8)Fig.2　MVDR spatial spectrum(H=8)

3DOA估計的SoC設計

3.1總體方案

圖3為基于FPGA的DOA的SoC構架，通過AXI總線實現軟核處理器Microblaze與片上外設的數據通信，其中，在FPGA的邏輯部分完成基于最大似然估計的協方差模塊的IP核、基于優化迭代算法的IP核的設計，而在軟核處理器上譜峰搜索，通過求極小值得到其譜峰，并在每搜尋一個角度時，調用迭代優化算法的IP核完成求逆的硬件加速，而其他外設如UART(串口)，則負責程序的調試，定時器計算程序運行時間。

圖3　基于FPGA的DOA SoC構架Fig.3　SoC architecture of DOA based on FPGA

3.2協方差模塊的IP核設計

根據式(7)，協方差矩陣可寫成：

基于AXI總線協方差模塊的IP核設計如圖4所示，快拍數設K=2q,q∈N，其中AXI Bram Control主要完成AXI總線地址與Dual RAM(雙端口RAM)地址的轉換。IP核的軟硬件協同工作如下：

1)處理器通過AXI總線讀取IP核的狀態寄存器中的busy位，判斷其是否為1，若是，則等待，直到busy=0，跳轉至下一步。

2)往寄存器K寫入快拍數，并使能控制寄存器的start位，使能有效后，雙端口RAM的R1數值將清零，AD采樣模塊開始工作，接收當前的快拍，并從雙端口RAM x的A口寫入，在下個周期從B口讀出送入MUY Writer模塊，經過相乘得到串行x(n)x(n)H，并與從雙端口RAM R1的B口讀出的數值R1(n-1)相加，得到新值。之后從A口回寫,覆蓋雙端口RAM的R1數值，雙端口RAM的R1主要功能是存儲累加的數值，直至K-1次累加后，送入桶形移位寄存器進行q位右移，移位后所得數據從A口寫入存放復數協方差的雙端口complex RAM中，最后根據式(15)，通過一個RAM Control模塊，將復數協方差RAM的數值轉換成實數，并寫入實數協方差雙端口real RAM中。

3)直到狀態寄存器的finish=1，通過總線讀取實數協方差雙端口RAM的數值。

圖4　協方差模塊的IP核設計Fig.4　Design of covariance module IP core

3.3基于迭代優化算法的IP核設計

圖5為迭代優化算法的實現框圖,以求解Ax=b,A∈RNN，其中，圖5的雙端口RAM r的初值為b，狀態機負責模塊的協調工作。整個IP核的軟硬件協同工作如下：

1)處理器通過AXI總線讀取IP核的狀態寄存器中的busy位，判斷其是否為1，若為1，則一直等待，直到busy=0，跳轉至下一步；

2)通過總線在雙端口RAM A寫入協方差矩陣、雙端口RAM r寫入方向矢量，并將雙端口RAM x清零；

3)寄存器M、H寫入迭代數值，寄存器a(k)寫入初始步徑，并使能控制寄存器的start位，使狀態機開始工作，直到狀態寄存器的finish=1，通過AXI總線將雙端口RAM x的值讀出。

狀態機的工作流程為：

1)狀態機一開始處于idle狀態，直到start信號為高，狀態機由idle狀態轉為state0，計數器m置0；

2)當狀態處于state0，判斷計數器m是否小于寄存器M，若是，則狀態跳轉到state1，計數器k置0，寄存器a(k)右移2位，否則，狀態跳轉到idle，狀態寄存器finish置1；

3)當狀態處于state1，判斷計數器k是否小于寄存器H，若是，則狀態跳轉到state2，計數器n置0，否則，狀態跳轉到state0狀態，計數器m加1；

4)當狀態處于state2，判斷計數器n是否小于N，若是，則跳轉到state3，否則，狀態跳轉到state1，計數器k加1；

5)當狀態處于state3，判斷done信號是否為高，若是，則狀態跳轉到state2，否則狀態一直維持在state3。

state3的工作主要是完成xn=xn+sign(rn)a(k)與r=r-sign(rn)a(k)A:,n的計算。如圖5所示，分別從雙端口RAM A和雙端口RAM r的B口讀取當前的A:,n與r，放入受控于當前a(k)的桶形移位寄存器得到a(k)A:,n，將a(k)A:,n和取反后的-a(k)A:,n送入二選一多路復用器，rn的最高位連接到二選一多路復用器的選擇端，當rn為正(最高位為0)，多路復用選擇的輸出為-a(k)A:,n，當rn為負(最高位為1)，選擇的輸出為a(k)A:,n，故多路復用器的輸出可表達為-sign(rn)a(k)A:,n，再將多路復用器的輸出與r相加得到新的r值，再從雙端口RAM r的B口回寫，覆蓋原來的值。

同理，rn的最高位連接到第2個多路復用器的選擇端，當rn為正，多路復用選擇的輸出為a(k)，當rn為負數，選擇的輸出為-a(k)，故多路復用器的輸出可表達為sign(rn)a(k)，最后將多路復用的輸出與從雙端口RAM x的B口讀取的xn相加得到新值，再從雙端口RAM的B口回寫，覆蓋原來的值。

圖5　基于迭代優化算法的IP核設計Fig.5　Design of IP core based on iterative optimization algorithm

從圖5可看出，影響IP核最高時鐘頻率的關鍵路徑是2個多路復用器的輸出與加法器、桶形移位寄存器輸出與第一個多路復用器的輸入，因此，需在上述關鍵路徑中插入觸發器D，以提高最高運行時鐘頻率。

基于迭代優化算法的IP核在Xilinx的Kintex 7 xc7k160t芯片上完成了驗證。其中，雙端口RAM A的元素用32 bit的有符號整數表示，雙端口RAM r和x的元素均用32 bit定點有符號數表示，小數位占20 bit。根據綜合布線結果，表1為迭代優化算法與其他算法的所耗硬件資源對比，其中，迭代優化算法的迭代參數為a(k)=2，M=8，H=8。從表1可看出，基于改進的迭代算法雖然所需處理時間并不一定是最優的，但是它能以更少的資源獲得更小的處理延時。同時，需要注意的是，由于迭代優化算法在實現的過程中采用的是一種串行流水線的結構，故無論矩陣的維數是多少，其所消耗的主要硬件資源差別不大，不同之處是處理所消耗的時間。

表1　優化迭代算法與其他算法的硬件資源對比

3.4DOA估計的軟硬件協同設計

波達角的估計可認為是極小值的搜尋，即在整個空間譜上找出所有的極大值[7]，再找出大于某個給定數值的極小值，則其對應的角度就是波達角。根據式(8)，MVDR空間譜表達式的分子為1，因此，可將DOA估計轉化為尋找a(θ)HR-1a(θ)在所有角度θ上的極小值，并滿足小于某個門限值Q。在軟件處理時，為了減少處理器的運算時間，預先定義一組8位有符號的常數組存放每個角度對應的方向矢量，該方向矢量是由a(θ)歸一化后的數值進行127倍放大，并經過式(15)轉換的實數方向矢量。DOA估計軟硬件協同的程序偽代碼如下：

1)定義數組c，存放空間譜的功率值。定義常數組存放方向矢量a(θ),θ=[-90°,90°]。

2)i=0,θ=-90°。調用協方差IP核求出協方差R。

4)i=1,θ=-89°,umin=q[0]。

while (θ≤90°)

{ 調用迭代算法IP核，計算當前x=R-1a(θ)；

if (c[i-1]≥c[i])

{umin=c[i]; }

i=i+1;θ=θ+1;

}

5) Q=2×umin,i=1;

while (i≤180°)

{ if (c[i]≤Q)

{ if (c[i-1]≥c[i])

{ y=c[i];z=i-90°;j=0; }

else

{j=j+1;

if (j= =1) 串口輸出變量y、z的

值；}

}

i=i+1;

}

偽代碼中的第4)步主要尋找空間譜的最小值，角度的步徑為1°，并且每對一個角度求逆時，均要調用迭代優化算法的IP核，以減少運算時間。第5)步主要完成空間譜極小值的搜尋，其門限值Q設為最小值的2倍，并從串口輸出變量z，該變量就是其波達角。

3.5硬件實驗

測試條件：四陣元的均勻線陣，線距d=10 m，接收機接收15.702 MHz的調幅信號，而在遠端5 km處發射一個同頻干擾源，發射功率為3 W，方位角約為5°。

通過陣列接收機將四陣元的信號下變頻至10 kHz處，并經過AD采樣，在FPGA進行基帶處理。根據改進坐標下降法的不同迭代參數得到表2，表2中M=8，H=8所對應MVDR的實測空間譜，如圖6所示。

表2　不同迭代參數的數值統計

圖6　基于MVDR的空間譜Fig.6　Spatial spectrum based on MVDR

圖6和表2所用數據是由片上系統的UART模塊輸出，而表2中耗時統計則由片上系統的定時器進行估算。由表2和圖6可知，2個信號的波達角為5°和43°，并且完成DOA估計所用時間約40 ms。同時，比較表2不同參數下所耗時間可知，由于求逆的迭代運算是用專門的硬件完成的，求逆矩陣需要的時間較短，故隨著迭代次數的增加，其所耗時間相對處理器不會明顯增加。

4結束語

為了完成DOA估計，簡化FPGA的硬件設計，提出了一種基于FPGA的SoC的設計方案。通過軟硬件協同設計，完成譜峰搜索算法實現。同時為了兼顧實時性，對坐標下降法進行改進，將求逆運算中的大量矩陣變換運算轉換為線性方程組的最優求解，并用移位的方式代替乘法運算。數值仿真結果表明了該算法計算結果的正確性，并在FPGA上完成了迭代優化算法的IP核實現。根據綜合布線后的結果，其能以較少的資源獲得較短的運算時間。為了驗證方案實際效果，搭建了整機測試平臺。實驗結果表明，基于SoC的設計方案能在40 ms內正確地測出信號的波達方向。在未來的波束成形芯片集成化中，該方案也具有一定參考價值。

參考文獻:

[1]楊小牛．軟件無線電技術與應用[M].北京：北京理工大學出版社，2010:459-461.

[2]WANG Qihui,WANG Lidong,KANG An,et al.DOA estimation of smart antenna signal based on MUSIC algorithm[J].Journal of Networks,2014,9(5):1309-1316.

[3]ROY R,KAILATH T.ESPRIT-Estimation of signal parameters via rotational invariance techniques[J].IEEE Transactions on Acoustics Speech and Signal Processing,1989,37(7):984-995.

[4]CAPON J.High-resolution frequency-wavenumber spectrum analysis[J].Proceedings of the IEEE,1967,57(8):1408-1418.

[5]ZHANG Jianfeng,CHOW P,LIU Hengzhu.An efficient FPGA implementation of QR decomposition using a novel systolic array architecture based on enhanced vectoring CORDIC[C]//International Conference on Field-Programmable Technology,2014:123-130.

[6]朱勇旭，吳斌，周玉梅，等.用于MIMO-OFDM系統QR分解的分布式脈動陣列處理算法[J].電子與信息學報，2012,34(8):1963-1973.

[7]KARKOOTI M,CAVALLARO J,DICK C.FPGA implementation of matrix inversion using QRD-RLS algorithm[C]//Asilomar Conference on Signals,Systems,and Computers,2005:1625-1629.

[8]WATKINS D S.Fundamentals of Matrix Computations[M].New York:John Wiley Sons,2002:559-570.

[9]LIU Jie.DCD algorithm: architectures,FPGA implementations and applications[D].York:University of York,2008:37-45.

[10]MA L,DICKSON K,MCALLISTER J,et al.QR Decomposition-based Matrix inversion for high performance embedded MIMO receivers[J].IEEE Transcations on Signal Processing,2011,59(4):1858-1867.

[11]YU Hanguang.FPGA-based implementation of QR decomposition[D].Arizona:Arizona State University,2014:50-80.

[12]魏嬋娟，張春水，劉健.一種基于Cholesky分解的快速矩陣求逆方法設計[J].電子設計工程，2014,22(1):159-164.

編輯：梁王歡

SoC implementation of MVDR spatial spectrum estimation based on iterative optimization algorithm

LIAO Bin1, LI Simin1,2, TANG Zhiling1,2

(1.School of Information and Communication Engineering, Guilin University of Electronic Technology, Guilin 541004, China;2.Key Laboratory of Cognitive Radio and Information Processing, Guilin University of Electronic Technology, Guilin 541004, China)

Abstract：Aiming at the large amount of calculation and time-consuming in matrix inversion of MVDR algorithm, a coordinate descent method based on iterative optimization is proposed. The matrix inversion is converted into optimum solution of a linear equation. Meanwhile, SoC design based on FPGA is adopted, the spectrum peaks are searched in the processor and the iterative optimization algorithm is implemented on the dedicated hardware for accelerating matrix inversion. The simulation result shows that the error between the calculated result and the actual result is very small. The hardware test shows that the SoC design can effectively search the peaks of the spectrum.

Key words：array signal processing; MVDR; DOA; matrix inversion; SoC

收稿日期：2016-01-16

基金項目：國家自然科學基金(61461013)

通信作者：李思敏(1963-)，男，江蘇蘇州人，教授，博士，研究方向為電磁場與電磁波、認知無線電。E-mail:siminl@guet.edu.cn

中圖分類號：TN911.72

文獻標志碼：A

文章編號：1673-808X(2016)02-0087-07

引文格式： 廖斌，李思敏，唐智靈.基于迭代優化算法的MVDR空間譜估計的SoC實現[J].桂林電子科技大學學報，2016,36(2)：87-93.