CN110110285B - 一种用于FPGA的并行Jacobi计算加速实现方法 - Google Patents

Abstract

本发明公开了一种用于FPGA的并行Jacobi计算加速实现方法。n×n维矩阵的数据输入到FPGA中利用并行Jacobi计算进行旋转变换处理，处理单元初始化，对角处理单元计算旋转角度对应的符号集并输出给非对角处理单元，对角处理单元元素更新，非对角处理单元元素更新，处理单元间元素交换，在对每个处理单元的元素进行更新后，将更新后的处理单元之间的元素进行交换。本发明能使用较少的FPGA资源，提高FPGA内部计算处理性能，能够有效地提高特征值分解在FPGA上实现的效率，在实际工程中有较高的应用价值。

Description

一种用于FPGA的并行Jacobi计算加速实现方法

技术领域

本发明属于信号处理技术领域的一种FPGA内部数据处理方法，涉及了一 种用于FPGA的并行Jacobi计算加速实现方法。

背景技术

雷达、无线通信、图像处理等诸多领域的许多算法都需要计算矩阵的特征 值。例如，特征值的计算是子空间类DOA(Direction of Arrival，到达角度)估 计算法和PCA(Principal Component Analysis，主成分分析)算法的关键步骤。

目前有大量计算特征值的算法，例如QR算法、LU分解算法、代数法等。 代数法由于求根步骤的复杂度随矩阵维度大大上升，不适合大规模矩阵求特征 值，而LU分解算法只能用于可逆矩阵求特征值。而且，尽管QR算法比串行Jacobi 计算计算特征值的速度更快，但是已有学者证明Jacobi计算比QR算法更精确。 Jacobi计算是通过一系列旋转将矩阵逐渐变换为一个近似对角矩阵的过程，矩阵 的对角元素即矩阵特征值。此外，Jacobi计算由于对实对称矩阵进行特征值分解 有其固有的并行性，使并行Jacobi计算(Jacobi计算的一种并行实现方法)在特 征值分解的FPGA实现中得到了普遍应用。

目前，已有一些并行Jacobi计算的加速研究，但是大多加速方法未能做到 在一个CORDIC算法周期内实现并行Jacobi计算的一步。现有的近似Jacobi计 算虽然可以在一个CORDIC算法周期内实现并行Jacobi计算的一步，但由于是 近似旋转，总共需要的旋转次数会增加，因此效果不佳。此外，由于FPGA总 LUT(查找表)资源有限，现有算法实现时均并未考虑FPGA中LUT资源的消 耗量。

发明内容

针对上述背景技术中存在的问题，本发明要解决的是提出了一种用于FPGA 的并行Jacobi计算加速实现方法，对并行Jacobi计算方法在FPGA内部实现设 计更优异的计算处理效果的方案，解决了FPGA内部数据处理较慢、资源消耗 多的技术问题，达到在一个CORDIC算法周期内实现并行Jacobi计算的仅一步 的目标小，且该FPGA片上资源消耗更小。

为了实现上述目的，本发明采用以下步骤的技术方案：

1)处理单元初始化

n×n维矩阵的数据输入到FPGA中利用并行Jacobi计算进行旋转变换处理， 并行Jacobi计算中采用CORDIC(坐标旋转数字计算)算法进行平面旋转，平 面旋转中建立二维xy坐标系；

FPGA(现场可编程门阵列)内分设有多个处理单元，多个处理单元阵列排 布，每个处理单元与自身相邻的处理单元通过数据接口连接，进行数据交互和 元素交换，将执行并行Jacobi计算的n×n维矩阵中各元素按以下公式分配到处 理单元P_ij中：

其中，P_ij表示第i行第j列的处理单元，a_2i,2j表示n×n维矩阵中第2i行第2j 列的元素，n表示矩阵的维度；

并且，下标i＝j的处理单元P_ij为对角处理单元，否则为非对角处理单元；处 理单元P_ij中下标满足2i＝2j和2i-1＝2j-1的元素为对角元素，否则为非对角元素；

由于n×n维矩阵为实对称矩阵，按照上述处理获得分配仅保留右上部分，左 下部分和右上部分以对角线对称。

2)对角处理单元计算旋转角度2θ对应的符号集并输出给非对角处理单元

用以下公式迭代求出CORDIC算法旋转角度2θ对应符号集{d_2θ,k}， k＝1,2,...,N，迭代总次数与CORDIC算法的迭代总次数相同：

θ₀＝2θ，tan(φ_k-1)＝2^-(k-1)，k＝1,2,...,N

其中，k表示迭代次数的序数，N表示迭代总次数，N取为FPGA所采用的 数据位数；α_k表示第k次迭代的第一符号参数，β_k表示第k次迭代的第二符号 参数，θ₀表示旋转角度初始值(即2θ)，θ_k表示经k次迭代后的剩余旋转角度， φ_k-1表示第k-1次迭代的角度参数，d_2θ,k表示第k次迭代下旋转角度2θ对应的符 号；

具体地，在符号计算模块中，d_2θ,k通过α_k-1和β_k-1的符号位进行异或运算得 到(符号位相同则d_2θ,k为1，符号位相反则d_2θ,k为-1)。α_k-1通过移位运算得到 α_k-12^k-1，β_k-1通过移位运算得到β_k-12^k-1。若d_2θ,k为1，则α_k通过α_k-12^k-1和β_k-1进行 减法运算得到，β_k通过β_k-12^k-1和α_k-1进行加法运算得到；若d_2θ,k为-1，则α_k通过 α_k-12^k-1和β_k-1进行加法运算得到，β_k通过β_{k- 1}2^k-1和α_k-1进行减法运算得到。

迭代计算开始，与对角处理单元中的非对角元素所对应的初始旋转角度为θ， 计算为：

α₀＝2a_pq，β₀＝a_pp-a_qq

其中，a_pq、a_qp分别表示对角处理单元中初始包含的两个非对角元素，且 a_qp＝a_pq，a_pp、a_pp分别表示对角处理单元中初始包含的对角元素；α₀表示初始 的第一符号参数，β₀表示初始的第二符号参数；

对角处理单元的对角元素a_pp和a_qq通过减法运算得到β₀＝a_pp-a_qq，非对角元 素a_pq通过移位运算得到α₀＝2a_pq，设在并行Jacobi计算中，与当前对角处理单元 中的非对角元素所对应的旋转角度为θ，β₀和α₀作为初始值送入符号集计算模 块，而符号计算模块通过迭代求出旋转角度2θ对应的符号集{d_2θ,k}。

最后，对角处理单元将自身计算获得的旋转角度2θ对应符号集{d_2θ,k}输出到 与自身处于同一行和同一列的非对角处理单元；

3)对角处理单元元素更新

由步骤2)中每次迭代求得的d_2θ,k作为CORDIC算法中第k次迭代的旋转符 号，代替传统CORDIC算法每次迭代后计算旋转符号的步骤，对第一待旋转坐 标(2a_pq，a_pp-a_qq)执行CORDIC算法以旋转角度2θ进行平面旋转；

步骤2)所有迭代完成后，将最终的平面旋转结果乘以第一补偿因子，得到 旋转后的y坐标，即y₁＝2a_pqsin2θ+(a_pp-a_qq)cos2θ，第一补偿因子用以下公式求出：

其中，C₁表示第一补偿因子；

然后对对角处理单元中的对角线元素用以下公式更新，并将非对角线元素 置0：

其中，a'_pp、a'_qq表示对角处理单元中更新后的两个对角线元素，y₁表示第 一待旋转坐标旋转后的y轴坐标；

4)非对角处理单元元素更新

4.1)非对角处理单元P_ij接收来自两个对角处理单元P_ii、P_jj输出的符号集， 表示为

和

分别表示第k次迭代下旋转角度2θ_i和旋转角 度2θ_j对应的符号，用以下公式分别计算两个符号

和

从而获得了两 个符号集

和

其中，

和

分别表示旋转角度θ_i+θ_j和旋转角度θ_i-θ_j对应的符号； 2θ_i和2θ_j分别表示两个对角处理单元P_ii、P_jj的非对角元素对应的旋转角度的二倍 角；

具体地，两符号集通过对每一对符号

的异或运算和数据选择器确 定旋转角度θ_l-θ_m的符号集，异或运算结果为1则取

作为

否则取0作 为

通过对每一对符号

的同或运算和数据选择器确定旋转角度 θ_l+θ_m的符号集，同或运算结果为1则取

作为

否则取0作为

4.2)

有{-1,0,1}三种取值，用以下公式计算由两个符号集

和

的前

个符号所有可能构成的符号组合对应的第二、第三补偿因子的取 值，一个符号组合是由

个符号构成，以各个不同符号组合对应的第二、第三 补偿因子取值建立查找表数据，以前

符号中各符号绝对值为查找地址，用 Block Memory(块随机存储器)生成查找表，查找表的地址位数取

数据 的深度

其中，C₂表示为第二补偿因子，C₃表示为第三补偿因子；

由于CORDIC算法迭代次数超过

(

向上取整)时，第二、第三补偿 因子与1的差值已经小于2^-N+1，而N位有符号定点数精度最高为2^-N+1，因此剩余 第二、第三补偿因子可以直接视为1，即无需补偿。

4.3)对于非对角处理单元，将非对角处理单元包含的四个元素表示为

由求得的

作为CORDIC算法中第k次迭代的旋转符号，对第二待旋转 坐标

执行CORDIC算法以旋转角度θ_i-θ_j进行平面旋转，将 平面旋转结果乘以第二补偿因子，第二补偿因子取值由步骤4.2)的查找表进行 查表取得，得到旋转后的坐标，表示为：

其中，x₂和y₂分别表示第二待旋转坐标旋转后的坐标；

由求得的

作为CORDIC算法中第k次迭代的旋转符号，对第三待旋转 坐标

执行CORDIC算法以旋转角度θ_i+θ_j进行平面旋转，将 平面旋转结果乘以第三补偿因子，第三补偿因子取值由步骤4.2)的查找表进行 查表取得，得到旋转后的坐标，表示为：

其中，x₃和y₃分别表示第三待旋转坐标旋转后的坐标；

4.4)然后采用以下公式对非对角处理单元中的元素进行更新：

其中，

和

分别表示非对角处理单元包含的四个元素；

具体的，x₂和y₃通过加法运算和移位运算，得到如公式所示的更新值

x₃和y₂通过加法运算和移位运算，得到如公式所示的更新值

x₃和y₂通过减 法运算和移位运算，得到如公式所示的更新值

x₂和y₃通过减法运算和移位 运算，得到如公式所示的更新值

5)处理单元间元素交换

在对每个处理单元的元素进行更新后，与之对称的矩阵元素也更新为相同 值，将更新后的处理单元之间的元素进行交换：

5.A)针对对角处理单元中的对角元素进行交换

设当前对角处理单元P_ii包含对角元素

和

然后：

针对对角元素

i表示对角处理单元的行列序数，若i＝1，则对角元素

不变，若i＝2，则对角元素

的值更换为对角元素

的值，若i>2，则对角 元素

的值更换为对角元素

的值；

针对对角元素

若

则对角元素

的值换为

的值；若

则对角元素

的值更换为对角元素

的值；

5.B)针对对角处理单元中的非对角元素和非对角处理单元中的元素进行交 换，均采用以下方式更换位置：将对角处理单元中的非对角元素和非对角处理 单元中的元素进行移动位置，可以跨处理单元地移动到其他处理单元，使元素 的行下标与移动后处于相同行的步骤5.A)交换后的对角元素的行号相同，且元 素的列下标与移动后处于相同列的步骤5.A)交换后的对角元素的列号相同；

6)交换后Jacobi计算将n×n维矩阵所有对角处理单元中的非对角元素都经 过更新一次，返回步骤2)进行下一处理和更新，不断重复以上更新过程使n×n 维矩阵的非对角元素逐渐收敛到0，达到预设收敛精度后结束更新，并行Jacobi 计算结束。

所述的n×n维矩阵为天线阵列采集到的或者图像降维前的数据的协方差矩 阵，为实对称矩阵。

所述步骤1)中，若n×n维矩阵的n为奇数，即奇数维的矩阵，则通过添加 第n+1列和第n+1行来将矩阵拓展成偶数维的矩阵，将添加的第n+1列和第n+1行 的元素数值全部取0。

本发明在步骤2)计算的第k个符号提供给步骤3)和步骤4)CORDIC算 法进行第k次迭代，因此所述步骤2)、步骤3)、步骤4)同时进行。

CORDIC算法计算结果经过简单组合即得到各处理单元元素更新值，所有 处理单元的元素都可以同时进行。本发明方法实现并行Jacobi计算的一步耗费 的时间仅为一个CORDIC周期，相比现有方法耗时的三个CORDIC周期大大减 少了计算时间，提高了计算性能。

本发明针对的n×n维矩阵的数据为天线阵列采集到的数据在进行DOA估计 时所用的，或者图像数据使用PCA算法进行降维时所使用的协方差矩阵。

本发明的有益效果是：

本发明采用特殊设计的线性组合方法取代现有并行Jacobi计算当中的双边 旋转方法，结合利用旋转角度符号集和两符号集的组合取代现有CORDIC算法 中计算旋转符号的步骤，提高了并行Jacobi计算的并行性，减少了并行Jacobi 计算中每步的计算时间，能够在一个CORDIC周期内实现并行Jacobi计算的一 步。

本发明有效地提高了并行Jacobi计算在硬件上的实现速度，能在一个 CORDIC算法周期内实现并行Jacobi的一步，弥补了传统方法的不足，耗时仅 为传统方法的三分之一。

本发明能使用较少的FPGA资源，提高FPGA内部计算处理性能，能够有 效地提高特征值分解在FPGA上实现的效率，在实际工程中有较高的应用价值。

附图说明

图1为本发明实施例对角处理单元结构图；

图2为本发明实施例非对角处理单元结构图；

图3为本发明实施例处理单元阵列结构图；

图4为本发明实施例计算方法流程图。

图5为本发明实施例MUSIC(多信号分类)算法的功率谱函数结果图。

具体实施方式

以下结合附图和具体实施例对本发明的实施作如下详述。

本发明的FPGA实现结构主要分为对角处理单元和非对角处理单元，对角 处理单元结构如图1所示，非对角处理单元结构如图2所示。处理单元阵列结 构图如图3所示。算法执行流程如图4所示。

本发明实施例及其实施过程如下：

本实例具体实施过程在Xilinx Virtex-7 XC7VX690T FPGA芯片上实现，实 现具体采用的是四元天线阵列采集无人机发射的无线信号，信号入射方向为0 度。根据四元天线接收的四组数据计算得到的一个4×4的实对称协方差矩阵，表 示为A。

采用16位定点数，对于

求 特征值，具体包括以下步骤：

(1)处理单元初始化。将R_r中各元素按分配到处理单元P_ij中。每个处理单 元与相邻的处理单元通过数据接口连接。下标满足i＝j的处理单元为对角处理单 元，否则为非对角处理单元。下标满足2i＝2j和2i-1＝2j-1的矩阵元素为对角元素， 否则为非对角元素。

(2)对角处理单元计算旋转角度对应的符号集并输出给非对角处理单元。 设对角处理单元中包含的非对角元素为a_pq，a_qp，且a_qp＝a_pq。设对角处理单元中 包含的对角元素为a_pp，a_pp。令α₀＝2a_pq，β₀＝a_pp-a_qq。设与当前对角处理单元中 的非对角元素所对应的旋转角度为θ。用迭代求出CORDIC算法旋转角度2θ对 应符号集d_2θ,k，k＝1,2,...,16。迭代次数与CORDIC算法迭代次数相同，取当前系 统采用的数据位数16。

(3)对角处理单元元素更新。用求出补偿因子。由步骤(2)中求得的d_2θ,k作为CORDIC算法中第k次迭代的旋转符号，代替传统CORDIC算法每次迭代 后计算旋转符号的步骤，对(2a_pq，a_pp-a_qq)执行CORDIC算法旋转2θ，结果乘以 补偿因子，得到旋转后的y坐标，即y₁＝2a_pqsin2θ+(a_pp-a_qq)cos 2θ，对角处理单元 中的对角线元素更新。并对非对角线元素置0。

(4)非对角处理单元元素更新。非对角处理单元P_ij接收来自两个对角处理 单元P_ii、P_jj输出的符号集，表示为

用分别计算

和

有{-1,0,1}三种取值，用计算出符号集

前16个符号所有可能取值 组合对应的补偿因子的取值，以补偿因子取值为查找表数据，以前16符号集中 各符号绝对值为查找地址，用Block Memory生成查找表。由于CORDIC算法迭 代次数超过8时，补偿因子与1的差值已经小于2^-7，而8位数据精度最高为2^-7， 因此剩余补偿因子可以直接视为1，即无需补偿。查找表的地址位数取8，数据 深度2⁸。本实例查找表如表1所示。

表1补偿值查找表

令当前非对角处理单元包含的矩阵元素为

以

作为 CORDIC算法中第k次迭代的旋转符号。对

执行CORDIC算 法旋转θ_l-θ_m，补偿因子取值由查表取得，对结果乘以补偿因子，得到旋转后的 坐标。

以

作为CORDIC算法中第k次迭代的旋转符号。对

执行CORDIC算法旋转θ_l+θ_m，补偿因子取值由查表取得，对结果乘以补偿因子， 得到旋转后的坐标。

非对角处理单元中元素更新。

(5)处理单元间元素交换。在对每个处理单元的元素进行更新后，与之对 称的矩阵元素也更新为相同的值，将更新后的元素和其他处理单元的元素进行 交换。

然后返回步骤2、3、4进行新一轮的计算和更新。经过3次交换后Jacobi 计算将矩阵所有非对角元素都经过对角处理单元更新一次，重复多次以上更新 过程使矩阵的非对角元素逐渐收敛到0，达到用户预设收敛精度后结束更新，并 行Jacobi计算结束。

具体结果如下：

第1轮：

更新后为

元素交换后为

第2轮：

更新后为

交换后为

第8轮：

更新后为

可见矩阵的非对角元素已经达到收敛条件(虽然并行Jacobi计算是一个使 对角元素趋近于0的算法，但是由于实际实现中采用有限位的定点数来表示小 数，故非对角元素可达到0，但也引入了误差)，此时A₈对角线上的元素即求得 的特征值。将所求得特征值用于信号DOA(到达角度)估计算法中，从图5可 见，MUSIC(多信号分类)算法的功率谱函数在0度有个峰值，可见本发明实 现了正确的功能。

本实施例从运行时间、FPGA资源消耗两个方面给出本发明实际应用的性能。

运行时间：由于数据采用16位定点数，因此CORDIC算法内部迭代共16 次，考虑结果补偿，CORDIC算法周期为17个FPGA时钟周期，考虑并行Jacobi 每步之间需进行元素交换占用1个时钟周期，因此本发明的并行Jacobi计算加 速实现方法实现并行Jacobi的一步共需要18个时钟周期。在本例中，设置收敛 条件为协方差矩阵的非对角元素最大值的绝对值小于0.001。经过8次迭代达到 收敛条件，用时144个时钟周期。在本例使用的时钟频率为250M，用时0.576 微秒。

资源消耗：实现本例的Verilog程序在Vivado 2017.1软件平台上进行综合, 结果表明本例共消耗LUT(查找表)2360个，消耗REG(寄存器)688个，分 别占总资源的0.54％和0.79‰，可见，该设计仅占用少量的FPGA资源。

传统实现并行Jacobi的方法在对角处理单元需要使用CORDIC算法周期求 解旋转角度,接着用该对角处理单元求得的旋转角度先后使用两次CORDIC算法 更新对角处理单元元素,共需要3个CORDIC算法周期，在非对角处理单元需要 等待对角处理单元求解旋转角度，接着也需连续使用两次CORDIC算法更新非 对角处理单元元素，两次旋转的角度分别为：同一行的对角处理单元传递的旋 转角度和同一列的对角处理单元传递的旋转角度旋转。各处理单元并行工作， 实现并行Jacobi的一步共需要至少3个CORDIC算法周期。而本发明使用的 CORDIC算法全部并行工作，仅需一次CORDIC算法周期。本发明与传统方法 的处理单元处理过程对比如表2所示。

表2本发明与现有并行Jacobi方法的处理单元处理过程对比

由此可见，本发明具有显著提高传统方法特征值求解速度的优势，在实际 工程需要快速实现特征值分解时具有较高应用价值。

本案由熟悉本领域技术的人员根据说明书和附图内容作出的等效结构变换， 均包含在本发明的专利范围内。

Claims (5)

1.一种用于FPGA的并行Jacobi计算加速实现方法，其特征在于：包括以下方面：

1)处理单元初始化

n×n维矩阵的数据输入到FPGA中利用并行Jacobi计算进行旋转变换处理，并行Jacobi计算中采用CORDIC算法进行平面旋转，平面旋转中建立二维xy坐标系；

FPGA内分设有多个处理单元，多个处理单元阵列排布，每个处理单元与自身相邻的处理单元通过数据接口连接，进行数据交互和元素交换，将执行并行Jacobi计算的n×n维矩阵中各元素按以下公式分配到处理单元P_ij中：

其中，P_ij表示第i行第j列的处理单元，a_2i,2j表示n×n维矩阵中第2i行第2j列的元素，n表示矩阵的维度；

并且，下标i＝j的处理单元P_ij为对角处理单元，否则为非对角处理单元；处理单元P_ij中下标满足2i＝2j和2i-1＝2j-1的元素为对角元素，否则为非对角元素；

2)对角处理单元计算旋转角度2θ对应的符号集并输出给非对角处理单元，

用以下公式迭代求出CORDIC算法旋转角度2θ对应符号集{d_2θ,k}，k＝1,2,...,N，迭代总次数与CORDIC算法的迭代总次数相同：

θ₀＝2θ，tan(f_k-1)＝2^-(k-1)，k＝1,2,...,N

其中，k表示迭代次数的序数，N表示迭代总次数，N取为FPGA所采用的数据位数；α_k表示第k次迭代的第一符号参数，β_k表示第k次迭代的第二符号参数，θ₀表示旋转角度初始值，θ_k表示经k次迭代后的剩余旋转角度，f_k-1表示第k-1次迭代的角度参数，d_2θ,k表示第k次迭代下旋转角度2θ对应的符号；

最后，对角处理单元将自身计算获得的旋转角度2θ对应符号集{d_2θ,k}输出到与自身处于同一行和同一列的非对角处理单元；

3)对角处理单元元素更新

由步骤2)中每次迭代求得的d_2θ,k作为CORDIC算法中第k次迭代的旋转符号，对第一待旋转坐标(2a_pq,a_pp-a_qq)执行CORDIC算法以旋转角度2θ进行平面旋转；

步骤2)所有迭代完成后，将最终的平面旋转结果乘以第一补偿因子，得到旋转后的y坐标，即y₁＝2a_pqsin2θ+(a_pp-a_qq)cos2θ，第一补偿因子用以下公式求出：

其中，C₁表示第一补偿因子；

然后对对角处理单元中的对角线元素用以下公式更新，并将非对角线元素置0：

其中，a'_pp、a'_qq表示对角处理单元中更新后的两个对角线元素，y₁表示第一待旋转坐标旋转后的y轴坐标；

4)非对角处理单元元素更新

4.1)非对角处理单元P_ij接收来自两个对角处理单元P_ii、P_jj输出的符号集，表示为

和

分别表示第k次迭代下旋转角度2θ_i和旋转角度2θ_j对应的符号，用以下公式分别计算两个符号

和

从而获得了两个符号集

和

其中，

和

分别表示旋转角度θ_i+θ_j和旋转角度θ_i-θ_j对应的符号；2θ_i和2θ_j分别表示两个对角处理单元P_ii、P_jj的非对角元素对应的旋转角度的二倍角；

4.2)用以下公式计算由两个符号集

和

的前

个符号所有可能构成的符号组合对应的第二、第三补偿因子的取值，一个符号组合是由

个符号构成，以各个不同符号组合对应的第二、第三补偿因子取值建立查找表数据，以前

符号中各符号绝对值为查找地址，用块随机存储器生成查找表，查找表的地址位数取

数据的深度

其中，C₂表示为第二补偿因子，C₃表示为第三补偿因子；

4.3)对于非对角处理单元，将非对角处理单元包含的四个元素表示为

由求得的

作为CORDIC算法中第k次迭代的旋转符号，对第二待旋转坐标

执行CORDIC算法以旋转角度θ_i-θ_j进行平面旋转，将平面旋转结果乘以第二补偿因子，第二补偿因子取值由步骤4.2)的查找表进行查表取得，得到旋转后的坐标，表示为：

其中，x₂和y₂分别表示第二待旋转坐标旋转后的坐标；

由求得的

作为CORDIC算法中第k次迭代的旋转符号，对第三待旋转坐标

执行CORDIC算法以旋转角度θ_i+θ_j进行平面旋转，将平面旋转结果乘以第三补偿因子，第三补偿因子取值由步骤4.2)的查找表进行查表取得，得到旋转后的坐标，表示为：

其中，x₃和y₃分别表示第三待旋转坐标旋转后的坐标；

4.4)然后采用以下公式对非对角处理单元中的元素进行更新：

其中，

和

分别表示非对角处理单元包含的四个元素；

5)处理单元间元素交换

在对每个处理单元的元素进行更新后，将更新后的处理单元之间的元素进行交换：

5.A)针对对角处理单元中的对角元素进行交换

设当前对角处理单元P_ii包含对角元素

和

然后：针对对角元素

i表示对角处理单元的行列序数，若i＝1，则对角元素

不变，若i＝2，则对角元素

的值更换为对角元素

的值，若i>2，则对角元素

的值更换为对角元素

的值；针对对角元素

若

则对角元素

的值换为

的值；若

则对角元素

的值更换为对角元素

的值；

5.B)针对对角处理单元中的非对角元素和非对角处理单元中的元素进行交换，均采用以下方式更换位置：将对角处理单元中的非对角元素和非对角处理单元中的元素进行移动位置，使元素的行下标与移动后处于相同行的步骤5.A)交换后的对角元素的行号相同，且元素的列下标与移动后处于相同列的步骤5.A)交换后的对角元素的列号相同；

6)交换后Jacobi计算将n×n维矩阵所有对角处理单元中的非对角元素都经过更新一次，返回步骤2)进行下一处理和更新，不断重复以上更新过程使n×n维矩阵的非对角元素逐渐收敛到0，达到预设收敛精度后结束更新，并行Jacobi计算结束。

2.根据权利要求1所述的一种用于FPGA的并行Jacobi计算加速实现方法，其特征在于：所述步骤2)中在迭代计算开始，与对角处理单元中的非对角元素所对应的初始旋转角度为θ，计算为：

α₀＝2a_pq，β₀＝a_pp-a_qq

其中，a_pq、a_qp分别表示对角处理单元中初始包含的两个非对角元素，且a_qp＝a_pq，a_pp、a_pp分别表示对角处理单元中初始包含的对角元素；α₀表示初始的第一符号参数，β₀表示初始的第二符号参数。

3.根据权利要求1所述的一种用于FPGA的并行Jacobi计算加速实现方法，其特征在于：所述的n×n维矩阵为天线阵列采集到的或者图像降维前的数据的协方差矩阵，为实对称矩阵。

4.根据权利要求1所述的一种用于FPGA的并行Jacobi计算加速实现方法，其特征在于：所述步骤1)中，若n×n维矩阵的n为奇数，则通过添加第n+1列和第n+1行来将矩阵拓展成偶数维的矩阵，将添加的第n+1列和第n+1行的元素数值全部取0。

5.根据权利要求1所述的一种用于FPGA的并行Jacobi计算加速实现方法，其特征在于：所述步骤2)、步骤3)、步骤4)同时进行。

Images