ESPRIT 算法广义逆矩阵求解的快速FPGA 实现

王卫江，张拓锋，蒋荣堃,2,3，李泽英，王晓华，谭志昕，薛丞博,2

（1. 北京理工大学 集成电路与电子学院，北京 100081；2. 北京理工大学 重庆创新中心，重庆 401120；3. 北京理工大学 重庆微电子中心，重庆 401332；4. 北京遥测技术研究所，北京 100094）

相较于多信号分类（multiple signal classification，MUSIC）算法[1]等其他空间谱估计算法，基于旋转不变子空间的信号参数估计（estimating signal parameters via rotational invariance techniques ，ESPRIT）算法由于省去了谱峰搜索这一步骤，减少了很大一部分计算量且仍旧能够保持良好的精确性和稳定性. 但由于ESPRIT 算法需要对复数矩阵进行特征值分解和广义特征值分解，依旧具有很大的计算量. 而此类高分辨测向系统一般都对实时性有较高要求，仅仅依靠软件实现可能难以满足现实应用的需求，因此使用硬件实现来完成系统的加速是时下空间谱估计领域一个热点[2-4].

ESPRIT 算法的基本原理是利用子空间的旋转不变性来求解信号的入射角度[5]，其中比较复杂的部分是利用最小二乘法来求解信号子空间的广义特征值，这一部分也可以转化成求解信号子空间矩阵的广义逆矩阵. 由于信号子空间矩阵为复数矩阵且阶数大，其广义逆矩阵的求解运算量较大，使用MATLAB 等软件很难达到实时性要求. 而FPGA 平台以其功耗低、计算速度快等优点使得其目前越来越多的被应用于矩阵的各种运算[6]，如定点、浮点运算、特征值求解、矩阵求逆、矩阵相乘等.

目前基于FPGA 的矩阵求逆设计大部分都只应用于非奇异矩阵，而针对一般矩阵的广义逆求解的硬件实现很少有人研究. 在现有的硬件广义逆求解研究中，SAHOO 等[7]提出了一种基于正交三角（QR）分解的广义逆矩阵求解算法，但该算法注重提高精确度，在速度上没有过高要求. 应俊等[8]提出使用基于坐标旋转数字计算机（coordinate rotation digital computer，CORDIC）算法来对奇异值（SVD）分解进行优化，从而降低硬件资源的消耗，但同样引入了其他计算单元，较为复杂. 陈晓东等[9]使用Cholesky 分解来计算矩阵求逆，该算法计算相较于SVD 分解和QR分解计算简单，但同样涉及大量开方运算，不利于硬件实现. 周金强等[10]对Cholesky 分解进行了改进，通过引入辅助矩阵的方法规避掉了开方运算，具有更好的实时性能，但只能处理共轭对称矩阵，局限性较大.

本文针对ESPRIT 算法中信号子空间矩阵的特点，提出了使用广义逆公式来构建一个一般复数矩阵的广义逆求解系统. 在FPGA 上进行了硬件实现，最后结合MATLAB 仿真结果对该系统的性能进行了分析. 经验证，该系统计算简单，相较于其他方法减少了硬件消耗，并且在保证结果精确度的同时有效减少了运算时间.

1 矩阵的广义逆求解

广义逆矩阵是逆矩阵的推广形式. 逆矩阵只针对非奇异矩阵，而奇异矩阵或者非方阵并不存在逆矩阵，因此引入广义逆矩阵的概念,在MATLAB 中可以使用函数B= pinv（A）来进行求解，结果返回一个与矩阵A转置同型的矩阵B且满足ABA=A，BAB=B，称矩阵B为矩阵A的广义逆.

1.1 广义逆矩阵的定义

1.2 广义逆矩阵的解法

对于矩阵的广义逆求解，一般有公式法、奇异值（SVD）分解法和正交三角（QR）分解法.

1.2.1 公式法

称为右逆；

如果A不满秩，则只能进行奇异值分解，使用SVD分解法.

1.2.2 SVD 分解法

因为酉矩阵的逆矩阵等于其共轭转置矩阵，对角矩阵的逆为对角线元素取倒数，则矩阵A的广义逆为

1.2.4 3 种方法的比较与选择

3 种方法中，公式法可以看成是对求解方程组Ax=b的推广，当矩阵A满秩（行满秩或列满秩）时都可以使用公式法计算广义逆矩阵；QR 分解也只应用在矩阵满秩的情况下，一般常用的QR 分解算法如Householder 变换、Givens 旋转和Schmidt 正交化等计算都较为复杂；当矩阵A不满秩时只能选用SVD 分解法通过奇异值求解广义逆矩阵，这一方法对于任何A∈Cm×n都适用. 但在计算矩阵奇异值分解时，随着矩阵的增大计算量会以指数型增长，使得系统的计算时间大大增加. MARKKANDAN 等[12]中对SVD分解和QR 分解算法的计算复杂度进行了详细的分析，QR 分解相较于SVD 分解计算更加简单.

在ESPRIT 算法中需要使用最小二乘法，最小二乘的本质就是一个超定方程组的求解问题，超定方程组的求解方式之一就是通过求解广义逆的形式，即求解信号子空间矩阵的广义逆. 由于信号子空间矩阵由若干个线性无关的信号子空间向量组成，因此必定为列满秩矩阵，因此可以使用式（4）直接进行求解. 该方法相较于SVD 分解和QR 分解复杂计算更少，更加适合用于硬件实现，因此本设计选择使用公式法求解广义逆矩阵.

2 公式法求解过程

公式法求解信号子空间矩阵广义逆步骤分为：

（1）输入信号子空间矩阵S∈Cm×n；

（3）计算 (SHS)-1， 即求解SHS的逆矩阵；

（4）矩阵 (SHS)-1与 矩阵SH相乘得到信号子空间矩阵的广义逆矩阵S+.

矩阵SH为矩阵S的共轭转置，因此相乘后的矩阵SHS为厄米特矩阵，其主对角线上的元素都是实数，矩阵中每一个第i行第j列的元素都与第j行第i列的元素的共轭相等. 针对方阵求逆的问题一般常用的方法为Cholesky 分解法、QR 分解法和LU 分解法.其中Cholesky 分解法虽然计算量较小，但只能对正定矩阵进行分解，适用范围过小； QR 分解虽然适用于任意型的矩阵，应用范围广泛，但 QR 分解的过程过于复杂，不易于实现硬件化；LU 分解虽然也对矩阵有一定要求，但处理的矩阵SHS必定为非奇异矩阵，满足LU 分解的充分条件，且计算格式简单，计算量小，易于硬件实现. MARKKANDAN 等[12]中同样对LU 分解的计算复杂度进行了分析，随着矩阵维度的增大，其所消耗的资源相较于QR 分解算法越来越小，因此本文使用LU 分解来求解矩阵SHS的逆矩阵.

2.1 LU 分解

设矩阵A∈Cn×n， 将矩阵A分解成一个单位下三角矩阵L和一个上三角矩阵U的乘积，即A=LU，该种分解被称为Doolittle 分解或LU 分解[13].

2.2 下三角矩阵求逆

使用LU 分解的矩阵求逆运算表达式为

求解出L和U之后，需要对单位下三角矩阵L和上三角矩阵U分别求逆，并求U-1和L-1的乘积. 采用初等变换法可以得到三角矩阵的求逆公式. 对于n阶下三角矩阵L，与n阶单位矩阵E组成增广矩阵（L|E）：

由于L、U矩阵都为可逆矩阵，根据矩阵转置的逆矩阵等于矩阵逆的转置矩阵这一性质：

可以先将上三角阵U转置成与L同型的下三角矩阵再求其逆矩阵，再将结果取转置就能得到U的逆矩阵，这样硬件实现时更加快速简洁.

2.3 矩阵乘法

3 FPGA 系统设计

基于上述公式法求解信号子空间矩阵广义逆的算法，将FPGA 硬件设计分为了3 个主要模块：矩阵乘法模块Ⅰ、LU 分解求逆模块和矩阵乘法模块Ⅱ.FPGA 系统的结构框图如图1 所示，第一个矩阵乘法模块接收初始信号子空间矩阵S的信息，计算共轭转置矩阵SH与原矩阵S的乘法并将结果发送到LU分解求逆模块；LU 分解求逆子模块将接收到的矩阵分解成L和U矩阵，对L和U矩阵分别求逆再相乘得到输入矩阵的逆矩阵，并将逆矩阵发送给第二个矩阵乘法模块；矩阵乘法模块Ⅱ计算逆矩阵(SHS)-1与矩阵SH相 乘，最终得到的结果就是广义逆矩阵S+.各模块中矩阵的存储使用的都是FPGA 中BRAM IP核的双口RAM，根据矩阵元素的行列信息来获得对应数据在RAM 中的地址，以此来完成数据的读写和运算.

广义逆求解系统设计的模块示意图如图1 所示.

图1 广义逆求解系统结构框图Fig. 1 Structure block diagram of generalized inverse solution system

在FPGA 中数据有定点数和浮点数两种表示方式. 浮点数能够表达的数值范围更广，但是与此同时计算的复杂度也更高并且误差会随着数据的增大而增大；定点数能够表示的数值范围相对较小但误差恒定不会随着数据的变化而变化. 尤其涉及到硬件实现时，浮点运算十分复杂、会占用大量的资源. 因此在本设计中的数据都使用定点数表示，数据使用二进制表示，包括符号位、整数位和小数位3 部分，涉及到的复数的加减乘除运算也都为定点运算.

3.1 矩阵乘法模块

矩阵乘法模块的主要运算单元由一系列乘累加组成，计算量较大、耗时也比较长. 在本次FPGA 硬件设计中，使用状态机来进行乘累加计算的控制，当状态机跳转到工作状态时，求解乘积矩阵AB的第i行j列元素，需要将矩阵A的第i行所有元素与矩阵B第j列的所有元素一一对应相乘并求和. 计算完毕后状态机将转入空闲状态直到下一个使能信号到来，开始计算乘积矩阵AB的下一个元素. 图2 为矩阵乘模块硬件结构.

图2 矩阵乘模块硬件结构Fig. 2 Hardware structure of matrix multiplication module

矩阵乘法模块Ⅰ的功能为计算信号子空间矩阵S的共轭转置矩阵SH与其自身的乘积. 由于矩阵SH是矩阵S的共轭转置，相乘后得到的矩阵SHS为厄米特矩阵，因此在计算时可以直接省去对应主对角线元素中虚部的计算过程，以此减少部分计算量，节约了相对应的资源. 该模块首先读取存储在RAM_S 和RAM_S_i 的信号子空间矩阵，其中RAM_S 存储信号子空间矩阵的实部，RAM_S_i 存储信号子空间矩阵的虚部. 经过乘累加计算后将得到的矩阵SHS的元素存储到RAM_SHS中.

3.2 LU 分解求逆模块

根据上文得到的式（12）来设计LU 分解求逆硬件实现的算法. 首先需要生成L、U矩阵的行列信息，从第一行第一列开始（i=1，j=1）计算矩阵L和U的第一个元素l11，u11. 计算完毕后列数加一继续计算第一行第二列的元素l12，u12，当列数增加到大于矩阵维数时行数加一，列数置1. 以此类推，直到行数也大于矩阵维数时便可以将矩阵三角分解后的L和U矩阵完全求解. 由于L和U分别是下三角矩阵和上三角矩阵，这两个矩阵对应位置上的元素一般总有一个是0，所以在计算元素前首先对当前的行列信息进行判断以此来减少不必要的计算：

图3 LU 分解求逆模块结构框图Fig. 3 Structure block diagram of LU decomposition inversion module

（1）如果行列相等（i=j），此时矩阵L元素li j=1，矩阵U元素uij不为0；

（2）如果行数大于列数（i＞j），此时矩阵L元素li j不为0，而矩阵U元素ui j=0；

（3）如果行数小于列数（i＜j），此时矩阵L元素li j=0， 而矩阵U元素uij不为0.

此外根据递推公式可以看出在计算矩阵L和U当前行列的元素时需要用到对应位置的原矩阵的元素以及此行列前求解的L、U矩阵的元素. 因此可以根据行列信息来产生地址值，如使用行数i乘以当前计算的元素位数s再加上列数j，以此来读取存储在BRAM 中的L、U矩阵和原矩阵SHS中所对应的元素.

4 仿真结果与分析

4.1 广义逆求解模块功能仿真

本次设计中软件程序的编写和仿真都使用MATLAB R2020a 在Intel(R) Core(TM) i5-3470 处 理 器 上完成，内存为8 GB、主频为3.20 GHz. 而系统的硬件编译和仿真都在Vivado 平台上进行，FPGA 选用的是Virtex-7 的xc7vx485tffg1157-1 开发板，工作频率为200 MHz. 使用MATLAB 随机生成用于测试的信号子空间矩阵，该矩阵为行列数分别为59 和15 的复矩阵.

4.2 FPGA 逻辑资源分析

广义逆求解模块的LUT（Look-Up-Table）、Block RAM、Registers 和DSP（digital signal processor）的资源占用情况如表1 所示. 根据表格可以看出LU 分解求逆模块需要的逻辑资源是最大的，而矩阵乘法模块需要的资源较少，仅为前者的1/4. 这是由于LU 分解求逆涉及到了大量的除法运算，相较于仅是大量乘累加运算的矩阵乘运算占用的资源更多.

表1 各模块所占用的逻辑资源Tab. 1 Logical resources of each module

4.3 运算时间分析

通过硬件仿真可以得出使用FPGA 的设计在求解矩阵广义逆时用时约为2.18 ms. 各模块具体的工作时间如表2 所示, 可以看出矩阵乘模块运行的时间较长，约占总体时间的80.7%.

表2 各模块运行时间对比Tab. 2 Operation time of each module

这是由于矩阵乘法模块需要计算大量的乘累加运算，需要较长的运算时间且矩阵越大需要的时间也就越长. 而在MATLAB 中使用pinv 函数对同样的矩阵进行广义逆求解，计算5 次取平均时间约为15.7 ms.可以得出当FPGA 工作在200 MHz 的时钟下时，求解广义逆矩阵的时间比软件端快了大约7.2 倍，可以节省大约86%的计算时间.

4.4 ESPRIT 算法结果分析

将FPGA 求出的广义逆矩阵重新输入到MATLAB 中完成ESPRIT 算法的后续步骤并且与在MATLAB 中使用pinv 函数求得信号子空间矩阵的广义逆的ESPRIT 算法最终的结果进行误差分析如表3，可以看出最终所得角度的平均误差大概为0.04°，满足工程的需求.

表3 ESPRIT 算法结果及误差分析Tab. 3 ESPRIT algorithm results and error analysis

4.5 与其他算法对比

应俊等[8]中使用了CORDIC 算法对矩阵进行SVD 分解，系统运行的总时间为547.6 μs. 使用行列数都为8 的输入矩阵在与此文献相同的工作频率250 MHz 下对本设计进行仿真，最终得到本系统运行的时间为139.4 μs. 可以看出本设计使用的算法相较于传统的SVD 分解法在计算广义逆矩阵上具有更好的时间性能.

周金强等[10]中使用了改进Cholesky 分解来求解广义逆矩阵. 使用行列数都为4 的信号子空间矩阵，对本设计进行仿真，在矩阵维度和工作频率都相同的条件下，两种算法的逻辑资源和系统运算时间如表4.

表4 逻辑资源与计算时间Tab. 4 Logical resources and computing time

由表4 可以看出本设计相较于改进的Cholesky算法消耗的逻辑资源更少但是计算时间更长. 但文献[10]中的算法仅能计算共轭对称矩阵，局限性较大，而本设计的算法能够计算任意矩阵的广义逆矩阵，具有更多的应用场景.

5 结 论

本文针对ESPRIT 算法中信号子空间矩阵的特点，利用广义逆矩阵的性质提出了一种新的矩阵广义逆求解算法. 这种算法比一般的QR 分解和SVD分解法计算更加简单，用时更短. 利用FPGA 技术对该广义逆矩阵求解系统进行了硬件实现，将整个系统分为了矩阵乘法模块Ⅰ、LU分解求逆模块和矩阵乘法模块Ⅱ共3 个模块，对每个模块分别进行了设计和测试，给出了总体设计框图. 并且使用Vivado对广义逆矩阵求解系统进行了硬件的仿真得到时序图和计算结果. 将其与MATLAB 仿真的计算结果和平均时间进行对比，发现使用FPGA 的设计在求解矩阵广义逆时可以节省约86%的计算时间，且最终所得角度的平均误差大概为0.04°. 此结果可以验证该硬件系统在满足一般工程精度需求的同时可以一定程度上减少计算所需时间，具有更好的实时性.