曾富紅, 司偉建, 曲志昱

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

?

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

曾富紅, 司偉建, 曲志昱

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

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

特征分解方法; Hermitian矩陣; 精度; 實(shí)時(shí)性

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

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

1 算法原理

1.1 核心思想

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

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

(1)

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

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

(2)

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

(3)

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

(4)

因而通過酉對(duì)角化,有

(5)

(6)

則有

(7)

1.2 具體公式推導(dǎo)

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

(8)

則有

(9)

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

(10)

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

(11)

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

(12)

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

(13)

同樣的計(jì)算方式可以得到:

(14)

(15)

當(dāng)i≠p,q時(shí),則有

(16)

(17)

(18)

(19)

(20)

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

(21)

將式(21)進(jìn)行變換得到

(22)

即有

(23)

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

(24)

(25)

(26)

(27)

2 MATLAB仿真分析

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

表1 特征值平均誤差

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

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

圖1 本文方法

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

圖2 “householder+一般QR”

圖3 “householder+單步QR”

圖4 雅克比方法

3 實(shí)時(shí)性分析

3.1 理論運(yùn)算量分析

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

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

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

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

表2 四種特征值分解方法在TMS320C6678上的運(yùn)行時(shí)間

4 結(jié) 論

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

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

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

(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雷達(dá)降維DOA估計(jì)[J]. 電子與信息學(xué)報(bào), 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 ] 尤國(guó)紅,邱天爽,夏楠,等. 基于均勻圓陣的擴(kuò)展循環(huán)MUSIC算法[J]. 通信學(xué)報(bào), 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ì)[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雷達(dá)降維DOA估計(jì)[J]. 電子與信息學(xué)報(bào), 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預(yù)處理及改進(jìn)秩損估計(jì)器的解相干和解耦合方法[J]. 沈陽大學(xué)學(xué)報(bào)(自然科學(xué)版), 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] 李慶揚(yáng),王能超,易大義. 數(shù)值分析[M]. 5版. 北京:清華大學(xué)出版社, 2008:261-266. (LI Q Y,WANG N C,YI D Y. Numerical analysis[M]. 5thedition. Beijing: Tsinghua University Press, 2008:261-266.)

[11] 杜鵑,馮思臣. 復(fù)矩陣的Givens變換及其QR分解[J]. 成都理工大學(xué)學(xué)報(bào), 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] 李小波,薛王偉,孫志勇. 一種求解復(fù)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]. 山東科學(xué), 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í)數(shù)化法[J]. 華東師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 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] 薛長(zhǎng)峰. Jacobi型方法的一些研究[J]. 鹽城工學(xué)院學(xué)報(bào), 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.)

【責(zé)任編輯: 肖景魁】

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

國(guó)家自然科學(xué)基金資助項(xiàng)目(61671168); 黑龍江省科學(xué)基金資助項(xiàng)目(QC2016085).

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

2095-5456(2016)06-0511-05

TN 911.7

A