文章快速检索     高级检索
  重庆邮电大学学报(自然科学版)  2020, Vol. 32 Issue (3): 434-440  DOI: 10.3979/j.issn.1673-825X.2020.03.013
0

引用本文 

应俊, 朱云鹏. 基于CORDIC矩阵奇异值分解的FPGA实现[J]. 重庆邮电大学学报(自然科学版), 2020, 32(3): 434-440.   DOI: 10.3979/j.issn.1673-825X.2020.03.013.
YING Jun, ZHU Yunpeng. FPGA implementation of matrix singular value decomposition based on CORDIC[J]. Journal of Chongqing University of Posts and Telecommunications (Natural Science Edition), 2020, 32(3): 434-440.   DOI: 10.3979/j.issn.1673-825X.2020.03.013.

基金项目

国家自然科学基金(61771082, 61871062, 61801065)

Foundation item

The National Natural Science Foundation of China(61771082, 61871062, 61801065)

作者简介

应俊(1976-),女,浙江杭州人,副教授,硕士,主要研究方向为电子新技术及其应用。E-mail:yingjun@cqupt.edu.cn; 朱云鹏(1992-),男,安徽合肥人,硕士研究生,主要研究方向为电力系统谐波检测。E-mail:yunpeng_zhu@foxmail.com

通讯作者

应俊 yingjun@cqupt.edu.cn.

文章历史

收稿日期: 2018-12-24
修订日期: 2019-12-16
基于CORDIC矩阵奇异值分解的FPGA实现
应俊1,2, 朱云鹏1,2     
1. 重庆邮电大学 光电工程学院,重庆 400065;
2. 重庆邮电大学 光通信与网络重点实验室,重庆 400065
摘要: 在矩阵的奇异值分解(singular value decomposition,SVD)过程中,随着矩阵维数的增加,SVD的计算量呈指数型增长,从而降低了算法运行的实时性。针对这个问题,基于Hestenes-Jacobi数值计算方法,提出了一种改进的基于坐标旋转数字计算机(coordinate rotation digital computer,CORDIC)的逻辑设计,该逻辑设计采用并行的全流水线设计思想,能够提高Jacobi平面旋转变换的运行速度,进而加快任意维矩阵奇异值分解的计算速度。分析了基于Hestenes-Jacobi方法的SVD的数值计算过程,介绍了CORDIC算法的基本原理,并具体说明了基于CORDIC算法的Jacobi平面旋转模块的设计,利用Verilog语言实现设计并验证,在现场可编程门阵列(field-programmable gate array,FPGA)上运行该逻辑设计单元,与Matlab软件的运行结果进行对比。实验测试结果表明,该结构能够减少计算时间,适应高速数据处理的要求。
关键词: 奇异值分解    现场可编程门阵列    单边雅克比变换    
FPGA implementation of matrix singular value decomposition based on CORDIC
YING Jun1,2 , ZHU Yunpeng1,2     
1. School of Optoelectronic Engineering, Chongqing University of Posts and Telecommunications, Chongqing 400065, P.R. China;
2. Key Laboratory of Optical Communication and Networks, Chongqing University of Posts and Telecommunications, Chongqing 400065, P.R. China
Abstract: In singular value decomposition (SVD) of a matrix, the matrix dimension increases along with the exponential increase in the computational complexity of SVD; that reduces the real-time performance of the algorithm. To alleviate this issue, this paper proposes an improved logic design for the coordinate rotation digital computer (CORDIC) based on the Hestenes-Jacobi numerical calculation method. The logic design adopts the parallel full pipeline design idea that improves the running speed of the Jacobi plane rotation transformation. Therefore, to accelerate the computational speed of singular value decomposition of an arbitrary dimensional matrix; this paper first analyzes the numerical calculation process of SVD based on Hestenes-Jacobi method, then introduces the basic principle of CORDIC algorithm, and specifies the design of Jacobi plane rotation module based on CORDIC algorithm. Furthermore, it also uses Verilog language to realize the design and verification of the logic. To validate the concept the logic design unit is implemented on the field programmable gate array (FPGA) and Matlab toolkit. The experimental test results show that the structure can reduce the calculation time and meet the requirements of high-speed data processing.
Keywords: singular value decomposition    field programmable gate array    one-sided Jacobi    
0 引言

奇异值分解(singular value decomposition, SVD)是矩阵分析中对矩阵进行对角化的数值方法,是线性代数中最有效且最有用的工具之一。在许多科学和工程领域中,SVD的计算都是必不可少的。在通信信号处理、信号与图像处理以及应用越来越广泛的大数据应用中,SVD算法都有着相当重要的研究与应用价值[1-4]

在过去的几十年中,人们一直致力于提高SVD的处理速度[5-6]和减少设计的资源消耗[7-8],如何快速高效地实现SVD已成为了一个重要的研究方向。文献[9]提出在图形处理器(graphics processing unit,GPU)上实现SVD,以提高处理速度。但是,利用GPU实现SVD存在线程同步和内存读取等问题,这会导致消耗大量的能量。文献[10]中提出了一种基于脉动阵列的架构用于SVD。但是,这种架构调度复杂,且仅适用于方阵,资源利用效率不高。

本文在已有研究的基础上,提出了一种基于坐标旋转数字计算机(coordinate rotation digital computer, CORDIC)算法的改进SVD硬件架构。该架构实现了基于定点数的Hestenes-Jacobi算法,可用于任意维矩阵的SVD。通过对Hestenes-Jacobi方法中旋转变换过程的分析,结合全流水的CORDIC运算单元,对该方法的实现过程进行改进,在提高处理速度的同时降低了资源消耗。

1 理论基础 1.1 奇异值分解理论

一个m×n的矩阵的奇异值分解可以表示为

$ {\mathit{\boldsymbol{A}}_{m \times n}} = {\mathit{\boldsymbol{U}}_{m \times m}} \times {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{m \times n}} \times \mathit{\boldsymbol{V}}_{n \times n}^{\rm{T}} $ (1)

(1) 式中:UV都是正交矩阵,它们分别满足UUT=IVVT=II是单位矩阵;Σ是对角矩阵,其对角线上的元素σi(i=1, 2…n)为矩阵A的奇异值,为非负数。

1.2 Hestenes-Jacobi方法

Hestenes发现,通过利用平面旋转变换的方法使2个向量正交与将矩阵的元素置为零这二者之间存在一定的等价关系[11]。在基于此项研究发现的基础上,Hestenes提出了一种单边Jacobi方法,也被称作Hestenes-Jacobi方法。

单边Jacobi方法的核心思想是通过对维度为m×n的矩阵A采用一系列的Jacobi平面旋转变换[12],对其进行正交化

$ \mathit{\boldsymbol{B}} = \mathit{\boldsymbol{A}}({\mathit{\boldsymbol{J}}_1}{\mathit{\boldsymbol{J}}_2}{\mathit{\boldsymbol{J}}_3} \cdots ) $ (2)

使得矩阵B中任意2列向量都满足关系biTbj=0。

Jacobi平面旋转变换如(4)式矩阵的结构

$ \begin{array}{l} \mathit{\boldsymbol{J}}(i,j,\theta ) = \left[ {\begin{array}{*{20}{c}} 1& \cdots &0& \cdots &0& \cdots &0\\ \vdots &{}& \vdots &{}& \vdots &{}& \vdots \\ 0& \cdots &c& \cdots &{ - s}& \cdots &0\\ \vdots &{}& \vdots &{}& \vdots &{}& \vdots \\ 0& \cdots &s& \cdots &c& \cdots &0\\ \vdots &{}& \vdots &{}& \vdots &{}& \vdots \\ 0& \cdots &0& \cdots &0& \cdots &0 \end{array}} \right]\begin{array}{*{20}{c}} {}\\ {}\\ i\\ {}\\ j\\ {}\\ {} \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{c}} \ \ \ \ \ \ \ \ \ \ \ \ {}&i&{}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j}&{}&{} \end{array} \end{array} $ (3)

(3) 式中:(i, j)表示消去元素在矩阵中的位置;c表示cosθs表示sinθθ称为旋转角。不难看出,JTJ=I。因此,Jacobi平面旋转变换为正交变换。因为Jacobi平面旋转变换矩阵本身即正交矩阵,令V=J1J2J3…,因此,V也是正交矩阵,即VT=V-1

正交化完毕后,对得到的矩阵B进行归一化,即

$ \mathit{\boldsymbol{B}} = \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }} = \mathit{\boldsymbol{U \boldsymbol{\varSigma} }} $ (4)

(4) 式中,Σ=diag(σ1, σ2σn-1, σn),σi2=biTbi

所以,矩阵A的奇异值分解可表示为

$ \mathit{\boldsymbol{A}} = \mathit{\boldsymbol{U \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{T}}} $ (5)

在对矩阵A进行正交化时,每一次的旋转变换过程中只有ij 2列元素受到影响。因此,一次正交变换可表示为

$ \left[ {\begin{array}{*{20}{c}} {{b_{i1}}}&{{b_{j1}}}\\ {{b_{i2}}}&{{b_{j2}}}\\ \vdots & \vdots \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{a_{i1}}}&{{a_{j1}}}\\ {{a_{i2}}}&{{a_{j2}}}\\ \vdots & \vdots \\ {{a_{im}}}&{{a_{jm}}} \end{array}} \right]\left( {\begin{array}{*{20}{c}} {{\rm{cos}}\theta }&{ - {\rm{sin}}\theta }\\ {{\rm{sin}}\theta }&{{\rm{cos}}\theta } \end{array}} \right) $ (6)

对2列向量进行正交化,即向量biTbjT正交

$ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{b}}_i^{\rm{T}}{\mathit{\boldsymbol{b}}_j} = (c\mathit{\boldsymbol{a}}_i^{\rm{T}} + s\mathit{\boldsymbol{a}}_j^{\rm{T}})(c{\mathit{\boldsymbol{a}}_j} - s{\mathit{\boldsymbol{a}}_i}) = }\\ {({c^2} - {s^2})\mathit{\boldsymbol{a}}_i^{\rm{T}}{\mathit{\boldsymbol{a}}_j} + cs(\mathit{\boldsymbol{a}}_j^{\rm{T}}{\mathit{\boldsymbol{a}}_j} - \mathit{\boldsymbol{a}}_i^{\rm{T}}{\mathit{\boldsymbol{a}}_i}) = 0} \end{array} $ (7)

经等式变换后,即可得到

$ \begin{array}{l} \frac{{{\rm{sin}}\theta {\rm{cos}}\theta }}{{{\rm{co}}{{\rm{s}}^2}\theta - {\rm{si}}{{\rm{n}}^2}\theta }} = \frac{1}{2}{\rm{tan}}2\theta = \frac{{\mathit{\boldsymbol{a}}_i^{\rm{T}}{\mathit{\boldsymbol{a}}_j}}}{{\mathit{\boldsymbol{a}}_j^{\rm{T}}{\mathit{\boldsymbol{a}}_j} - \mathit{\boldsymbol{a}}_i^{\rm{T}}{\mathit{\boldsymbol{a}}_i}}}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{1}{2}{\rm{ta}}{{\rm{n}}^{ - 1}}\left( {\frac{{2\mathit{\boldsymbol{a}}_i^{\rm{T}}{\mathit{\boldsymbol{a}}_j}}}{{\mathit{\boldsymbol{a}}_j^{\rm{T}}{\mathit{\boldsymbol{a}}_j} - \mathit{\boldsymbol{a}}_i^{\rm{T}}{\mathit{\boldsymbol{a}}_i}}}} \right) \end{array} $ (8)

根据(8)式,可求得旋转角度值,即可通过该角度值计算相应的正余弦函数值,对矩阵列向量对元素实现正交变换。

1.3 CORDIC算法

CORDIC算法的思想最早在1959年被提出,其主要是通过将初始向量按照一个特定的角度序列不断旋转迭代以逼近预期的目标向量,每次旋转变换之后,对坐标值进行更新[13]。当迭代次数足够多时,就认为当前向量的横纵坐标值是所求目标角度的正余弦函数值输出,图 1即为平面坐标旋转模型。

图 1 平面旋转变换 Fig.1 Plane rotation transformation

图 1,直角坐标系中的点(x0, y0)逆时针转动θ角度,即可得另一个点(x1, y1),两坐标点的关系为

$ \begin{array}{*{20}{l}} {{x_1} = {x_0}{\rm{cos}}\theta - {y_0}{\rm{sin}}\theta }\\ {{y_1} = {x_0}{\rm{sin}}\theta + {y_0}{\rm{cos}}\theta } \end{array} $ (9)

通过提取并去除因数cosθ,方程可写成下面的形式,即

$ \begin{array}{*{20}{l}} {{{\hat x}_1} = {x_0} - {y_0}{\rm{tan}}\theta }\\ {{{\hat y}_1} = {y_0} + {x_0}{\rm{tan}}\theta } \end{array} $ (10)

(9) 式的变换被称为伪旋转变换,其旋转的角度是正确的,但是变换前后向量的模值增加了cos-1θ倍。

CORDIC方法的核心是旋转角度θ的变换[14],通过将旋转变换的角度θ分解为一系列旋转角度ai的和,即$\sum\limits_{i = 0}^{n - 1} {{d_i}} {\alpha _i}\left( {{d_i} = \pm 1, i = 0, 1, 2, \cdots n} \right)$,其中,角 度ai =arctan(2-i),完成角度θ的旋转逼近。引入第3个方程,称为角度累加器,用来在每次迭代过程中追踪累加的旋转角度。在每次旋转后,(10)式的变换关系变为

$ \left\{ {\begin{array}{*{20}{l}} {{{\hat x}_{i + 1}} = {x_i} - {d_i}{y_i}{2^{ - i}}}\\ {{{\hat y}_{i + 1}} = {y_i} + {d_i}{x_i}{2^{ - i}}}\\ {{z_{i + 1}} = {z_i} - {d_i}{a_i}} \end{array}} \right. $ (11)

(11)式中,符号di是一个判决算子,其值为±1,用于确定旋转的方向是顺时针还是逆时针。判决算子值的确定见表 1。选择不同的工作模式,即判决算子的确定方式,可以得到不同的函数值的计算。

表 1 CORIC算法的操作模式 Tab.1 CORDIC algorithm mode of operation

将CORDIC变换中的起点(x0, y0)设为(aiT, ajT)中的元素,终点(x1, y1)设为(biT, bjT)中的元素时,对比(6)式和(9)式可发现,CORDIC变换可以实现Jacobi平面旋转变换,因此,可以用基于CORDIC变换的硬件来实现Jacobi平面旋转变换,实现矩阵奇异值分解。

2 奇异值分解架构设计

本文奇异值分解的Hestenes-Jacobi算法具体流程如算法1。Hestenes-Jacobi方法是一种迭代方法,通过对每次迭代后任意2列向量的内积值,判断算法的收敛性,当达到设定的收敛阈值时,即认为算法达到收敛,TOL为设定的收敛阈值。由于矩阵A中任意2列向量进行正交化时,与其余列向量之间并无数据依赖关系。因此,算法1中Addr函数为控制多对列向量索引的生成,实现并行运算。Addr函数以后步骤即为正交化过程,根据(8)式求解正交化过程的角度值。利用求得的角度值,根据(6)式,对列向量对进行正交化运算。

算法1   Hestenes-Jacobi Method。

Input: matrix A

Output: UV

U=A;

V=I;

While con>TOL do

  for m=1→n-1 do

    Addr(a, b);

    for k=0→(n/2)-1

      i=a[k]; j=b[k]

      α=Ui*Ui;

      β=Uj*Uj;

      γ=Ui*Uj;

      θ=0.5tan-1(2γ/(β-α));

      Ui=Ui*cosθ-Uj*sinθ;

      Uj=Ui*sinθ+Uj*cosθ;

      Vi=Vi*cosθ-Vj*sinθ;

      Vj=Vi*sinθ+Vj*cosθ;

      con=max(con, abs(γ));

    end for

  end for

end while

其中,输入变量即为需要分解的矩阵A,首先计算旋转角度θ,经过多次正交旋转后,得到输出值,分别为左奇异矩阵U及右奇异矩阵V

2.1 奇异值分解

具体的奇异值分解实现电路的顶层设计如图 2图 2中的COM_Top模块为奇异值分解的Hestenes-Jacobi正交化变换单元,实现列对向量的正交化运算过程,本例设计了4个COM_Top模块,进行并行运算。电路采用2个双端口RAM单元,分别用于存储左奇异向量和右奇异向量,Addr单元即对应算法1中Addr函数,并行读取矩阵A的列向量,输出作为RAM单元的地址。Converge模块为收敛判决模块,用于判断运算是否达到预期的收敛精度。Norm模块为归一化模块,用于对正交变换后的矩阵进行归一化,依次输出奇异值。状态机作为控制单元,提供同步和控制信号,在整个分解过程中每个模块执行不同的任务。

图 2 奇异值分解电路顶层设计 Fig.2 Top-level design of singular value decomposition circuit

Hestenes-Jacobi方法的核心思想是实现一系列的列向量对的正交化变换。因此,本文主要是在对COM_Top模块以及其子模块的具体实现进行改进。因此,下面针对COM_Top的具体设计进行说明。COM_Top模块由3个部分组成:向量乘积运算模块,旋转角度运算模块以及旋转变换运算模块,如图 3

图 3 COM_Top模块硬件结构 Fig.3 COM_Top module hardware structure

COM_Top模块的运算过程如下:当初始矩阵数据存入双口RAM中完毕后,数据在状态机的控制下进入向量乘积运算运算单元,计算求得2个列向量的内积值,将向量的内积值输出至旋转角度运算单元。根据(8)式可知,旋转变换矩阵的三角函数角度值可通过向量的3个内积值计算得到。在计算得到角度值θ后,直接将角度值以及需要正交化的列向量送入旋转变换运算单元进行正交化。

输入分别为矩阵A的列向量以及单位向量。通过对单位向量的旋转变换对应式V=J1J2J3…,最后得到矩阵V。由于旋转变换过程相同,因此,本文仅以对矩阵A的变换过程为例进行说明。

2.2 向量乘积运算单元

图 4中,向量乘积模块主要是对矩阵列对向量的内积进行运算,分别为aiTaiaiTajajTaj

图 4 向量乘积运算模块硬件结构 Fig.4 Vector_Mul module hardware structure

当状态进入向量乘积运算时,Addr单元为双端口RAM生成相应的地址,以向量乘积运算提供操作数据,其提供的数据分别是矩阵A的第i列和第j列向量。在状态机的控制下,依次送入列向量的每个元素进行乘积和累加运算。在一组列向量完成运算后,输出内积值。

2.3 CORDIC运算单元

图 5为流水线方案的CORDIC硬件架构图。为达到设置流水线方案,在每一级迭代运算结束后插入寄存器单元,寄存中间的运算结果,供下一级运算使用,以提高整个模块的运行效率。CORDIC算法的用处很多,可以应用在众多运算单元中。

图 5 CORDIC模块硬件结构 Fig.5 CORDIC module hardware structure

COM_Top模块中旋转角度运算模块和旋转变换运算模块则是以CORDIC算法为基础,如表 1,根据具体的需求,可以设计成相应的运算单元。

旋转角度运算模块作用为计算正交化变换的角度值θ。本文通过利用CORDIC算法实现求解反正切函数运算单元,求解角度值。根据表 1,具体的设置如下:当选择操作模式为矢量模式时,即设置角度z0的值为0,x0的值设置为aiTai-ajTajy0的值设置为2aiTaj,即可实现求解角度值θ

旋转变换运算模块即为旋转变换运算单元,实现列向量对的正交化。由(6)式和(9)式可知,列向量的正交化变换过程与CORDIC算法的旋转变换运算过程保持一致。因此,针对旋转变换运算模块,本文采用CORDIC运算单元进行实现。由表 1可知,具体的设置如下:选择操作模式为旋转模式,设置初始角度z0θx0y0分别为(ij)列向量对中对应的元素。通过对列向量中所有元素进行运算,即可实现(ij)列向量对的正交化。

通过利用全流水线方案的CORDIC运算单元实现旋转角度运算模块以及旋转变换运算模块,避免了算法1中求解正余弦函数之后的乘法运算过程,降低了资源消耗。同时,由于减少了乘法运算过程,大大降低了算法的时延,提高了算法的运行效率。

3 仿真分析

本文分别在中央处理器(central processing unit,CPU)和FPGA 2种不同的硬件平台上实现Hestenes-Jacobi算法。FPGA采用在Altera DE2-115开发板的Cyclone IV EP4CE115F29C7N上实现,工作频率为250 MHz。而CPU选择的是Intel(R) CoreTM i5-8500处理器,主频为3.0 GHz,利用MATLAB软件进行奇异值分解运算,与本文设计的电路结构性能进行对比。选择用于测试的矩阵维数为8×8,数据用32位定点数表示,且测试数据全部由MATLAB生成。为了验证模块的正确性,利用ModelSim软件作为仿真验证平台,对COM_Top运算单元中的主要模块进行仿真与分析。

3.1 向量乘积运算模块功能仿真

本次设定的矩阵A维数为8×8,因此,针对向量乘积运算模块的仿真验证,通过编写test_benth测试文件,将矩阵A中2个8×8的列向量中的元素依次送入向量乘积运算模块进行运算,对该运算单元的功能测试结果如图 6

图 6 向量乘积运算单元仿真 Fig.6 Vector_Mul unit simulation

图 6中,AB分别表示为2个列向量的元素,Vector_MulAA, Vector_MulAB和Vector_MulBB则为计算得到的列向量的内积,done表示为该组列向量对的内积计算完成,否则即为内积运算的中间值。输入数据和输出结果之间有一个时钟周期延时,且在该组列向量对计算完成后,直接开始下组列向量的内积运算。

3.2 CORDIC运算模块功能仿真

针对CORDIC运算单元的验证,本文选择矢量模式,实现求解反正切函数,对其进行仿真、验证。本次设置的流水线深度为12级,设置输入为正切函数值,输出即为求得的角度值。仿真结果如图 7

图 7 CORDIC运算单元仿真 Fig.7 CORDIC unit simulation

图 7中,y_in是输入的函数值信号,result即为输出信号:角度值,而result_valid表示的是输出信号有效的状态信号。result信号的第1个输出值与y_in的第1个输入数据之间间隔12个周期的延时,第2个输出结果与第1个输出之间则存在1个时钟延时。分析其原因,因为设置的流水线深度为12级,所以运行完整条流水线需要12个周期。同时,由于流水线展开的CORDIC算法,大大提高了模块的整体运算效率。

3.3 奇异值分析

完成基于FPGA的奇异值分解硬件设计及仿真后,对该设计进行奇异值分解得到的奇异值与MATLAB进行分解运算得到的奇异值进行对比,结果如表 2

表 2 奇异值分解计算结果 Tab.2 SVD calculation result

表 2中可看出,基于FPGA硬件结构分解得到的奇异值与在MATLAB中运算得到的奇异值并不完全相同,分析其原因是由于数据的精度不同造成的。由于本文使用的是32位的定点数,数据在浮点-定点转换以及在进行乘法运算时,由于定点截断导致运算存在误差,但误差相较小,可忽略不计。

3.4 系统运行时间分析

表 3列出了同一算法在FPGA平台以及CPU平台下奇异值分解算法的运行时间,尽管FPGA的工作频率远低于CPU的主频,FPGA的运算耗时却远小于CPU的耗时,大大提高了算法的运行速率。

表 3 SVD运行时间对比 Tab.3 SVD calculation time comparison  
4 结论

针对奇异值分解过程中计算复杂、实时性差的问题,提出了一种改进的基于CORDIC算法的实现方案。通过对平面旋转变换运算进行分析,提出利用全流水线的CORDIC运算单元实现,去除乘法运算过程,实现列向量的旋转变换。在此基础上进行实验验证,实验结果验证,该设计方案合理、有效,具有一定的工程应用价值。

参考文献
[1]
郭强.并行Jacobi方法求解矩阵奇异值的研究[D].苏州: 苏州大学, 2011.
GUO Q. Research on Parallel Jacobi Method for SVD Problem[D]. Shuzhou: Soochow University, 2011. http://d.wanfangdata.com.cn/thesis/Y1991195
[2]
张倩倩, 詹永照, 林涵阳. 基于分块和NSCT-SVD的彩色图像数字水印算法[J]. 江苏大学学报(自然科学版), 2017, 38(5): 582-589.
ZHANG Q Q, ZHAN Y Z, LIN H Y. Digital Watermarking algorithm for Color Images based on Blocks and NSCT-SVD[J]. Journal of Jiangsu University(Natural Science Edition), 2017, 38(5): 582-589.
[3]
余成波, 张林, 龙曦. 基于FPGA数字PLL谐振频率的跟踪研究[J]. 重庆理工大学学报(自然科学), 2019, 33(4): 141-146.
YU C B, ZHANG L, LONG X. Tracking Research of Digital PLL Resonant Frequency Based on FPGA[J]. Journal of Chongqing University of Technology (Natural Science), 2019, 33(4): 141-146.
[4]
邓翔宇, 刘增力. 基于改进的MCA和K-SVD的图像稀疏表示去噪算法[J]. 四川大学学报(自然科学版), 2016, 53(4): 774-780.
DENG X Y, LIU Z L. Image denoising via the sparse representation based on improved MCA and K-SVD[J]. Journal of Sichuan University (Natural Science Edition), 2016, 53(4): 774-780.
[5]
LAHABAR S, NARAYANAN P J. Singular value decomposition on GPU using CUDA[C]//2009 IEEE International Symposium on Parallel & Distributed Processing. Rome, Italy: IEEE, 2009: 1-10.
[6]
SNOPCE H, SPAHIU I. Parallel Computation of the SVD[C]//Proceedings of the 5th European conference on European computing conference. Paris, France: WSEAS, 2011: 456-461.
[7]
BRAVO I, MAZO M, LÁZARO J L, et al. Novel HW architecture based on FPGAs oriented to solve the eigen problem[J]. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 2008, 16(12): 1722-1725. DOI:10.1109/TVLSI.2008.2001939
[8]
RAHMATI M, SADRI M S, NAEINI M A. FPGA based singular value decomposition for image processing applications[C]//2008 International Conference on Application-Specific Systems, Architectures and Processors. Leuven, Belgium: IEEE, 2008: 185-190.
[9]
唐吉卓.基于GPU平台的SVD并行计算研究与实现[D].成都: 电子科技大学, 2014.
TANG J Z. Research and Implementation of Parallel Computing of SVD based on GPU[D]. Chengdu: University of Electronic Science and Technology of China, 2014.
[10]
BRENT R P, LUK F T, VAN L C. Computation of the singular value decomposition using mesh-connected processors[R]. Ithaca, NY, USA: Cornell University, 1982.
[11]
MA W, KAYE M E, LUKE D M, et al. An fpga-based singular value decomposition processor[C]//2006 Canadian Conference on Electrical and Computer Engineering. Ottawa, Ont., Canada: IEEE, 2006: 1047-1050.
[12]
马亚峰.基于FPGA的矩阵奇异值分解加速方案的设计与实现[D].北京: 北京交通大学, 2017.
MA Y F. Design and Implementation of Singular Value Decomposition Acceleration Scheme[D]. Beijing: Beijing Jiaotong University, 2017. http://cdmd.cnki.com.cn/Article/CDMD-10004-1017053745.htm
[13]
CHEN L, HAN J, LIU W, et al. Algorithm and design of a fully parallel approximate coordinate rotation digital computer (CORDIC)[J]. IEEE Transactions on Multi-Scale Computing Systems, 2017, 3(3): 139-151. DOI:10.1109/TMSCS.2017.2696003
[14]
GARRIDO M, KÄLLSTRÖM P, KUMM M, et al. CORDIC Ⅱ: a new improved CORDIC algorithm[J]. IEEE Transactions on Circuits and Systems Ⅱ: Express Briefs, 2016, 63(2): 186-190. DOI:10.1109/TCSII.2015.2483422