METHOD OF REALIZING ACCELERATED PARALLEL JACOBI COMPUTING FOR FPGA

Abstract

The invention discloses a method of realizing accelerated parallel Jacobi computing for an FPGA. Data of a n×n-dimensional matrix are input to the FPGA, and a rotation transformation process is carried out by using parallel Jacobi computing. Processors are initialized. A diagonal processor computes a symbol set corresponding to a rotation angle and outputs the symbol set to a non-diagonal processor. Elements of the diagonal processor are updated. Elements of the non-diagonal processor are updated. Elements between the processors are exchanged. After the elements of the respective processors are updated, the updated elements between the processors are exchanged. The invention requires less FPGA resources while yields a higher internal computational processing performance of the FPGA. Accordingly, the invention is capable of facilitating the efficiency of realizing eigenvalue decomposition in the FPGA and is highly applicable in actual processing.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention relates to an internal data processing method for an FPGA, and particularly relates to a method of realizing accelerated parallel Jacobi computing for an FPGA.

2. Description of Related Art

Many algorithms in various fields, such as radars, wireless communication, image processing, etc., need to compute the eigenvalues of a matrix. For example, the computation of eigenvalues is a key step in sub-space-based direction-of-arrival (DOA) estimation algorithms and principle component analysis (PCA) algorithms.

Currently, algorithms for computing a large number of eigenvalues include, for example, QR algorithms, LU decomposition algorithms, algebraic methods, etc. Since the complexity of root extraction in algebraic methods increases as the dimensions of the matrix increase, algebraic methods are not suitable for obtaining eigenvalues from a large-scale matrix. Meanwhile, LU decomposition algorithms are only suitable for obtaining eigenvalues from an invertible matrix. Besides, while QR algorithms have been proven as faster than serial Jacobi computing in the computation for eigenvalues, studies have nonetheless shown that Jacobi computing exhibits a higher accuracy than that of QR algorithms. Jacobi computing refers to a process of gradually converting a matrix into a nearly-diagonal matrix through a series of rotations. The diagonal elements in the matrix are the eigenvalues of the matrix. In addition, owing to the inherent parallelism of Jacobi computing resulting from its decomposition of eigenvalues from a real symmetric matrix, parallel Jacobi computing (a parallel method for realizing Jacobi computing) has been broadly applied for eigenvalue decomposition in FPGA-related applications.

Currently, some studies have been undertaken in hope to accelerate parallel Jacobi computing. However, most of these acceleration methods are unable to realize one step of parallel Jacobi computing in one CORDIC algorithm cycle. While the conventional approximate Jacobi computing is capable of realizing one step of parallel Jacobi computing in one CORDIC algorithm cycle, the performance is not entirely satisfactory. This is because the rotations are approximate rotations, which increase the total number of rotations required. Besides, while the total lookup table (LUT) resources in an FPGA are limited, the consumption of the LUT resources in the FPGA is not considered when the conventional algorithms are put into practice.

SUMMARY OF THE INVENTION

To address the above issues, the invention proposes a method of realizing accelerated parallel Jacobi computing for an FPGA and designs a solution which brings forth favorable computation processing efficacies when the parallel Jacobi computing method is realized in the FPGA. With the method, the technical issues that internal data processing in an FPGA is slow and resource-consuming are addressed, and the objective of realizing one step of parallel Jacobi computing within one CORDIC algorithm cycle is achieved. Besides, the resource consumption in the FPGA is reduced.

To achieve the objectives, the invention provides a technical solution with the following steps.

(1) Initializing Processors:

Data of a n×n-dimensional matrix are input into the FPGA and a rotation transformation process using the parallel Jacobi computing is carried out. A coordinate rotation digital computer (CORDIC) algorithm is adopted in Jacobi computing to carry out a planar rotation, and a two-dimensional X-Y coordinate system is established in the planar rotation.

A plurality of processors are provided in a field programmable gate array (FPGA). The processors are arranged in an array. Each of the processors is connected with an adjacent processor via a data interface to exchange data and elements, and each element in the n×n-dimensional matrix for carrying out the parallel Jacobi computing is assigned to a processor P_ij according to a formula as follows:

$P_{\mathrm{ij}}=\left(\begin{array}{cc}{a}_{2i-1,2j-1}& a_{2i-1,2j}\\ a_{2i,2j-1}& a_{2i,2j}\end{array}\right),i\le j,j=1,2,\dots ,\frac{n}{2},$

wherein P_ij represents a processor in an i^th row and a j^th column, a_2i,2j represents an element in a 2i^th row and a 2j^th column in the n×n-dimensional matrix, and Y represents dimensionality of the matrix.

A processor P_ij whose subscripted symbol satisfies i=j is a diagonal processor and a processor P_ij whose subscripted symbol does not satisfy i=j is a non-diagonal processor. In the processors P_ij, an element whose subscripted symbol satisfies 2i=2j and 2i−1=2j−1 is a diagonal element, and an element whose subscripted symbol does not satisfy 2i=2j and 2i−1=2j−1 is a non-diagonal element.

The n×n-dimensional matrix is a real symmetric matrix. Therefore, with the assignment of the above process, only the upper right portion is preserved, and the lower left portion and the upper right portion are symmetric to each other with respect to the diagonal.

(2) Computing a Symbol Set Corresponding to a Rotation Angle 2θ by the Diagonal Processor and Outputting the Symbol Set to the Non-Diagonal Processor:

A symbol set {d_2θ,k}, k=1, 2, . . . , N which corresponds to a rotation angle 2θ of the CORDIC algorithm, is obtained through iterations by using a formula as follows. A total number of the iterations is the same as a total number of iterations of the CORDIC algorithm:

$\tan \left({\theta }_k\right)=\frac{{\alpha }_k}{{\beta }_k}=\tan \left({\theta }_{k-1}-d_{2\theta ,k}{\varphi }_{k-1}\right)=\frac{\tan {\theta }_{k-1}-d_{\theta ,k}\tan {\varphi }_{k-1}}{1+d_{\theta ,k}\tan {\theta }_{k-1}\tan {\varphi }_{k-1}}=\frac{{\alpha }_{k-1}{2}^{k-1}-d_{\theta ,k}{\beta }_{k-1}}{{\beta }_{k-1}{2}^{k-1}+d_{\theta ,k}{\alpha }_{k-1}},d_{2\theta ,k}=\{\begin{array}{c}-1,\frac{{\alpha }_{k-1}}{{\beta }_{k-1}}<0\\ 1,\frac{{\alpha }_{k-1}}{{\beta }_{k-1}}\ge 0\end{array},{\theta }_0=2\theta ,\tan \left({\varphi }_{k-1}\right)=2^{-\left(k-1\right)},k=1,2,\dots ,N,$

wherein k represents an ordinal number of an iteration, N represents the total number of the iterations and is set as a data bit number adopted by the FPGA, α_k represents a first symbol parameter of a k^th iteration, β_k represents a second symbol parameter of the k^th iteration, θ₀ represents a rotation angle initial value, that is, 2θ, θ_k represents a residual rotation angle through k times of iterations, ϕ_k−1 represents an angle parameter of a (k−1)^th iteration, and d_2θ,k represents a symbol corresponding to the rotation angle 2θ at the k^th iteration.

Specifically, in a symbol computing module, d_2θ,k is obtained by performing an exclusive OR operation on the sign bits of α_k−1 and β_k−1, where d_2θ,k is 1 if the sign bits are the same, and d_2θ,k is −1 if the sign bits are opposite. α_k−12^k−1 is obtained from α_k−1 through a shift operation, and β_k−12^k−1 is obtained from β_k−1 through a shift operation. If d_2θ,k is 1, α_k is obtained by performing a subtract operation on α_k−12^k−1 and β_k−1, and β_k is obtained by performing an add operation on β_k−12^k−1 and α_k−1. If d_2θ,k is −1, α_k is obtained by performing an add operation on α_k−12^k−1 and β_k−1, and β_k is obtained by performing a subtract operation on β_k−12^k−1 and α_k−1.

Iterative computing starts, and an initial rotation angle corresponding to the non-diagonal elements in the diagonal processor is θ. The computation is as follows:

$\tan \left(2\theta \right)=\frac{{\alpha }_0}{{\beta }_0},$

a₀=2a_pq, β₀=a_pp−a_qq,wherein a_pq, a_qp respectively represent two non-diagonal elements initially included in the diagonal processor, a_qp=a_pq, a_pp, and a_pp respectively represent diagonal elements initially included in the diagonal processor, α₀ represents an initial first symbol parameter, and β₀ represents an initial second symbol parameter.

The diagonal elements a_pp and a_qq of the diagonal processor obtain β₀=a_pp−a_qq through a subtract operation. The non-diagonal element a_pq obtains α₀=2a_pq in a shift operation. It is set that, in the parallel Jacobi computing, the rotation angle corresponding to the non-diagonal element in the current diagonal processor is θ, and β₀ and α₀ are input, as initial values, into a symbol set computing module, and the symbol computing module obtains a symbol set {d_2θ,k} corresponding to the rotation angle 2θ through iterations.

The diagonal processor outputs the rotation angle 2θ obtained through computing carried out by itself and and the corresponding symbol set {d_2θ,k} to a non-diagonal processor on the same row and a non-diagonal processor on the same column.

(3) Updating Elements of the Diagonal Processor:

The CORDIC algorithm is carried out on first to-be-rotated coordinates (2a_pq,a_pp−a_qq) by using d_2θ,k obtained in each of the iterations in Step (2) as a rotation symbol of the k^th iteration in the CORDIC algorithm. This process replaces the step of calculating a rotation symbol after each iteration in the conventional CORDIC algorithm. Accordingly, a planar rotation is carried out by using the rotation angle 2θ.

After all the iterations in Step (2) are completed, a final planar rotation result is multiplied by a first compensation factor to obtain rotated y coordinates, that is, y₁=2a_pq sin 2θ+(a_pp−a_qq)cos 2θ. The first compensation factor is obtained according to a formula as follows:

$C_1=\prod _{k=1}^N\cos \left({\varphi }_{k-1}\right),$

wherein C₁ represents the first compensation factor.

Diagonal elements in the diagonal processor are updated by using a formula as follows, and non-diagonal elements are set to 0:

$a_{\mathrm{pp}}^{\prime }=\frac{a_{\mathrm{qq}}+a_{\mathrm{pp}}+y_1}{2},a_{\mathrm{qq}}^{\prime }=\frac{a_{\mathrm{qq}}+a_{\mathrm{pp}}-y_1}{2},$

wherein a′_pp, a′_qq represent two updated diagonal elements in the diagonal processor, y₁ represents a rotated Y-axis coordinate of the first to-be-rotated coordinates.

(4) Updating Elements of the Non-Diagonal Processor:

(4.1) The non-diagonal processor P_ij receives symbol sets output from two diagonal processors P_ii, P_jj and represented as {d_2θi,k}, {d_2θj,k}. d_2θi,k and d_2θj,k respectively represent symbols corresponding to a rotation angle 2θ₁ and a rotation angle 2θ_j at the k^th iteration, and two symbols d_θi+θj,k and d_θi−θj,k are respectively computed by using formulae as follows to obtain two symbol sets {d_θi+θj,k} and {d_θi−θj,k}:

d _{θ i +θ j ,k}=½(d_2θi,k+d_2θj,k), k=1,2, . . . ,N,

d _{θ i −θ j ,k}=½(d_2θi,k−d_2θj,k), k=1,2, . . . ,N,

wherein d_θi+θj,k and d_θi−θj,k respectively represent symbols corresponding to a rotation angle θ_i+θ_j and a rotation angle θ_i−θ_j, and 2θ₁ and 2θ_j respectively represent double angles of rotation angles corresponding to the non-diagonal elements of the two diagonal processors P_ii and P_jj.

Specifically, by subjecting each pair of symbols d_2θi,k, d_2θj,k to an exclusive OR and a data selector, the symbol set for the rotation angle θ_l-θ_m is determined. If the result of the exclusive OR operation is 1, d_2θl,k is taken as d_θl−θm,k, otherwise 0 is taken as d_θl−θm,k. By subjecting each pair of symbols d_2θl,k, d_2θj,k to an exclusive NOR and a data selector, the symbol set for the rotation angle θ_l+θ_m is determined. If the result of the exclusive NOR operation is 1, d_2θl,k is taken as d_θl+θm,k, otherwise 0 is taken as d_θl+θm,k.

(4.2) d_θi±θj,k has three values, i.e., {−1,0,1}. Values of a second compensation factor and a third compensation factor corresponding to all possible symbol combinations formed by first

$\lceil \frac{N}{2}\rceil$

symbols in the two symbol sets {d_θi+θj,k} and {d_θi−θj,k} are computed by using formulae as follows, one symbol combination being formed by

$\lceil \frac{N}{2}\rceil$

symbols, so as to establish lookup table data by using the values of the second compensation factor and the third compensation factor corresponding to the respective, different symbol combinations. The absolute value of each symbol in the first

$\lceil \frac{N}{2}\rceil$

symbols serves as a lookup address, and a lookup table is generated by using a block memory (block random access memory). An address bit number of the lookup table is set as

$\lceil \frac{N}{2}\rceil ,$

and a data depth is

$2^{\lceil \frac{N}{2}\rceil };$

$C_2=\stackrel{\stackrel{k=1}{\lceil \frac{N}{2}\rceil }}{\Pi }\cos \left(d_{{\theta }_i-{\theta }_j,k}{\varphi }_{k-1}\right),d_{{\theta }_i-{\theta }_j,k}\in \left\{-1,0,1\right\},C_3=\stackrel{\stackrel{k=1}{\lceil \frac{N}{2}\rceil }}{\Pi }\cos \left(d_{{\theta }_i+{\theta }_j,k}{\varphi }_{k-1}\right),d_{{\theta }_i+{\theta }_j,k}\in \left\{-1,0,1\right\},$

wherein C₂ represents the second compensation factor, and C₃ represents the third compensation factor.

When the number of iterations of the CRODIC algorithm exceeds

$\lceil \frac{N}{2}\rceil$

(rounding up

$\frac{N}{2}),$

the difference between the second and third compensation factors and 1 is less than 2^−N+1, and the accuracy of N-bit signed fixed-point number is at most 2^−N+1. Therefore, the remaining second and third compensation factors may be directly considered as 1, i.e., no compensation is required.

(4.3) For the non-diagonal processor, four elements included in the non-diagonal processor are represented as

$\left(\begin{array}{cc}{a}_{p_1{q}_1}& a_{p_1{q}_2}\\ a_{p_2{q}_1}& a_{p_2{q}_2}\end{array}\right).$

The obtained d_θi−θj,k is used as the rotation symbol of the k^th iteration in the CORDIC algorithm, and the CORDIC algorithm is carried out on second to-be-rotated coordinates (a_p1q1+a_p2q2,a_p1q2−a_p2q1) to carry outa planar rotation by using the rotation angle θ_i-θ_j. A planar rotation result is multiplied by the second compensation factor whose value is obtained by accessing the lookup table of Step (4.2) to obtain rotated coordinates represented as:

$\{\begin{array}{c}{x}_2=\left(a_{p_1{q}_1}+a_{p_2{q}_2}\right)\cos \left({\theta }_i-{\theta }_j\right)-\left(a_{p_1{q}_2}-a_{p_2{q}_1}\right)\sin \left({\theta }_i-{\theta }_j\right)\\ y_2=\left(a_{p_1{q}_2}-a_{p_2{q}_1}\right)\cos \left({\theta }_i-{\theta }_j\right)+\left(a_{p_1{q}_1}+a_{p_2{q}_1}\right)\sin \left({\theta }_i-{\theta }_j\right)\end{array}\hspace{1em}$

wherein x₂ and y₂ respectively represent rotated coordinates of the second to-be-rotated coordinates.

The obtained d_θi+θj,k is used as the rotation symbol of the k^th iteration in the CORDIC algorithm. The CORDIC algorithm is carried out on third to-be-rotated coordinates (a_p1q2+a_p2q1,a_p1q1−_p2q2) to carry out a planar rotation by using the rotation angle θ_i+θ_j. A planar rotation result is multiplied by the third compensation factor whose value is obtained by accessing the lookup table of Step (4.2) to obtain rotated coordinates represented as:

$\{\begin{array}{c}{x}_3=\left(a_{p_1{q}_2}+a_{p_2{q}_1}\right)\cos \left({\theta }_i+{\theta }_j\right)-\left(a_{p_1{q}_1}-a_{p_2{q}_2}\right)\sin \left({\theta }_i+{\theta }_j\right)\\ y_3=\left(a_{p_1{q}_1}-a_{p_2{q}_2}\right)\cos \left({\theta }_i+{\theta }_j\right)+\left(a_{p_1{q}_2}+a_{p_2{q}_1}\right)\sin \left({\theta }_i+{\theta }_j\right)\end{array}\hspace{1em}$

wherein x₃ and y₃ respectively represent rotated coordinates of the third to-be-rotated coordinates.

(4.4) Adopting formulae as follows to update the elements in the non-diagonal processor:

a′ _{p 1 q 1} =½(x₂+y₃),

a′ _{p 1 q 2} =½(x₃+y₂),

a′ _{p 2 q 1} =½(x₃−y₂),

a′ _{p 2 q 2} =½(x₂−y₃),

wherein a′_p1q1, a′_p1q2, a′_p2q1, and a′_p2q2 respectively represent four elements included in the non-diagonal processor.

Specifically, by performing an add operation and a shift operation on x₂ and y₃, an updated value a′_p1q1 as shown in the formula is obtained. By performing an add operation and a shift operation on x₃ and y₂, an updated value a′_p1q2 as shown in the formula is obtained. By performing a subtract operation and a shift operation on x₃ and y₂, an updated value a′_p2q1 as shown in the formula is obtained. By performing a subtract operation and a shift operation on x₂ and y₃, an updated value a′_p2q2 as shown in the formula is obtained.

(5) Exchanging the Elements Between the Processors:

After the elements of the respective processors are updated, the matrix elements symmetric thereto are also updated to the same values. The updated elements between the processors are exchanged.

(5.A) The diagonal elements in the diagonal processor are exchanged.

It is assumed that a current diagonal processor P_ii includes a diagonal element a_pipi and a diagonal element a_qiqi.

For the diagonal element i represents a diagonal processor row/column ordinal number, if i=1, the diagonal element a_pipi is not changed, if i=2, a value of the diagonal element a_pipi is changed to a value of a diagonal element a_qi+1qi+1, and if i>2, the value of the diagonal element a_pipi is changed to a value a_pi+1pi+1 of a diagonal element.

For the diagonal element a_qiqi, if

$i<\frac{n}{2},$

a value of the diagonal element a_qiqi is changed to a value of a_qi+1qi+1, and if

$i=\frac{n}{2},$

the value of the diagonal element a_qiqi is changed to the value of the diagonal element a_pipi.

(5.B) The non-diagonal elements in the diagonal processor and elements in the non-diagonal processor are exchanged by changing positions according to the following. Positions of the non-diagonal elements in the diagonal processor and the elements in the non-diagonal processor are shifted. It is noted that the elements may be shifted out of a processor to another processor. Accordingly, a subscripted row symbol of an element is the same as a row number of a diagonal element shifted to the same row after the exchanging of Step (5.A), and a subscripted column symbol of the element is the same as a column number of a diagonal element shifted to the same column after the exchanging of Step (5.A).

(6) Updating the non-diagonal elements in all the diagonal processors in the n×n-dimensional matrix by the Jacobi computing after the exchanging, returning to Step (2) for another round of processing and updating, repeating the updating until the non-diagonal elements in the n×n-dimensional matrix gradually converge to 0, finishing the updating when a predetermined convergence accuracy is met, and ending the parallel Jacobi computing.

The n×n-dimensional matrix is a covariance matrix of data collected by an antenna array or data before image dimensionality reduction, and is a real symmetric matrix.

In Step (1), if n in the n×n-dimensional matrix is an odd number, i.e., a matrix with odd-numbered dimensionality, the matrix is expanded into a matrix with even-number dimensionality by adding a n+1^th column and a n+1^th row, and element values of the added n+1^th column and n+1^th row are all set to 0.

In the invention, the k^th symbol computed in Step (2) is provided to the CORDIC algorithms of Steps (3) and (4) for the k^th iteration. Therefore, Steps (2), (3), and (4) are carried out simultaneously.

The updated values of the elements in the respective processors can be obtained by simply combining the computing results of the CORDIC algorithms. Therefore, operations for the elements of all the processors can be carried out simultaneously. According to the method of the invention, the time consumed for realizing one step of the parallel Jacobi computing is only one CORDIC cycle. Compared with the conventional method which consumes three CORDIC cycles, the computing time is significantly reduced, and the computing performance is facilitated.

The data of the n×n-dimensional matrix of the invention may be data collected by an antenna array for DOA estimation, or a covariance matrix used for reducing the dimensionality of image data by using a PCA algorithm.

The beneficial effects of the invention include the following.

The invention adopts a specifically designed linear combination method to replace the bilateral rotation method in the conventional parallel Jacobi computing. By using the symbol set of the rotation angle together with the combination of two symbol sets to replace the step of computing the rotation symbols in the CORDIC algorithm, the parallelism of the parallel Jacobi computing is facilitated, the computing time for each step in the parallel Jacobi computing is reduced, and one step of the parallel Jacobi computing can be realized in one CORDIC cycle.

The invention effectively facilitates the speed of realizing the parallel Jacobi computing in hardware. By realizing one step of the parallel Jacobi in one CORDIC algorithm cycle, the deficiency of the conventional method is alleviated. In addition, the time required is only one-third of that required by the conventional method.

The invention requires less FPGA resources while yields a higher internal computational processing performance of the FPGA. Accordingly, the invention is capable of facilitating the efficiency of realizing eigenvalue decomposition in the FPGA and is highly applicable in actual processing.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram illustrating a framework of a diagonal processor according to an embodiment of the invention.

FIG. 2 is a schematic diagram illustrating a framework of a non-diagonal processor according to an embodiment of the invention.

FIG. 3 is a schematic diagram illustrating a framework of a processor array according to an embodiment of the invention.

FIG. 4 is a flowchart illustrating a computing method according to an embodiment of the invention.

DESCRIPTION OF THE EMBODIMENTS

In the following, details of the invention will be described with reference to the accompanying drawings in combination with exemplary embodiments of the invention.

The framework realized in an FPGA of the invention mainly includes a diagonal processor and a non-diagonal processor. The framework of the diagonal processor is as shown in FIG. 1, and the framework of the non-diagonal processor is as shown in FIG. 2. A framework of a processor array is as shown in FIG. 3. A flowchart for executing a computing method is as shown in FIG. 4.

The embodiment of the invention and the implementation process thereof are described in the following.

The specific implementation processes of the embodiment are realized in a Xilinx Virtex-7 XC7VX690T FPGA chip. Specifically, wireless signals emitted by a collector drone with a four-element antenna array is adopted in the implementation, and the signal incident direction is 0 degrees. A 4×4 real symmetric covariance matrix obtained through the computation carried out according to four data sets received by the four-element antenna is represented as A.

16-bit fixed-point number is adopted to obtain the eigenvalues under a condition that

$A=\left(\begin{array}{cccc}237.4904& 231.6229& -104.5409& 24.3696\\ 231.6229& 541.7360& -78.0729& -10.1869\\ -104.5409& -78.0729& 273.6585& -34.0290\\ 24.3696& -10.1869& -34.0290& 170.5949\end{array}\right).$

Specifically, the following steps are included.

(1) Initializing the processors: The respective elements in R_r are assigned to the processors P_ij. Each processor is connected with the adjacent processor via a data interface. A processor whose subscripted symbol satisfies i=j is defined as a diagonal processor, and a processor whose subscripted symbol does not satisfy such condition is defined as a non-diagonal processor. A matrix element whose subscripted symbol satisfies 2i=2j and 2i−1=2j−1 is defined as a diagonal element, and a matrix element whose subscripted symbol does not satisfy such condition is defined as a non-diagonal element.

(2) Computing a symbol set corresponding a rotation angle by the diagonal processor and outputting the symbol set to the non-diagonal processor: It is assumed that the non-diagonal elements included in the diagonal processor are a_pq, a_qp, and a_qp=a_pq. It is assumed that the diagonal elements included in the diagonal processor are a_pp and a_pp. It is also assumed that α₀=2a_pq, β₀=a_pp−a_qq. It is assumed that the rotation angle corresponding to the non-diagonal element in the current diagonal processor is θ. A symbol set d_2θ,k, k=1, 2 . . . , 16, which corresponds to a rotation angle 2θ of the CORDIC algorithm is obtained through iterations. The number of iterations is the same as the number of iterations in the CORDIC algorithm, and the data bit number, i.e., 16, adopted in the current system is used.

(3) Updating the elements of the diagonal processor: A compensation factor is obtained. d_2θ,k obtained in Step (2) is used as the rotation symbol of the k^th iteration in the CORDIC algorithm. This process replaces the step of computing the rotation symbol after each iteration in the conventional CORDIC algorithm. The CORDIC algorithm is executed to rotate (2a_pq,a_pp−a_qq) by 2θ, and the result is multiplied by the compensation factor to obtain rotated y coordinates, i.e., y₁=2a_pq sin 2θ+(a_pp−a_qq) cos 2θ, and the diagonal elements in the diagonal processor are updated. In addition, the non-diagonal elements are set to 0.

(4) Updating the elements of the non-diagonal processor: The non-diagonal processor P_ij receives the symbol sets output from the two diagonal processors P_ii, P_jj, and the symbol sets are represented as d_2θi,k, d_2θj,k k=1, 2, . . . , 32. d_θi+θj,k and d_θi−θj,k are respectively computed.

d_θi±θj,k has three values, i.e., {−1,0,1}. The values of the compensation factors corresponding to all the possible value combinations for the first 16 symbols of the symbol set d_θi±θj,k are computed. The values of the compensation factors are used as lookup table data. The absolute values of the respective symbols of the first 16 symbols in the symbol set are used as lookup addresses, and a lookup table is generated by using a block memory. When the number of iterations of the CRODIC algorithm exceeds 8, the difference between the compensation factor and 1 is less than 2⁻⁷, and the accuracy of 8-bit data is at most 2⁻⁷. Therefore, the remaining compensation factors may be directly considered as 1, i.e., no compensation is required. 8 is set as the address bit number of the lookup table, and the data depth is 2⁸. The lookup table of the example is as shown in Table 1.

TABLE 1
Lookup Table for Compensation Value
Address   Compensation value
0000 0000 1
0000 0001 0.99996948
.         .
.         .
.         .
1111 1111 0.60726543

It is assumed that the matrix elements included in the current non-diagonal processor are

$\left(\begin{array}{cc}{a}_{p_1{q}_1}& a_{p_1{q}_2}\\ a_{p_2{q}_1}& a_{p_2{q}_2}\end{array}\right),$

and d_θl−θm,k is used as the rotation symbol of the k^th iteration in the CORDIC algorithm. A CORDIC algorithm rotation θ_l-θ_m is carried out on (a_p1q1+a_p2q2,a_p1q2−a_p2q1). A compensation factor is obtained from the lookup table. The result is multiplied by the compensation factor, thereby deriving the rotated coordinates.

d_θl−θm,k is used as the rotation symbol of the k^th iteration in the CORDIC algorithm. A CORDIC algorithm rotation θ_l+θ_m is carried out on (a_p1q2+a_p2q1, a_p1q1−a_p2q2). A compensation factor is obtained from the lookup table. The result is multiplied by the compensation factor, thereby deriving the rotated coordinates.

The elements of the non-diagonal processor are updated.

(5) Exchanging the elements between the processors: After the elements of the respective processors are updated, the matrix elements symmetric thereto are also updated to the same values. The updated elements are exchanged with the elements of other processors.

Then, the flow returns to Steps 2, 3, and 4 again for another round of computation and update. After three times of exchange, all the non-diagonal elements in the matrix have been updated once by the diagonal processors through Jacobi computing. Through multiple times of update, the non-diagonal elements in the matrix gradually converge to 0. The update ends after a predetermined convergence accuracy set by the user is met, and the parallel Jacobi computing ends.

The specific results are as follows:

First Round:

$P_{11}=\left(\begin{array}{cc}237.4904& 231.6229\\ 231.6229& 541.7360\end{array}\right),P_{22}=\left(\begin{array}{cc}273.6585& -34.0290\\ -34.0290& 170.5949\end{array}\right),P_{12}=\left(\begin{array}{cc}-104.5409& 16.3696\\ -78.0729& -10.1869\end{array}\right),$

and the following is rendered after the update:

${\stackrel{\_}{A}}_1=\left(\begin{array}{cccc}112.5& 0& -58.0625& 2.5625\\ 0& 666.75& -112.9375& -35.1875\\ -58.0625& -112.9375& 283.875& 0\\ 2.5625& -35.1875& 0& 160.375\end{array}\right),$

and the following is rendered after the element exchange:

$A_1=\left(\begin{array}{cccc}112.5& 2.5625& 0& -58.0625\\ 2.5625& 160.375& -35.1875& 0\\ 0& -35.1875& 666.75& -112.9375\\ -58.0625& 0& -112.9375& 283.875\end{array}\right).$

Second Round:

$P_{11}=\left(\begin{array}{cc}112.5& 2.5625\\ 2.5625& 160.375\end{array}\right),P_{22}=\left(\begin{array}{cc}666.75& -112.9375\\ -112.9375& 283.875\end{array}\right),P_{12}=\left(\begin{array}{cc}0& -58.0625\\ -35.1875& 0\end{array}\right),$

and the following is rendered after the update:

${\stackrel{\_}{A}}_2=\left(\begin{array}{cccc}112.375& 0& 17.125& -55.4375\\ 0& 160.5& -33.0625& -12.25\\ 17.125& -33.0625& 697.625& 0\\ -55.4375& -12.25& 0& 253\end{array}\right),$

and the following is rendered after the element exchange:

$A_2=\left(\begin{array}{cccc}-112.375& -55.4375& 0& 17.125\\ -55.4375& 253& -12.25& 0\\ 0& -12.25& 697.625& -33.0625\\ 17.125& 0& -33.0625& 160.5\end{array}\right).$

Eighth Round:

$P_{11}=\left(\begin{array}{cc}92.5& -0.0625\\ -0.0625& 700.1875\end{array}\right),P_{22}=\left(\begin{array}{cc}273.5& -0.0625\\ -0.0625& 157.3125\end{array}\right),P_{12}=\left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right),$

and the following is rendered after the update:

${\stackrel{\_}{A}}_8=\left(\begin{array}{cccc}92.5& 0& 0& 0\\ 0& 700.25& 0& 0\\ 0& 0& 273.4375& 0\\ 0& 0& 0& 157.3125\end{array}\right)$

As shown above, the non-diagonal elements of the matrix has met the convergence condition (while parallel Jacobi computing is an algorithm that approximates the diagonal elements to 0, a fixed-point number with a limited number of bits is used in actual implementation to represent decimals, so even though the values of the non-diagonal elements may reach 0, errors are also introduced). At this time, the elements on the diagonal of Ā₈ are the eigenvalues that are obtained. By applying the obtained eigenvalues to a signal direction-of-arrival (DOA) estimation algorithm, as shown in the following diagram, the power spectrum function of a multiple signal classification (MUSIC) algorithm reaches the peak value at 0 degrees, which shows that the invention realizes a proper function.

In the embodiment, the performance of the invention in the actual application is demonstrated in two aspects, i.e., operation time, FPGA resource consumption.

Operation time: Since 16-bit fixed-point number is set for the data, the number of times of internal iterations is 16 in the CORDIC algorithm. Considering the result compensation, there are 17 FPGA clock cycles for the CORDIC algorithm cycle. Considering also that the elements need to be exchanged between steps of the parallel Jacobi, which takes one clock cycle, it takes a total of 18 clock cycles for the method of realizing the accelerated parallel Jacobi computing of the invention to realize one step of the parallel Jacobi. In the embodiment, the convergence condition is set as the absolute value of the maximum value of the non-diagonal elements in the covariance matrix being less than 0.001. The convergence condition is met after 8 iterations, which takes 144 clock cycles. The clock frequency used in the example is set as 250M, which takes 0.576 milliseconds.

Resource consumption: A Verilog program realizing the example is integrated on the Vivado 2017.1 software platform. The results shows that the example consumes 2360 lookup tables (LUTs) and 688 registers (REGs), which respectively takes up 0.54% and 0.79% of the total resources. According to the above, the design consumes only a limited amount of FPGA resources.

According to the conventional method for realizing parallel Jacobi, the diagonal processors require a CORDIC algorithm cycle to obtain rotation angles, and then require two successive CORDIC algorithms using the rotation angles obtained by the diagonal processors to update the elements of the diagonal processors. In other words, the conventional method requires a total of three CORDIC algorithm cycles. The non-diagonal processor needs to wait for the diagonal processor to obtain the rotation angle. Then, the non-diagonal processor also requires two consecutive CORDIC algorithms to update the elements of the non-diagonal processor. The angles of two rotations are respectively the rotation angle transmitted from the diagonal processor on the same row and the rotation angle transmitted from the diagonal processor on the same column. The respective processors operate in parallel, and it requires at least three CORDIC algorithm cycles to realize one step of parallel Jacobi. According to the CORDIC algorithm used in the invention, the operations are carried out in parallel and only require only one CORDIC algorithm cycle. The comparison between the invention and the conventional method in terms of the processes of the processors is as shown in Table 2.

TABLE 2
Comparison on Processes of Processors between the Invention and the
Conventional Parallel Jacobi Method
	Conventional method	Present invention
	Diagonal	Non-diagonal	Diagonal	Non-diagonal
Time	processor	processor	processor	processor
First	Compute a	Wait	Compute a	Compute a symbol
CORDIC	rotation		symbol set	set corresponding
algorithm	angle		corresponding	to an angle sum
cycle			to a double	and an angle
			angle of the	difference of two
			rotation angle	rotation angles and
			and carry out	carry out a rotation
			a rotation
Second	First	First rotation
CORDIC	rotation
algorithm
cycle
Third	Second	Second
CORDIC	rotation	rotation
algorithm
cycle

Based on the above, the invention significantly facilitates the speed of obtaining eigenvalues over the conventional method and is applicable when eigenvalue decomposition needs to be carried out quickly in actual processing.

The equivalent structural changes made by those skilled in the art based on the contents of the description and drawings of the invention shall be comprehensively covered in the scope of the invention.

Claims

1. A method of realizing accelerated parallel Jacobi computing for an FPGA, comprising: step (1) initializing processors:inputting data of a n×n-dimensional matrix into the FPGA and carrying out a rotation transformation process using the parallel Jacobi computing, wherein a CORDIC algorithm is adopted in the parallel Jacobi computing to carry out a planar rotation, and a two-dimensional X-Y coordinate system is established in the planar rotation,a plurality of processors are provided in the FPGA, the processors are arranged in an array, each of the processors is connected with an adjacent processor via a data interface to exchange data and elements, and each element in the n×n-dimensional matrix for carrying out the parallel Jacobi computing is assigned to a processor if of the processors according to a formula as follows:

2. The method of realizing the accelerated parallel Jacobi computing for the FPGA as claimed in claim 1, wherein in the step (2), an initial rotation angle corresponding to the non-diagonal elements in the diagonal processor when iterative computing starts is θ, and a computation is as follows:

3. The method of realizing the accelerated parallel Jacobi computing for the FPGA as claimed in claim 1, wherein the n×n-dimensional matrix is a covariance matrix of data collected by an antenna array or data before image dimensionality reduction, and is a real symmetric matrix.

4. The method of realizing the accelerated parallel Jacobi computing for the FPGA as claimed in claim 1, wherein in the step (1), if n in the n×n-dimensional matrix is an odd number, the n×n-dimensional matrix is expanded into a matrix with even-numbered dimensionality by adding a n+1th column and a n+1th row, and element values of the added n+1th column and n+1th row n+1th are all set to 0.

5. The method of realizing the accelerated parallel Jacobi computing for the FPGA as claimed in claim 1, wherein the step (4) comprises: step (4.1) receiving, by the non-diagonal processor Pij, symbol sets output from two diagonal processors Pii, Pjj and represented as {d2θi,k}, {d2θj,k}, wherein d2θi,k and d2θj,k respectively represent symbols corresponding to a rotation angle 2θi and a rotation angle 2θj at the kth iteration, and two symbols dθi+θj,k and dθi−θj,k are respectively computed by using formulae as follows to obtain two symbol sets {dθi+θj,k} and {dθi−θj,k}: dθi+θj,k=½(d2θi,k+d2θj,k), k=1,2, . . . ,N, dθi−θj,k=½(d2θi,k−d2θj,k), k=1,2, . . . ,N, wherein dθi+θj,k and dθi−θj,k respectively represent symbols corresponding to a rotation angle θi+θj and a rotation angle θi−θj, and 2θi and 2θj respectively represent double angles of rotation angles corresponding to the non-diagonal elements of the two diagonal processors Pii and Pjj;step (4.2) computing values of a second compensation factor and a third compensation factor corresponding to all possible symbol combinations formed by first

6. The method of realizing the accelerated parallel Jacobi computing for the FPGA as claimed in claim 1, wherein in the step (5), exchanging updated elements between the processors after the elements of each processor are updated, and the step (5) further comprises: step (5.A) exchanging the diagonal elements in the diagonal processor, wherein it is assumed that a current diagonal processor Pii comprises a diagonal element apipi and a diagonal element aqiqi, and, for the diagonal element apipi, i represents a diagonal processor row/column ordinal number, if i=1, the diagonal element apipi is not changed, if i=2, a value of the diagonal element apipi is changed to a value of a diagonal element aqi−1qi−1, and if i>2, the value of the diagonal element apipi is changed to a value api−1pi−1 of a diagonal element, and for the diagonal element aqiqi, if

7. The method of realizing the accelerated parallel Jacobi computing for the FPGA as claimed in claim 1, wherein the steps (2), (3), and (4) are carried out simultaneously.