可配置Cholesky分解矩阵求逆的FPGA实现

胡铁乔，张毛毛，李阳波

（中国民航大学天津市智能信号与图像处理重点实验室，天津 300300）

可配置Cholesky分解矩阵求逆的FPGA实现

胡铁乔，张毛毛，李阳波

（中国民航大学天津市智能信号与图像处理重点实验室，天津 300300）

在阵列信号抗干扰处理算法中，功率倒置法作为简单有效的干扰抑制手段被广泛应用，本算法的关键是求出对应协方差矩阵的逆矩阵。鉴于FPGA的速度优势，目前阵列信号抗干扰处理算法多采用FPGA实现。本研究针对常用阵列信号抗干扰系统的实现，设计出一种基于FPGA的Cholesky矩阵分解模块，可适用于3～8阶正定Hermitian矩阵的求逆，配置简单方便。测试结果表明，在不影响计算精度的情况下，该方法节省了FPGA中的硬件资源。

Cholesky分解；矩阵求逆；阵列信号抗干扰；FPGA实现

目前，阵列天线抗干扰比较常用的算法包括功率倒置法、自适应波束形成法和盲自适应波束形成法等。由于功率倒置法不需知道信号的波达方向，只在干扰来向上形成零陷，因此该方法得到了广泛应用[1]。功率倒置法的关键是计算协方差矩阵的逆矩阵。

通过矩阵分解实现矩阵求逆的运算量小且易于实现，本文采用矩阵分解方法得到协方差矩阵的逆矩阵。常用的矩阵分解方法有：Cholesky分解、LU分解和QR分解。其中，Cholesky分解适用于共轭对称正定矩阵；LU分解适用于顺序主子式不为0的任意矩阵；QR分解适用于任意矩阵。在功率倒置法的算法中，协方差矩阵为共轭对称特性，所以采用Cholesky分解方法求解。

文献[2-3]采用基于Cholesky分解的矩阵求逆算法，其算法实现中使用较多的乘法器和除法器，这样大大提高了运算速度。文献[4]同样采用基于Cholesky分解的矩阵求逆算法，虽然也是可配置，但配置的矩阵维度较少，只能配置2、3、4三种维度的矩阵，且其算法在计算不同维度矩阵时消耗相同的资源。

目前，常用的阵列信号抗干扰处理系统使用的阵元数目多为4、5和7阵元，本文针对常用的阵列信号处理算法，对协方差矩阵求逆提出了基于FPGA的Cholesky分解，其适用于3～8阶的协方差矩阵求逆算法，鉴于协方差矩阵计算耗时较多，算法对时间的敏感度不强，于是采用分时复用计算，节省了硬件资源。

1 基于Cholesky分解求逆方法

根据Cholesky分解定理可得[4-5]：对角线元素为实数的对角矩阵D和对角线元素为1的下三角矩阵L及其共轭转置矩阵 LH，A=LDLH，即

其中：A为n×n正定Hermitian矩阵；W为A的逆矩阵；L为上三角矩阵；D为对角矩阵；C和Y均为下三角矩阵。

第1步由矩阵A分解得到矩阵L和矩阵C，令C=LD，即可得

由式（2）根据矩阵乘法公式可得，矩阵C的第一列与矩阵A的第一列元素相同，即得矩阵C的第一列。然后由矩阵C第一列的元素计算矩阵L的第一列元素，即

然后通过交叉求解按列由A求出矩阵C和矩阵L的其他列元素，具体计算公式为

第2步求矩阵C的逆矩阵Y。由于W为A的逆矩阵，AW=CLHW=E，令 LHW=Y，得 CY=E，且 Y为下三角矩阵，即

由式（6）根据矩阵乘法公式，按照对角线元素优先求取，可得矩阵Y的对角线元素为

然后计算与对角线元素平行元素依次求取，具体计算公式为

第3步由矩阵L和矩阵Y求矩阵W，即

由矩阵性质可知，共轭对称矩阵的逆矩阵也为共轭对称矩阵，即矩阵W为共轭对称矩阵，只需求解其一半元素即可获得矩阵W，从下至上依次计算，具体计算公式为

由以上3步即可完成一个正定Hermitian矩阵的求逆过程[6-7]。

2 可配置矩阵求逆的硬件实现

2.1 总体设计结构

本文提出的针对3～8阶正定Hermitian矩阵求逆的可配置Cholesky分解模块电路，包括矩阵输入、输出模块，矩阵运算模块和控制模块，其中矩阵运算模块包括矩阵A分解（矩阵C、L交叉计算）、矩阵C的求逆电路（矩阵Y的计算）、逆矩阵W的计算。根据式（4）、式（7）和式（10）的运算可知，矩阵的求逆过程需进行大量的加法、乘法和除法运算，而乘法器和除法器会占用大量的计算资源，所以在实现过程中需在“资源-速度”上做出权衡，系统采用了分时复用技术。硬件结构示意图如图1所示（图中实线箭头表示数据的计算流程，虚线箭头表示控制信号的流程）。

图1 硬件实现框图Fig.1 Hardware implementation block diagram

矩阵A的实部和虚部与使能信号进入Cholesky求逆系统，使能信号进入控制单元完成计数器清0和控制信号复位，矩阵求逆计算由控制单元控制矩阵A分解，交叉求解得到矩阵L和矩阵C；然后将矩阵L缓存，同时计算矩阵C的逆矩阵Y；矩阵Y计算完成后由控制单元控制矩阵Y和矩阵L计算得到逆矩阵W，最后输出逆矩阵和计算完成信号。

2.2 控制模块

该模块的功能是完成对Cholesky分解求逆模块的运算控制。本模块信号包括使能信号、输入矩阵维度、可配置Count_Chol计数器和控制字。图2是5×5矩阵的一个除法器的分时复用控制图。

图2 除法器分时复用控制图Fig.2 Divider time-sharing control chart

当使能信号输入，计数器Count_Chol从0开始计数。当计数器计数到某一数值时，产生信号对除法器进行赋值。控制字分为3种，分别为矩阵C和矩阵L交叉计算的控制字Scl_1、Scl_2…；矩阵C的逆矩阵Y计算的控制字Sy_1、Sy_2…；以及矩阵W计算的控制字Sw_1、Sw_2…。控制字符号的数字则代表计算的某一行或某一列，如Scl_1代表开始计算式（3）（即矩阵C的第二列），Sy_1则代表开始计算式（6）（即矩阵Y的对角线元素）。与此相次，复乘法器也通过相同计数器进行赋值。

以计算5阶矩阵的逆矩阵中除法器的赋值为例，如图2所示，根据式（3）和式（4）计算矩阵C的各列元素，当计数器 Count=Scl_1+divP（i-2）（divP 表示两次除法之间的间隔）时，即 Count=[0、5、10、15]时生成使能信号给除法器赋值c21/a11、c31/a11…；当Count=Scl_2+divP（i-3）时，即 Count=[20、25、30]时生成使能信号给除法器赋值c32/a22、c42/a22…，依此类推。根据式（7）计算矩阵Y的对角线元素，当Count=Sy_1+divP（i-2）时，即 Count=[75、80、85、90、95]时生成使能信号给除法器赋值1/c11、1/c22…，乘法器的赋值与此类似。

对于不同维度的矩阵，其求逆的计算过程相似，但计算的元素数目不同，通过改变控制字来控制计算过程。根据矩阵的阶数来计算不同阶矩阵的控制字，且每两个控制字之间相互联系。以5阶矩阵为例，介绍其中两个控制字之间的计算关系。Scl_3=Scl_2+divP（ChNO-1），其中 ChNO 为矩阵阶数，divP（ChNO-1）则表示计算矩阵C的第三列和矩阵L的第二列所需要的时延。Sy_1=Scl_3+divP（ChNO-2），其中ChNO为矩阵阶数。

不同维度矩阵的控制字数目也不同，随着矩阵的维度递减，需要的控制信号数目也递减。以8×8矩阵为例，整个计算过程需要21个控制字（Scl_1～Scl_6、Sy_1～Sy_8和 Sw_1～Sw_7），而对于 5× 5矩阵，整个计算过程则只需要12个控制信号（Scl_1～Scl_3、Sy_1～Sy_5和Sw_1～Sw_4），以此类推。这样做不仅简单灵活，而且节省了硬件资源。

2.3 矩阵运算模块

矩阵计算模块中采用分时复用的方法，由于系统处理的信号均为复信号，所以本系统占用了2个除法器（实部和虚部）和3个复乘模块，数据通过控制模块发出控制信号作为各个乘法器和除法器的控制输入使能。其中矩阵C、矩阵Y和矩阵W在计算过程中分别复用了一个复乘模块，矩阵L、矩阵C和矩阵Y在计算过程中共同复用了除法器模块。

由于FPGA计算浮点数时资源消耗较大，且计算速度较慢，所以在算法实现的过程中并保证精度的情况下采用移位为定点数计算。除法器和乘法器的运算结果均移位为整数。对于不同的系统要求，可通过修改移位位宽来满足不同的精度要求。

在算法实现过程中的复乘模块，要实现两个复数的共轭相乘，若用常规方法计算，根据

可知，该模块需占用4个乘法器和2个加法器。本文所用计算法为

采用修改后的方法，虽然加法器的数目增加了4个，但是乘法器只用了3个，这样虽增加了逻辑资源占用，但是节省了更多的运算资源。

对于不同阶矩阵，采用相同的计算模块。通过不同的控制字控制运算模块的赋值，这样可做到矩阵阶数越低，用时越短。若计算8×8矩阵的逆矩阵，需620个时钟周期，而如果计算3×3矩阵，则只需230个时钟周期。

3 资源占比及实验结果分析

完成算法设计后，验证算法的可行性，本文所使用的硬件编程语言为Verilog HDL，硬件编程环境为Xilinx ISE 13.4软件。并在Chip Scope Pro Analyzer中验证运算结果。

测试数据为7×7协方差矩阵，输入矩阵A数据位宽为32位，输出矩阵W数据位宽为20位。

3.1 资源占比

由于本方法采用的FPGA芯片为Xilinx Virtex-5 XC5VSX95T，此芯片含有58 880个片寄存单元（Slice Registers）和 58 880 个可编程逻辑单元（Slice LUTs），以及640个DSP核。表1为不同维度矩阵求逆所占用的资源对比。

表1 硬件资源占用情况对比（占用数目）Tab.1 Comparison of coccupied hardware resource

从表1中可以看出，计算的矩阵维度越少，系统占用的片寄存单元和可编程逻辑单元越少。

3.2 结果验证

结果验证方法采用FPGA硬件运算的结果与Matlab调用函数求逆的结果相比较。为了测试对比，inv_w为直接调用Matlab inv函数的结果，并扩大230倍的取整值，即

硬件计算结果采用按行输出，结果如图3所示，图中wr（1～7）_mux表示逆矩阵每一行的实部，wi（1～7）_mux表示逆矩阵每一行的虚部，Count_Chol即为控制计数器。对比inv_w和图3可以看出，硬件实现结果与Matlab计算结果相吻合。

图3 Cholesky矩阵求逆计算结果Fig.3 Calculation result of Cholesky matrix inversion

4 结语

本文基于阵列天线抗干扰功率倒置法的实现，综合考虑常用阵元数目，提出了一种基于Cholesky的协方差矩阵求逆算法，可灵活配置3～8阶协方差矩阵，降低了矩阵求逆的实现难度。在FPGA实现过程中，采用分时复用乘法器和除法器，节省了FPGA的计算资源。FPGA计算结果与Matlab直接调用函数计算结果相一致。另外，更多维度的矩阵求逆也可参照本算法实现。

[1]吴仁彪,卢 丹,李 婵.用于高动态GPS接收机的微分约束最小功率抗干扰方法[J].中国科学,2011,41(8):968-977.

[2]CHU Xuezheng,KHALED B,JOHN T.Rapid Prototyping of an Improved Cholesky Decomposition Based MIMO Detector on FPGAs[R].AHS,NASA/ESA Conference on,2009.

[3]魏婵娟,张春水,刘 健.一种基于Cholesky分解的快速矩阵求逆方法设计[J].电子设计工程,2014,22(1):159-164.

[4]潘 晓,许有云,甘小莺.基于Cholesky分解的可配置矩阵求逆FPGA实现[J].信息技术,2009,5(11):141-143.

[5]张贤达.矩阵分析与应用[M].北京:清华大学出版社,2004:225-227.

[6]李 涛,张忠培.矩阵求逆的FPGA实现[J].通信技术,2010,43(11):147-149.

[7]郭 磊.矩阵运算的硬件加速技术研究[D].长沙:国防科学技术大学,2010.

（责任编辑：杨媛媛）

FPGA implementation of matrix inversion based on configurable Cholesky decomposition

HU Tieqiao,ZHANG Maomao,LI Yangbo
（Intelligent Signal and Image Processing Key Lab of Tianjin,CAUC,Tianjin 300300,China）

As a simple and effective interference suppression algorithm,power inversion algorithm is widely used in various arraysignalprocessingalgorithms.Thekeyofthealgorithmistofindtheinversematrixofcorrespondingcovariance matrix.Given the high speed of FPGA,array signal anti-jamming processing algorithm is mostly accomplished by FPGA at present.In order to realize the commonly used system of array signal anti-jamming system,a Cholesky matrix decomposition based on FPGA is designed,which can be applied to the inversion of 3～8 positive definite Hermitian matrixes.Test results show that this algorithm can save hardware resource without affecting calculation precision.

Cholesky decomposition;matrix inversion;array signal anti-jamming;FPGA implementation

TN919；V243.1

：A

：1674－5590（2017）04－0007－04

2017-02-22；

：2017-04-09

：国家重点研发计划（2016YFB0502402）

胡铁乔（1970—），男，河南洛阳人，副教授，工学硕士，研究方向为基于FPGA、DSP高速处理平台的阵列信号处理系统.