曾富紅, 司偉建, 曲志昱

(哈爾濱工程大學 信息與通信工程學院, 黑龍江 哈爾濱 150001)

?

基于Hermitian矩陣的特征分解算法

曾富紅, 司偉建, 曲志昱

(哈爾濱工程大學 信息與通信工程學院, 黑龍江 哈爾濱 150001)

研究了基于多重信號分類(MUSIC)算法的特征分解算法,即關(guān)于Hermitian矩陣的特征分解.該算法的運算過程均是對低階矩陣進行運算,并且其中將復數(shù)運算轉(zhuǎn)化成了實數(shù)運算,不僅降低了運算的復雜度,而且所得到的特征值精度較高.該算法的此特性使得其易于在數(shù)字信號處理器(DSP)中實現(xiàn),應用于實際工程中具有較高的實時性.通過計算機仿真實驗以及DSP實現(xiàn)驗證了本算法的有效性以及實時性.

特征分解方法; Hermitian矩陣; 精度; 實時性

空間譜估計技術(shù)是陣列信號處理中的一個重要分支,其估計精度較高,可接近克拉美-羅界(Cramer-Rao Bound, CRB)[1].以多重信號分類(multiple signal classification, MUSIC)算法[2-5]和旋轉(zhuǎn)不變子空間(estimation of signal parameters via rotational invariance techniques, ESPRIT)算法[6-8]為代表的子空間分解類算法因其精確的角度估計和超分辨性能而被廣泛用于波達方向(direction of arrival, DOA)估計[9].其中,算法中所涉及到的噪聲子空間或信號子空間是通過對接收數(shù)據(jù)協(xié)方差矩陣進行特征分解所得到的.在實際工程中應用此算法進行DOA估計時,一方面要考慮估計角度的精度,另一方面需要考慮算法程序的運行時間.若在算法程序中的各個模塊都能保證運算的精度以及實時性,那么將很容易達到這兩點.特征分解作為此類算法中至關(guān)重要的部分,對其進行研究得到精度高并且實時性好的特征分解方法將能提高算法的測向效率.

“householder變換+一般QR”特征分解方法為在文獻[10]中記錄的方法的基礎上將其從實數(shù)域拓展到了復數(shù)域內(nèi),為先將Hermitian矩陣通過householder變換成對稱三對角矩陣,然后再對其使用QR方法進行特征分解[11];“householder變換+單步QR”特征分解方法也是如此,不同之處在于使用了位移的QR算法,能夠加快收斂速度;雅克比方法為參考文獻[12]將Hermitian矩陣轉(zhuǎn)化為實對稱矩陣,然后應用jacobi方法進行特征分解.本文所介紹的特征分解方法----二階主子陣實數(shù)化法與jacobi方法[13]相結(jié)合的特征分解方法不僅精度高而且實時性好.本算法的主要思想是在每個迭代中,將矩陣的一個二階主子陣用酉對角陣相似變換成二階的實對稱矩陣,然后再利用jacobi方法將此二階實對稱矩陣對角化.通過不斷循環(huán),不斷選取不同的二階主子陣,應用jacobi方法進行對角化,直到整個矩陣最終近似為對角陣[14-15].此種特征值分解方法與首先利用復Hermitian矩陣的實數(shù)部分與虛數(shù)部分構(gòu)成實對稱矩陣,然后應用jacobi方法對此實對稱矩陣進行特征分解求取特征值與特征向量的方法不同.本特征分解算法直接對復Hermitian矩陣進行特征分解,不需要將n階復Hermitian矩陣先轉(zhuǎn)換成2n階實對稱矩陣,如此便避免了在更高階矩陣上進行特征分解,減少了計算量,從而可以提高實時性.

1 算法原理

1.1 核心思想

設待特征分解的復Hermitian矩陣為A.此特征值分解方法的核心思想是得到合適的酉矩陣U使得UHAU=Λ=diag(λ1,λ2,…,λn).

令A1=A,依次構(gòu)造正交相似矩陣序列

(1)

其中,U(pk,qk)是平面(pk,qk)上的一個旋轉(zhuǎn)矩陣.得出酉矩陣U(pk,qk)的具體推導過程如下所示.

設在Hermitian矩陣A=B+i C中,ap q=bp q+icp q≠0,(p

(2)

根據(jù)反對角線上的元素作酉對角陣,有

(3)

在式(3)中,φ=angle(ap q),其中e-jφ=cosφ-jsinφ,故有

(4)

因而通過酉對角化,有

(5)

(6)

則有

(7)

1.2 具體公式推導

由上述所構(gòu)造的酉對角矩陣以及jacobi方法中的平面旋轉(zhuǎn)矩陣模型可令

(8)

則有

(9)

則在整個算法的矩陣變換過程中,所用到的旋轉(zhuǎn)矩陣U的表達式如式(10)所示(矩陣U(pk,qk)旁邊所標注的p和q表示的是第p行、q行和第p列、q列):

(10)

由式(1)中Ak與Ak+1的關(guān)系,根據(jù)文獻[16]所介紹的方法推導介紹從Ak到Ak+1的元素計算公式.

(11)

則可得到Ak+1中任意s行t列的元素為

(12)

在進行UHAkU=Ak+1運算之后,將所得到的矩陣Ak+1的元素與Ak的元素進行對比,同時又根據(jù)平面旋轉(zhuǎn)矩陣的性質(zhì),可以知道在Ak+1中,只有第p,q行與第p,q列發(fā)生了變化,因而只需計算Ak+1的第p,q行與第p,q列元素.Ak+1中的各元素的表達式如下(其中p

(13)

同樣的計算方式可以得到:

(14)

(15)

當i≠p,q時,則有

(16)

(17)

(18)

(19)

(20)

將式(20)兩邊同時除以cos2θ,可得

(21)

將式(21)進行變換得到

(22)

即有

(23)

通過式(23)則可以求得θ的值,進而求得cosθ以及sinθ的值.再由φ=angle(ap q),可由式(4)求得e-jφ=cosφ-jsinφ的值,從而可以得到酉矩陣U(pk,qk),將U(pk,qk)記為Uk,每次變換、迭代都有相應的旋轉(zhuǎn)矩陣.假設在m次迭代之后,使得原矩陣變換為對角陣,即為

(24)

(25)

(26)

(27)

2 MATLAB仿真分析

(1) 仿真條件:在MATLAB中,隨機產(chǎn)生20個8階的復Hermitian矩陣,對每個復Hermitian矩陣分別采用本文特征值分解方法以及MATLAB中的eig函數(shù)進行特征分解,然后分別將所得到的8個特征值從大到小排序并編號,將運用本文特征值分解方法所得到的特征值與相對應特征值順序的eig函數(shù)所產(chǎn)生的特征值進行相減,得到差值,即為絕對誤差,每個特征值編號對應一個絕對誤差.在將20個復Hermitian矩陣都用上述方式進行特征分解并求絕對誤差之后,對應每個特征值編號都有20個絕對誤差.然后再對每個特征值編號所對應的20個絕對誤差求均值,即為平均誤差,得到各個特征值序號對應的特征值的平均誤差如表1所示.

表1 特征值平均誤差

由表1可知,應用本文特征分解方法對Hermitian矩陣進行特征分解所得到的特征值的平均誤差的數(shù)量級在10-9左右,此數(shù)量級的誤差在實際工程中可以忽略不計,因而本文方法與Matlab中的特征分解函數(shù)eig函數(shù)的特征分解精度相當,在精度方面得到了保證.

(2) 將前述對本文特征值分解方法與“householder變換+一般QR”特征分解方法、householder變換+單步QR以及雅克比方法進行Matlab仿真實驗所得到的圖整理如圖1～圖4所示.

圖1 本文方法

比較圖1～圖4可以看到,運用本文特征分解方法對Hermitian矩陣進行特征分解所得到的特征值的精度最高,且遠遠高于其余幾種特征分解方法.

圖2 “householder+一般QR”

圖3 “householder+單步QR”

圖4 雅克比方法

3 實時性分析

3.1 理論運算量分析

假設待特征分解的hermitian矩陣為n階方陣.“householder變換與一般QR方法相結(jié)合”算法中(n-1)次householder變換共需要8n3(n-1)+(3n-2)次實數(shù)乘法,m次QR迭代共需要8mn3次實數(shù)乘法,計算得到特征向量矩陣共需要n3(m+n-1)次實數(shù)乘法,整個算法所需的運算量跟迭代次數(shù)相關(guān);而“householder變換與單步QR方法相結(jié)合”算法中的householder變換部分與前述算法的運算量一致,但是QR迭代部分在前述算法基礎上引入了單步位移,能夠加快收斂速度,因而算法運算量要少于前述算法,并且迭代次數(shù)是由設定門限值所決定的,故兩算法的運算量只能進行定性比較,而無法進行準確的定量比較;雅克比方法將對n階hemitian矩陣的特征分解轉(zhuǎn)化為對2n階實對稱矩陣的特征分解,也是設置了門限來確定整個迭代過程的迭代次數(shù),特征分解過程從復數(shù)域到實數(shù)域的變換雖然能夠在一定程度上減少算法的運算量,但矩陣的階數(shù)大大增加,又增加了計算量,并且算法總的運算量是與門限值所設置的大小是相關(guān)的,因而雅克比方法的運算量也無法進行定量統(tǒng)計;本文方法是遍歷選取非對角線元素非零的二階主子陣實數(shù)化為實對稱矩陣,再進行對角化,參與運算的都為二階矩陣,并且大部分為實值運算,因而整體的運算量較小,但由于算法最終的運算量與遍歷次數(shù)以及非對角線上元素為零的個數(shù)有關(guān),本文方法中所設置的遍歷次數(shù)為4次,已經(jīng)能夠達到較高的精度了,因而其運算量也無法進行準確的定量統(tǒng)計,只是從定性分析來說,本文方法的運算量最低.

3.2 硬件實現(xiàn)耗時對比分析

對本文特征分解方法進行DSP實現(xiàn), 分析其實時性,在本文中測試所用的DSP芯片為TMS320C6678. TMS320C6678是一款高性能的支持定點/浮點的多核DSP, 其總共含有8個DSP內(nèi)核,每個內(nèi)核的主頻可以達到1.25 GHz,支持16GFLOPs,而功耗僅為10 W.TMS320C6678的存儲器包括32 K字節(jié)的L1P存儲器/每核、32 K字節(jié)的L1D存儲器/每核、512 K字節(jié)的二級存儲器/每核. 在實際測試過程中,其每個內(nèi)核的主頻為1 GHz. 以其他三種特征分解方法為對比,進行實時性比較. 將此三種方法也進行了DSP實現(xiàn).通過使用軟件CCS 5.5上的計時工具, 分別得到本文以及其他三種特征值分解算法的C語言程序在TMS320C6678上運行的總周期個數(shù),然后轉(zhuǎn)化為以ms記時,如表2所示.

由表2中可以看到, 本文方法在DSP上實現(xiàn)時所耗費的時間最少,僅為0.348 405 ms,比其他幾種特征分解方法所耗費的時間要少得多,實時性最好.

表2 四種特征值分解方法在TMS320C6678上的運行時間

4 結(jié) 論

(1) 為在工程中應用子空間分解類算法對DOA進行快速精確估計,對特征分解這一部分進行仔細研究,得到本文中的二階主子陣與jacobi方法相結(jié)合的特征分解方法,該方法不僅精度高,而且實時性好;

(2) 由Matlab仿真結(jié)果可以看到,本文方法與eig函數(shù)的特征分解精度幾乎一致,滿足了工程上對于精度方面的要求;

(3) 由在TMS320C6678開發(fā)板上的測試結(jié)果可以看到,本文方法在眾多特征分解方法中實時性最高,算法程序運行時間很短,滿足了工程上對于實時性的要求;

(4) 該算法魯棒性好,性能穩(wěn)定,適用于各低階或高階的Hermitian矩陣的特征分解.

[ 1 ] YANG Z F,WANG Y Q,WANG L. Research and simulation of spatial spectrum estimation algorithm[C]∥2010 International Conference on Image Analysis and Signal Processing, 2010:532-536.

[ 2 ] 王偉,王曉萌,李欣,等. 基于MUSIC算法的L型陣列MIMO雷達降維DOA估計[J]. 電子與信息學報, 2014,36(8):1954-1959. (WANG W,WANG X M,LI X,et al. Reduced-dimensional DOA estimation based on MUSIC algorithm in MIMO radar with L-shaped array[J]. Journal of Electronics and Information Technology, 2014,36(8):1954-1959.

[ 3 ] 尤國紅,邱天爽,夏楠,等. 基于均勻圓陣的擴展循環(huán)MUSIC算法[J]. 通信學報, 2014,35(2):9-15. (YOU G H, QIU T S, XIA N,et al. Novel extended cyclic MUSIC algorithm based on uniform circular array[J]. Journal on Communications, 2014,35(2):9-15.)

[ 4 ] ZHANG X,CHEN C,LI J,et al. Blind DOA and polarization estimation for polarization-sensitive array using dimension reduction MUSIC[J]. Multidimensional Systems and Signal Processing, 2014,25(1):67-82.

[ 5 ] 劉帥,閆鋒剛,金銘,等. 基于四元數(shù)MUSIC的錐面共形陣列極化-DOA聯(lián)合估計[J]. 系統(tǒng)工程與電子技術(shù), 2016,38(1):1-7. (LIU S,YAN F G,JIN M,et al. Joint polarization-DOA estimation for conical conformal array based on quaternion MUSIC[J]. Systems Engineering and Electronics, 2016,38(1):1-7.)

[ 6 ] 梁浩,崔琛,余劍. 基于ESPRIT算法的十字型陣列MIMO雷達降維DOA估計[J]. 電子與信息學報, 2016,38(1):80-89. (LIANG H,CUI C,YU J. Reduced-dimensional DOA estimation based on ESPRIT algorithm in monostatic MIMO radar with cross array[J]. Journal of Electronics and Information Technology, 2016,38(1):80-89.)

[ 7 ] LIU Z,XU T. Source localization using a non-cocentered orthogonal loop and dipole (NCOLD) array[J]. Chinese Journal of Aeronautics, 2013,26(6):1471-1476.

[ 8 ] YANG H C. ESPRIT-based coherent source localization with forward and backward vectors[J]. IEEE Transactions on Signal Processing, 2010,58(12):6416-6420.

[ 9 ] 劉昕卓,吳迪,司偉建,等. Toeplize預處理及改進秩損估計器的解相干和解耦合方法[J]. 沈陽大學學報(自然科學版), 2015,27(5):405-409. (LIU X Z,WU D,SI W J,et al. The decorrelation and decoupling method by toeplize pretreatment and improving the rank loss estimator[J]. Journal of Shenyang University(Natural Science), 2015,27(5):405-409.)

[10] 李慶揚,王能超,易大義. 數(shù)值分析[M]. 5版. 北京:清華大學出版社, 2008:261-266. (LI Q Y,WANG N C,YI D Y. Numerical analysis[M]. 5thedition. Beijing: Tsinghua University Press, 2008:261-266.)

[11] 杜鵑,馮思臣. 復矩陣的Givens變換及其QR分解[J]. 成都理工大學學報, 2011,38(6):693-696. (DU J,FENG S C. Givens transform and QR decomposition of complex matrix[J]. Journal of Chengdu University of Technology, 2011,38(6):693-696.)

[12] 李小波,薛王偉,孫志勇. 一種求解復Hermite矩陣特征值的方法[J]. 數(shù)據(jù)采集與處理, 2005,20(4):403-406. (LI X B,XUE W W,SUN Z Y. High quality method for solving eigenvalues of complex Hermite matrix[J]. Journal of Data Acquisition and Processing, 2005,20(4):403-406.)

[13] 于正文,尹慶莉. 求特征值的Jacobi方法[J]. 山東科學, 2011,24(6):19-21. (YU Z W,YIN Q L. Jacobi method foreigenvalues calculation[J]. Shandong Science, 2011,24(6):19-21.)

[14] 征道生. Hermite矩陣特征值問題的2階主子陣實數(shù)化法[J]. 華東師范大學學報(自然科學版), 1996(3):1-6. (ZHENG D S. An algorithm to realification two order principal submatrix for Hermitian matrix eigenproblems[J]. Journal of East China Normal University(Natural Science), 1996(3):1-6.)

[15] 薛長峰. Jacobi型方法的一些研究[J]. 鹽城工學院學報, 2000,13(1):11-17. (XUE C F. Some studies of Jacobi method[J]. Journal of Yancheng Institute of Technology, 2000,13(1):11-17.)

[16] 劉婷婷,劉俊卿,張健楠. 基于Jacobi方法的Hermitian矩陣特征分解算法[J]. 電子科技, 2010,23(12):60-61. (LIU T T,LIU J Q,ZHANG J N. The eigendecomposition algorithm of Hermitian matrix based on the method of Jacobi[J]. Electronic Science and Technology, 2010,23(12):60-61.)

【責任編輯: 肖景魁】

Eigendecomposition Algorithm Based on Hermitian Matrix

ZengFuhong,SiWeijian,QuZhiyu

(College of Information and Communication Engineering, Harbin Engineering University, Harbin 150001, China)

The eigendecomposition algorithm is studied, which is about the eigendecomposition of the Hermitian matrix and based on MUSIC algorithm. Most of the operations about this algorithm are low order matrix operations and some plural operations are transformed to real operations. It can reduce the complexity of computation and the eigenvalues have high precision. The feature of the proposed algorithm makes it easier to conduct DSP implement and when it is applied to practical engineering; it has the high real-time performance. At last, through the computer simulation experiment and the DSP implementation, the effectiveness and real-time performance of the algorithm are verified.

eigendecomposition method; Hermitian matrix; precision; real time

2016-01-13

國家自然科學基金資助項目(61671168); 黑龍江省科學基金資助項目(QC2016085).

曾富紅(1993-),女,湖南衡陽人,哈爾濱工程大學碩士研究生; 司偉建(1971-),男,黑龍江哈爾濱人,哈爾濱工程大學教授.

2095-5456(2016)06-0511-05

TN 911.7

A