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

曾富红, 司伟建, 曲志昱

(哈尔滨工程大学 信息与通信工程学院, 黑龙江 哈尔滨 150001)



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

曾富红, 司伟建, 曲志昱

(哈尔滨工程大学 信息与通信工程学院, 黑龙江 哈尔滨 150001)

研究了基于多重信号分类(MUSIC)算法的特征分解算法,即关于Hermitian矩阵的特征分解.该算法的运算过程均是对低阶矩阵进行运算,并且其中将复数运算转化成了实数运算,不仅降低了运算的复杂度,而且所得到的特征值精度较高.该算法的此特性使得其易于在数字信号处理器(DSP)中实现,应用于实际工程中具有较高的实时性.通过计算机仿真实验以及DSP实现验证了本算法的有效性以及实时性.

特征分解方法; Hermitian矩阵; 精度; 实时性

空间谱估计技术是阵列信号处理中的一个重要分支,其估计精度较高,可接近克拉美-罗界(Cramer-Rao Bound, CRB)[1].以多重信号分类(multiple signal classification, MUSIC)算法[2-5]和旋转不变子空间(estimation of signal parameters via rotational invariance techniques, ESPRIT)算法[6-8]为代表的子空间分解类算法因其精确的角度估计和超分辨性能而被广泛用于波达方向(direction of arrival, DOA)估计[9].其中,算法中所涉及到的噪声子空间或信号子空间是通过对接收数据协方差矩阵进行特征分解所得到的.在实际工程中应用此算法进行DOA估计时,一方面要考虑估计角度的精度,另一方面需要考虑算法程序的运行时间.若在算法程序中的各个模块都能保证运算的精度以及实时性,那么将很容易达到这两点.特征分解作为此类算法中至关重要的部分,对其进行研究得到精度高并且实时性好的特征分解方法将能提高算法的测向效率.

“householder变换+一般QR”特征分解方法为在文献[10]中记录的方法的基础上将其从实数域拓展到了复数域内,为先将Hermitian矩阵通过householder变换成对称三对角矩阵,然后再对其使用QR方法进行特征分解[11];“householder变换+单步QR”特征分解方法也是如此,不同之处在于使用了位移的QR算法,能够加快收敛速度;雅克比方法为参考文献[12]将Hermitian矩阵转化为实对称矩阵,然后应用jacobi方法进行特征分解.本文所介绍的特征分解方法----二阶主子阵实数化法与jacobi方法[13]相结合的特征分解方法不仅精度高而且实时性好.本算法的主要思想是在每个迭代中,将矩阵的一个二阶主子阵用酉对角阵相似变换成二阶的实对称矩阵,然后再利用jacobi方法将此二阶实对称矩阵对角化.通过不断循环,不断选取不同的二阶主子阵,应用jacobi方法进行对角化,直到整个矩阵最终近似为对角阵[14-15].此种特征值分解方法与首先利用复Hermitian矩阵的实数部分与虚数部分构成实对称矩阵,然后应用jacobi方法对此实对称矩阵进行特征分解求取特征值与特征向量的方法不同.本特征分解算法直接对复Hermitian矩阵进行特征分解,不需要将n阶复Hermitian矩阵先转换成2n阶实对称矩阵,如此便避免了在更高阶矩阵上进行特征分解,减少了计算量,从而可以提高实时性.

1 算法原理

1.1 核心思想

设待特征分解的复Hermitian矩阵为A.此特征值分解方法的核心思想是得到合适的酉矩阵U使得UHAU=Λ=diag(λ1,λ2,…,λn).

令A1=A,依次构造正交相似矩阵序列

(1)

其中,U(pk,qk)是平面(pk,qk)上的一个旋转矩阵.得出酉矩阵U(pk,qk)的具体推导过程如下所示.

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

(2)

根据反对角线上的元素作酉对角阵,有

(3)

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

(4)

因而通过酉对角化,有

(5)

(6)

则有

(7)

1.2 具体公式推导

由上述所构造的酉对角矩阵以及jacobi方法中的平面旋转矩阵模型可令

(8)

则有

(9)

则在整个算法的矩阵变换过程中,所用到的旋转矩阵U的表达式如式(10)所示(矩阵U(pk,qk)旁边所标注的p和q表示的是第p行、q行和第p列、q列):

(10)

由式(1)中Ak与Ak+1的关系,根据文献[16]所介绍的方法推导介绍从Ak到Ak+1的元素计算公式.

(11)

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

(12)

在进行UHAkU=Ak+1运算之后,将所得到的矩阵Ak+1的元素与Ak的元素进行对比,同时又根据平面旋转矩阵的性质,可以知道在Ak+1中,只有第p,q行与第p,q列发生了变化,因而只需计算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,每次变换、迭代都有相应的旋转矩阵.假设在m次迭代之后,使得原矩阵变换为对角阵,即为

(24)

(25)

(26)

(27)

2 MATLAB仿真分析

(1) 仿真条件:在MATLAB中,随机产生20个8阶的复Hermitian矩阵,对每个复Hermitian矩阵分别采用本文特征值分解方法以及MATLAB中的eig函数进行特征分解,然后分别将所得到的8个特征值从大到小排序并编号,将运用本文特征值分解方法所得到的特征值与相对应特征值顺序的eig函数所产生的特征值进行相减,得到差值,即为绝对误差,每个特征值编号对应一个绝对误差.在将20个复Hermitian矩阵都用上述方式进行特征分解并求绝对误差之后,对应每个特征值编号都有20个绝对误差.然后再对每个特征值编号所对应的20个绝对误差求均值,即为平均误差,得到各个特征值序号对应的特征值的平均误差如表1所示.

表1 特征值平均误差

由表1可知,应用本文特征分解方法对Hermitian矩阵进行特征分解所得到的特征值的平均误差的数量级在10-9左右,此数量级的误差在实际工程中可以忽略不计,因而本文方法与Matlab中的特征分解函数eig函数的特征分解精度相当,在精度方面得到了保证.

(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方法相结合”算法中(n-1)次householder变换共需要8n3(n-1)+(3n-2)次实数乘法,m次QR迭代共需要8mn3次实数乘法,计算得到特征向量矩阵共需要n3(m+n-1)次实数乘法,整个算法所需的运算量跟迭代次数相关;而“householder变换与单步QR方法相结合”算法中的householder变换部分与前述算法的运算量一致,但是QR迭代部分在前述算法基础上引入了单步位移,能够加快收敛速度,因而算法运算量要少于前述算法,并且迭代次数是由设定门限值所决定的,故两算法的运算量只能进行定性比较,而无法进行准确的定量比较;雅克比方法将对n阶hemitian矩阵的特征分解转化为对2n阶实对称矩阵的特征分解,也是设置了门限来确定整个迭代过程的迭代次数,特征分解过程从复数域到实数域的变换虽然能够在一定程度上减少算法的运算量,但矩阵的阶数大大增加,又增加了计算量,并且算法总的运算量是与门限值所设置的大小是相关的,因而雅克比方法的运算量也无法进行定量统计;本文方法是遍历选取非对角线元素非零的二阶主子阵实数化为实对称矩阵,再进行对角化,参与运算的都为二阶矩阵,并且大部分为实值运算,因而整体的运算量较小,但由于算法最终的运算量与遍历次数以及非对角线上元素为零的个数有关,本文方法中所设置的遍历次数为4次,已经能够达到较高的精度了,因而其运算量也无法进行准确的定量统计,只是从定性分析来说,本文方法的运算量最低.

3.2 硬件实现耗时对比分析

对本文特征分解方法进行DSP实现, 分析其实时性,在本文中测试所用的DSP芯片为TMS320C6678. TMS320C6678是一款高性能的支持定点/浮点的多核DSP, 其总共含有8个DSP内核,每个内核的主频可以达到1.25 GHz,支持16GFLOPs,而功耗仅为10 W.TMS320C6678的存储器包括32 K字节的L1P存储器/每核、32 K字节的L1D存储器/每核、512 K字节的二级存储器/每核. 在实际测试过程中,其每个内核的主频为1 GHz. 以其他三种特征分解方法为对比,进行实时性比较. 将此三种方法也进行了DSP实现.通过使用软件CCS 5.5上的计时工具, 分别得到本文以及其他三种特征值分解算法的C语言程序在TMS320C6678上运行的总周期个数,然后转化为以ms记时,如表2所示.

由表2中可以看到, 本文方法在DSP上实现时所耗费的时间最少,仅为0.348 405 ms,比其他几种特征分解方法所耗费的时间要少得多,实时性最好.

表2 四种特征值分解方法在TMS320C6678上的运行时间

4 结 论

(1) 为在工程中应用子空间分解类算法对DOA进行快速精确估计,对特征分解这一部分进行仔细研究,得到本文中的二阶主子阵与jacobi方法相结合的特征分解方法,该方法不仅精度高,而且实时性好;

(2) 由Matlab仿真结果可以看到,本文方法与eig函数的特征分解精度几乎一致,满足了工程上对于精度方面的要求;

(3) 由在TMS320C6678开发板上的测试结果可以看到,本文方法在众多特征分解方法中实时性最高,算法程序运行时间很短,满足了工程上对于实时性的要求;

(4) 该算法鲁棒性好,性能稳定,适用于各低阶或高阶的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 ] 尤国红,邱天爽,夏楠,等. 基于均匀圆阵的扩展循环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 ] 刘帅,闫锋刚,金铭,等. 基于四元数MUSIC的锥面共形阵列极化-DOA联合估计[J]. 系统工程与电子技术, 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] 李庆扬,王能超,易大义. 数值分析[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]. 数据采集与处理, 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阶主子阵实数化法[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