MUSIC算法的GPU高效實現

朱新峰，陸小霞，黎仁剛，王國海，王彥凌

(1.揚州大學，江蘇 揚州 225001；2.中國船舶重工集團公司第七二三研究所，江蘇 揚州 225101；3.中電科航空電子有限公司，四川 成都 611731;4.浙江理工大學，浙江 杭州 310018)

0 引 言

空間譜搜索[1]是陣列信號處理中的一個重要研究方向，在雷達、通信、聲吶等眾多領域有著極為廣闊的應用前景。目前對于空間譜搜索主要用信號子空間類方法[2]來處理。信號子空間類方法主要由2步組成:一是對接收數據協方差矩陣進行正交分解來估計信號(或噪聲)子空間;二是進行譜峰搜索來提取波達方向估值。無論是正交分解還是譜值計算都面臨著巨大的運算量。在信息技術的不斷更新中，傳統測向技術[3]已經無法滿足當下的需求，因此以MUSIC算法為代表的超分辨測向算法應運而生[4]。然而MUSIC算法中接收矩陣的子空間分解部分計算量龐大，影響了計算效率。為了高效實現MUSIC算法，因此考慮將MUSIC算法放在GPU下去實現，GPU可以提供數十倍乃至于上百倍于CPU的性能，能夠明顯提升運行的速度。

1 MUSIC原理和運算流程

1.1 MUSIC原理

在經典窄帶MUSIC算法中，通常要求入射信號源數目小于陣元數，即P

設快拍數為N，根據式(1)所示接收數據模型，可求得陣列接收數據的協方差矩陣Rxx為：

Rxx= E[x(t)xH(t)]=

AE[s(t)sH(t)]AH+E[n(t)nH(t)]=

ARssAH+σ2I

(1)

式中：E[*]表示求算數平均；Rss表示原始信號的協方差矩陣，當信號間相互獨立時，Rss為P維對角矩陣，對角線上值為對應位置信號的平均功率；σ2為噪聲功率；I為M維單位矩陣。

因為相互獨立的假設前提，不同陣元噪聲間、信號與噪聲間的協方差值接近0，所以接收數據的協方差矩陣可近似為信號矩陣的主對角線上加上噪聲功率。對Rxx進行奇異值分解，得到特征向量矩陣U、V和對應特征值組成的對角矩陣Σ：

Rxx=UΣV

(2)

式中：Σ對角線呈遞減排列時，可根據值大小將特征向量矩陣分為信號子空間Us與噪聲子空間Un。

(3)

U=[UsUn]

(4)

奇異值分解所得的特征向量兩兩正交，所得信號子空間與噪聲子空間必然正交，且信號子空間與入射信號陣列流形張成的空間一致，所以陣列流形中任意導向矢量與噪聲子空間正交。

aiUn=0，i=1,2,…,P

(5)

但在實際情況中，由于噪聲的存在，以及入射信號復包絡近似不變的假設，導致估計結果存在一定偏差，導向矢量與噪聲子空間不可能完全正交。

因此，通過掃描不同入射角構造各來波方向上的陣列流形，求取式(3)模值平方，取得極小值的點即為來波方向，如式(6):

(6)

式中：argmin{*}表示求*最小值時掃描變量的值。

為了更加直觀地表示來波方向，且使接近的2個峰值更易于分辨，繪制空間譜圖時通常取倒數，如式(7):

(7)

式中：Pmusic表示空間譜。

根據MUSIC算法的基本原理，繪出MUSIC算法的一般流程，如圖1所示。

圖1 MUSIC算法流程圖

1.2 算法運行時間資源分配

以五元非均勻陣列為目標陣形，入射方向分別為10°，30°，60°，信噪比為10 dB，高斯白噪聲為500的線性調頻信號作為入射信號，MUSIC算法計算出的空間譜如圖2所示。

圖2 5元陣空間譜

圖2中，3個最大譜峰處的方向角坐標值即為來波方向，算法能夠準確估算出3個入射信號的方向，具有優越的多信號同時測向能力。

在時間上，算法中陣元數為9，使用信噪比是60，中心頻率是18e9，帶寬是20e6，采樣頻率是40e6，根據上面數據計算出采樣時間、快拍持續時間、頻率變化率，并在Matlab平臺下進行仿真，如圖3所示。

圖3 九元陣方向角空間譜

使用的總時間為2.148 926 s，并且單獨測試了計算占用時間、接收矩陣、協方差計算、奇異值分解(SVD)，譜峰搜索分別穩定在0.001 462 s，0.002 166 s，0.003 893 s，1.963 615 s，各占大約0.07%,0.14%,0.21%,0.58%。可以看出其中SVD和譜峰搜索占了大部分的計算時間。本文著重進行SVD和譜峰搜索部分的優化，以提高計算效率。

2 使用GPU提升MUSIC算法運算效率

2.1 MUSIC算法分解

2.1.1 接收數據

寬帶信號與窄帶信號[5]是一個相對概念，不同系統中的寬帶信號與窄帶信號具有不同的區分標準，且寬帶信號與窄帶信號的處理方法也有較大差距。

窄帶信號是指信號帶寬遠遠小于信號波前跨越陣列最大口徑所需要的時間的倒數，即：

(8)

式中：L為陣列最大口徑；F和λ為信號中心頻率和該頻率對應的波長。

2.1.2 奇異值分解

基于單邊Jacobi[6]旋轉的SVD算法,相對于雙邊、單邊的計算量小，并且容易并行實現。單邊Jacobi方法直接對原矩陣A進行單邊正交旋轉，A可以是任意矩陣。AJ1J2…JK=B，同樣每一輪的迭代都要使A的任意兩列都正交一次，迭代退出的條件是B的任意兩列都正交。單邊Jacobi旋轉有這么一個性質：旋轉前若‖ap‖>‖aq‖，則旋轉后依然是‖ap‖>‖aq‖；反之亦然。把矩陣B每列的模長提取出來，B=U·Σ，把J1J2J3…Jk記為V，則A=V·Σ·VT。

2.1.3 譜峰搜索

應用于MUSIC算法的譜峰搜索方法多種多樣，主要有柵格化多細粒度搜索，應用于特殊陣列的免搜索或快速搜索算法，以及最優化求解的智能優化算法。根據系統資源與實現需求的不同，譜峰搜索的選擇也不盡相同[5]。

由于本文系統的實時性要求較高，選取了三級細粒度直接搜索的方法。分別采用5°，1°，0.25°對全方位，一級點[-10°,10°]，二級點[-2°,2°]作為柵格大小，如圖4所示。

圖4 常用譜峰搜索方法

2.2 算法的GPU高效實現

為了高效實現MUSIC算法，將MUSIC算法中部分復雜運算放在GPU下去實現，GPU可以提供數十倍乃至于上百倍CPU的性能，將接收信號的模型、協方差的計算、SVD放入GPU里去并行計算，能夠明顯提升運行的速度。

3 算法仿真

3.1 改進SVD

為了提高復數SVD的計算效率，將獲取列范數，歸并排序，矩陣排序，單邊Jacobi以及奇異值的排序放在GPU里去進行并行計算。CPU與GPU上的運行時間對比如圖5所示。

圖5 CPU與GPU上的運行時間對比

由圖5可以看出，在GPU里的運行時間最終穩定在0.00 255 s，與原先在CPU里的運行時間0.003 893 s相比速度明顯加快。

下面在譜峰搜索部分GPU的實現中進行了使用共享內存、分配鎖頁內存和增加CUDA流操作3個方面的優化。

3.2 使用共享內存

首先將信噪比、中心頻率、帶寬、采樣頻率、采樣時間、快拍持續時間、頻率變化率設置為全局變量，減少大量的內存獲取請求，從而減少時間的交互。

然后再將譜峰搜索中的計數器設置為共享內存，在調用過程中便不需要進行繁瑣的交換數據。并且將調用__syncthread()函數來完成線程同步，這個函數的調用將確保線程塊中的每個線程都執行完__syncthreads()前面的語句，才會執行下一條語句。這樣可以確保對共享內存進行讀取之前，想要寫入的操作已經完成。優化前后的時間對比如圖6所示。

圖6 GPU上使用共享內存前后的時間對比

從圖6可以看出，沒有使用共享內存時，GPU的時間穩定在350 ms左右，使用共享內存后時間穩定在315 ms左右，有了執行時間降低了35 ms左右，性能提升10%左右。

3.3 分配鎖頁內存

PCI-E數據總線傳輸的數據始終是鎖頁內存中的數據，若應用程序沒有對鎖頁內存進行分配，則CUDA驅動會在后臺為程序分配。這之間必然存在分頁內存與鎖頁內存中不必要的拷貝操作，所以，通過分配鎖頁內存來減少這些拷貝操作，提高效率[8]。

在應用程序中，將主機端分配內存的代碼塊:

double*host_n=(double *)malloc(size_of_n_in_bytes);

替換為：

double *host;

CUDA_CALL(cudaMallocHost(&host_n,size_of_n_in_bytes));

并將最后釋放內存的代碼：

Free(host_n);

替換為：

CUDA_CALL(host_n)。

測試運行時間，仿真結果如圖7所示。

圖7 GPU上分頁內存與鎖存內存的時間對比

從圖7可以看出，使用分頁內存時運行時間穩定在315 ms，使用鎖頁內存時穩定在310 ms左右，時間減少5 ms，性能上提升了1.6%，提升的空間并沒有前面使用的共享內存優化那么明顯。

3.4 增加CUDA流操作

CUDA流[8]表示一個GPU操作隊列，并且該隊列中的操作將以指定的順序執行。我們可以在流中添加一些操作，如核函數啟動、內存復制等。將這些操作添加到流的順序也就是他們的執行順序。將每個流視為GPU上的一個任務，并且這些任務可以并行執行。

首先是創建一個流：

cudaStream_t aes_async_stream;

CUDA_CALL(cudaStreamCreate(&aes_async_stream));

相對應地會有銷毀流：

CUDA_CALL(cudaSteamDestroy(aes_async_stream));

然后將需要的事件壓入流中：

CUDA_CALL(cudaEventRecord(start_time,aes_async_stream));

CUDA_CALL(cudaMemcpyAsync(dev_xt,host_xt,size_of_xt_in_bytes,cudaMemcpyHostToDevice,aes_async_stream));

運行時間如圖8所示。

圖8 GPU上使用CUDA流前后的時間對比

從圖8可以看出，在沒有使用CUDA流之前，時間穩定在310 ms左右，使用CUDA流之后，時間穩定在284 ms左右，時間減少26 ms，性能上提升9.15%。

4 結束語

本文針對MUSIC算法的高效GPU實現，經過復數SVD的GPU并行計算，在譜峰搜索部分GPU實現中使用共享內存、分配鎖頁內存、增加CUDA流操作等幾個方面的改進，整個計算過程減少56 ms的時間。總的來說，性能提升了20%的效率。實驗結果表明，性能的提升能夠符合預想的結果。