QR decomposition matrix inversion optimization method and system

Abstract

The application discloses a QR decomposition matrix inversion optimization method and system. The method comprises the following steps: receiving signal data, performing FFT preprocessing and data downsampling; adopting double-buffering optimization parallel processing, storing the FFT signal data based on a double-buffering FIFO data storage structure; calculating the covariance matrix of the FFT signal data; adopting Givens rotation to gradually eliminate the non-diagonal elements of the covariance matrix, performing Givens rotation calculation, and adopting an optimized CORDIC calculation method to calculate the Givens rotation factor; adopting a systolic array architecture to accelerate the Givens rotation calculation process; and outputting the final QR decomposition result according to the Givens rotation calculation result. The method optimizes data flow, improves calculation throughput, reduces hardware resource occupation, provides a low-power and high-performance hardware implementation scheme for efficient matrix calculation, and has wide engineering application value.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of communication signal processing, and particularly relates to a QR decomposition matrix inversion optimization method and system. BACKGROUND

[0002] In GNSS anti-jamming, MIMO communication, radar signal processing and other applications, QR decomposition is a key algorithm for realizing beam forming, interference suppression and matrix preprocessing. The existing FPGA / ASIC scheme mainly adopts the following two methods for QR decomposition:

[0003] One is QR decomposition based on direct matrix inversion, which adopts Gram-Schmidt orthogonalization and calculates the matrix through column-by-column projection orthogonalization. This method calculates matrix inversion by solving linear equations, has high complexity, causes heavy hardware calculation burden, and is difficult to meet the real-time application requirements. At the same time, the Gram-Schmidt method is easy to introduce numerical error, causing numerical instability problem and affecting signal processing accuracy.

[0004] The other is QR decomposition based on Givens rotation + traditional division calculation, which first converts the matrix into an upper triangular matrix by Givens rotation, and then uses hardware multipliers and dividers to calculate the rotation factor. This method has large calculation amount including square root operation and division operation, is difficult to optimize, needs multiple multiplication units (MAC), has high FPGA / ASIC resource occupation, has high power consumption, and is difficult to meet the low-power chip requirements. At the same time, the Givens rotation factor needs to be calculated sequentially, which is difficult to realize parallel optimization, resulting in low throughput.

[0005] Therefore, the existing QR decomposition methods have defects such as high calculation complexity, large hardware resource consumption and low throughput. SUMMARY

[0006] In view of the deficiencies of the prior art, the purpose of the present application is to provide a QR decomposition matrix inversion optimization method, which optimizes the calculation complexity, data throughput and hardware resource utilization of QR decomposition, and is suitable for high-performance matrix calculation scenarios such as GNSS anti-jamming, MIMO communication and radar signal processing.

[0007] The second purpose of the present application is to provide a system for realizing the QR decomposition matrix inversion optimization method.

[0008] The present application provides a QR decomposition matrix inversion optimization method, comprising the following steps:

[0009] S1. receiving signal data, performing FFT preprocessing and data downsampling to obtain FFT signal data;

[0010] S2. Double-buffered optimized parallel processing is adopted, based on a double-buffered FIFO data storage structure, to store FFT signal data;

[0011] S3. Calculate the covariance matrix of the FFT signal data and use it as the Givens rotation input data;

[0012] S4. The Givens rotation is performed by gradually eliminating the off-diagonal elements of the covariance matrix, and the Givens rotation factor is calculated using an optimized CORDIC calculation method. The Givens rotation calculation process is accelerated by using a Systolic array architecture.

[0013] S5. Based on the Givens rotation calculation results, output the final QR decomposition results.

[0014] Step S1 includes the following steps:

[0015] Receive signal data, perform FFT calculations to obtain frequency point data;

[0016] Principal component analysis was performed on the frequency data to obtain the effective frequency data;

[0017] The effective frequency data is downsampled to obtain FFT signal data.

[0018] The principal component analysis specifically refers to:

[0019] Calculate the signal power spectrum of frequency point data Use the following formula to calculate: ;in, The complex spectrum value of the frequency point data obtained by FFT calculation; The complex number is the conjugate of the complex spectral value of the frequency point data;

[0020] Set the energy threshold again Select the frequency point where the cumulative energy percentage reaches the energy threshold, and express it using the following formula: ;in, This represents the final number of effective frequency points selected.

[0021] The downsampling ratio of the data is calculated using the following formula: .

[0022] When the source of the signal data to be processed is a GNSS signal or a Doppler frequency offset signal, the signal energy is unevenly distributed across different frequency bands. A non-uniform downsampling optimization scheme is adopted to replace the principal component analysis + downsampling scheme. The specific steps of the replaced step S1 are as follows:

[0023] Receive signal data, perform FFT calculations to obtain frequency point data;

[0024] Frequency-specific downsampling of frequency data is performed, and the non-uniform downsampling function within the frequency range is expressed by the following formula: ;in, The preset principal component threshold of the signal;

[0025] Set low-frequency and high-frequency ranges. In the low-frequency range, use the first frequency point retention ratio to retain frequencies, and in the high-frequency range, use the second frequency point retention ratio to retain frequencies; the first frequency point retention ratio is greater than the frequency point retention ratio.

[0026] Finally, the data downsampling ratio is calculated and the data is downsampled.

[0027] Specifically, in practical applications, the low-frequency region is 0~500MHz, with 80% of the frequency points retained; the high-frequency region is greater than 500MHz, with 20% of the frequency points retained.

[0028] In step S2, the double-buffered FIFO data storage structure includes buffer A and buffer B;

[0029] The dual-buffered optimized parallel processing is as follows: First, buffer A is used to store the currently obtained FFT signal data; when buffer A is insufficient, the process switches to buffer B to store the currently obtained FFT signal data, while simultaneously starting the calculation of step S3 on the FFT signal data already stored in buffer A; when buffer B is insufficient and the subsequent calculation of the FFT signal data already stored in buffer A is completed, buffer A is cleared, and the process switches back to buffer A to store the currently obtained FFT signal data, while simultaneously starting the calculation of step S3 on the FFT signal data already stored in buffer B; the two buffers are alternately switched for storage and reading calculation according to the above steps to achieve parallel processing.

[0030] In step S3, the covariance matrix of the FFT signal data is calculated using the following formula: ;in, This is FFT signal data; This is the conjugate transpose of the FFT signal data; The number of points in the FFT signal data; Let be the covariance matrix of the FFT signal data.

[0031] In step S4, Givens rotation is used to gradually eliminate the off-diagonal elements of the covariance matrix. The rotation matrix is ​​represented by the following formula:

[0032]

[0033] in, Let be the Givens rotation matrix for rotating between the i-th and j-th rows, where i is the row index of the element to be eliminated and j is the column index of the element to be eliminated. is the rotation angle calculated from the element to be eliminated; c is the first rotation factor, and s is the second rotation factor;

[0034] The CORDIC iterative optimization method is used to calculate the first twitch factor c and the second twitch factor s using a systolic array architecture, specifically as follows:

[0035] systolic array uses The structure contains several Processing Element (PE) units that perform Givens rotation calculations in parallel.

[0036] The systolic array calculation process is as follows: the first column PE calculates the Givens rotation of the first column data, eliminates the second diagonal elements of the first column, and flows the calculation result to the next column PE to continue eliminating the second diagonal elements of the second column, and so on, processing the subsequent columns in turn, and finally transforming the input matrix into an upper triangular matrix.

[0037] The CORDIC iterative optimization method is specifically as follows:

[0038] Pre-store the first rotation factor c and the second rotation factor s of several angles, and construct a lookup table;

[0039] When the required accuracy of the calculation result is lower than the angle interval in the pre-stored lookup table, the first rotation factor c and the second rotation factor s are obtained directly by looking up the table.

[0040] When the required accuracy of the calculation result is higher than the angle interval in the pre-stored lookup table, the error part is calculated by performing a binary angle approximation method in combination with the pre-stored lookup table to obtain the first rotation factor c and the second rotation factor s.

[0041] The finite iterative calculation also includes an error control strategy. The error is calculated, and when the error is less than a set threshold, the iterative calculation stops, and the final values ​​of c and s are taken as the results. The error e is calculated using the following formula: .

[0042] The present invention also provides a system for implementing the QR decomposition matrix inversion optimization method, including a data receiving and preprocessing module, a double-buffered FIFO data storage module, a covariance matrix calculation module, a Givens rotation calculation module, and a QR decomposition result output module;

[0043] The data receiving and preprocessing module receives signal data, performs FFT preprocessing and data downsampling to obtain FFT signal data, and uploads the data to the double-buffered FIFO data storage module;

[0044] The double-buffered FIFO data storage module adopts double-buffered optimized parallel processing. Based on the double-buffered FIFO data storage structure, it stores the received FFT signal data and uploads the data to the covariance matrix calculation module.

[0045] The covariance matrix calculation module calculates the covariance matrix of the FFT signal data based on the received data, uses it as the Givens rotation input data, and uploads the data to the Givens rotation calculation module.

[0046] The Givens rotation calculation module uses Givens rotation to gradually eliminate the off-diagonal elements of the covariance matrix based on the received data, performs Givens rotation calculation using a Systolic array architecture, calculates the Givens rotation factor using an optimized CORDIC calculation method, and uploads the data to the QR decomposition result output module.

[0047] The QR decomposition result output module outputs the final QR decomposition result based on the Givens rotation calculation result.

[0048] This invention discloses an optimized method and system for QR decomposition matrix inversion, which optimizes data flow, improves computational throughput, and reduces hardware resource consumption. It provides a low-power, high-performance hardware implementation scheme for efficient matrix computation and has broad engineering application value. Attached Figure Description

[0049] Figure 1 This is a schematic flowchart of the method of the present invention;

[0050] Figure 2 This is a schematic diagram of the system of the present invention. Detailed Implementation

[0051] This invention provides an optimization method for QR decomposition matrix inversion, the flowchart of which is shown below. Figure 1 As shown, it includes the following steps:

[0052] S1. Receive signal data, perform FFT preprocessing and data downsampling to obtain FFT signal data;

[0053] Step S1 includes the following steps:

[0054] Receive signal data, perform FFT calculations to obtain frequency point data;

[0055] Principal component analysis was performed on the frequency data to obtain the effective frequency data;

[0056] The effective frequency data is downsampled to obtain FFT signal data.

[0057] The principal component analysis specifically refers to:

[0058] Calculate the signal power spectrum of frequency point data Use the following formula to calculate: ;in, The complex spectrum value of the frequency point data obtained by FFT calculation; The complex number is the conjugate of the complex spectral value of the frequency point data;

[0059] Set the energy threshold again Select the frequency point where the cumulative energy percentage reaches the energy threshold, and express it using the following formula: ;in, This represents the final number of effective frequency points selected.

[0060] The downsampling ratio of the data is calculated using the following formula: .

[0061] When the source of the signal data to be processed is a GNSS signal or a Doppler frequency offset signal, the signal energy is unevenly distributed across different frequency bands. A non-uniform downsampling optimization scheme is adopted to replace the principal component analysis + downsampling scheme. The specific steps of the replaced step S1 are as follows:

[0062] Receive signal data, perform FFT calculations to obtain frequency point data;

[0063] Frequency-specific downsampling of frequency data is performed, and the non-uniform downsampling function within the frequency range is expressed by the following formula: ;in, The preset principal component threshold of the signal;

[0064] Set low-frequency and high-frequency ranges. In the low-frequency range, use the first frequency point retention ratio to retain frequencies, and in the high-frequency range, use the second frequency point retention ratio to retain frequencies; the first frequency point retention ratio is greater than the frequency point retention ratio.

[0065] Finally, the data downsampling ratio is calculated and the data is downsampled.

[0066] Specifically, in practical applications, the low-frequency region is 0~500MHz, with 80% of the frequency points retained; the high-frequency region is greater than 500MHz, with 20% of the frequency points retained.

[0067] Using FFT data downsampling reduces computation by 50%.

[0068] S2. Double-buffered optimized parallel processing is adopted, based on a double-buffered FIFO data storage structure, to store FFT signal data;

[0069] In step S2, the double-buffered FIFO data storage structure includes buffer A and buffer B;

[0070] The dual-buffered optimized parallel processing is as follows: First, buffer A is used to store the currently obtained FFT signal data; when buffer A is insufficient, the process switches to buffer B to store the currently obtained FFT signal data, while simultaneously starting the calculation of step S3 on the FFT signal data already stored in buffer A; when buffer B is insufficient and the subsequent calculation of the FFT signal data already stored in buffer A is completed, buffer A is cleared, and the process switches back to buffer A to store the currently obtained FFT signal data, while simultaneously starting the calculation of step S3 on the FFT signal data already stored in buffer B; the two buffers are alternately switched for storage and reading calculation according to the above steps to achieve parallel processing.

[0071] By using a double-buffered FIFO data storage structure, pipelined CORDIC computation is realized, which not only improves parallelism by 2 times, but also reduces memory access and power consumption.

[0072] S3. Calculate the covariance matrix of the FFT signal data and use it as the Givens rotation input data;

[0073] In step S3, the covariance matrix of the FFT signal data is calculated using the following formula: ;in, This is FFT signal data; This is the conjugate transpose of the FFT signal data; The number of points in the FFT signal data; Let be the covariance matrix of the FFT signal data.

[0074] S4. The Givens rotation is performed by gradually eliminating the off-diagonal elements of the covariance matrix, and the Givens rotation factor is calculated using an optimized CORDIC calculation method. The Givens rotation calculation process is accelerated by using a Systolic array architecture.

[0075] In step S4, Givens rotation is used to gradually eliminate the off-diagonal elements of the covariance matrix. The rotation matrix is ​​represented by the following formula:

[0076]

[0077] in, Let be the Givens rotation matrix for rotating between the i-th and j-th rows, where i is the row index of the element to be eliminated and j is the column index of the element to be eliminated. is the rotation angle calculated from the element to be eliminated; c is the first rotation factor, and s is the second rotation factor;

[0078] The CORDIC iterative optimization method is used to calculate the first twitch factor c and the second twitch factor s using a systolic array architecture, specifically as follows:

[0079] systolic array uses The structure contains several Processing Element (PE) units that perform Givens rotation calculations in parallel.

[0080] The systolic array calculation process is as follows: the first column PE calculates the Givens rotation of the first column data, eliminates the second diagonal elements of the first column, and flows the calculation result to the next column PE to continue eliminating the second diagonal elements of the second column, and so on, processing the subsequent columns in turn, and finally transforming the input matrix into an upper triangular matrix.

[0081] This invention adopts Parallel computing using a systolic array structure improves throughput by 3 times.

[0082] The CORDIC iterative optimization method is specifically as follows:

[0083] Pre-store the first rotation factor c and the second rotation factor s of several angles, and construct a lookup table;

[0084] When the required accuracy of the calculation result is lower than the angle interval in the pre-stored lookup table, the first rotation factor c and the second rotation factor s are obtained directly by looking up the table.

[0085] When the required accuracy of the calculation result is higher than the angle interval in the pre-stored lookup table, the error part is calculated by performing a binary angle approximation method in combination with the pre-stored lookup table to obtain the first rotation factor c and the second rotation factor s.

[0086] The finite iterative calculation also includes an error control strategy. The error is calculated, and when the error is less than a set threshold, the iterative calculation stops, and the final values ​​of c and s are taken as the results. The error e is calculated using the following formula: .

[0087] This invention employs a method that combines lookup tables with CORDIC calculation, reducing the number of CORDIC iterations by 50%.

[0088] S5. Based on the Givens rotation calculation results, output the final QR decomposition results.

[0089] The present invention also provides a system for implementing the QR decomposition matrix inversion optimization method, the schematic diagram of which is shown below. Figure 2 As shown, it includes a data receiving and preprocessing module, a double-buffered FIFO data storage module, a covariance matrix calculation module, a Givens rotation calculation module, and a QR decomposition result output module;

[0090] The data receiving and preprocessing module receives signal data, performs FFT preprocessing and data downsampling to obtain FFT signal data, and uploads the data to the double-buffered FIFO data storage module;

[0091] The double-buffered FIFO data storage module adopts double-buffered optimized parallel processing. Based on the double-buffered FIFO data storage structure, it stores the received FFT signal data and uploads the data to the covariance matrix calculation module.

[0092] The covariance matrix calculation module calculates the covariance matrix of the FFT signal data based on the received data, uses it as the Givens rotation input data, and uploads the data to the Givens rotation calculation module.

[0093] The Givens rotation calculation module uses Givens rotation to gradually eliminate the off-diagonal elements of the covariance matrix based on the received data, performs Givens rotation calculation using a Systolic array architecture, calculates the Givens rotation factor using an optimized CORDIC calculation method, and uploads the data to the QR decomposition result output module.

[0094] The QR decomposition result output module outputs the final QR decomposition result based on the Givens rotation calculation result.

[0095] The method of the present invention will be further described below with reference to an embodiment:

[0096] In this embodiment, a public dataset, a standard random 4×4 complex matrix set, with 1000 samples, is selected as the test input. First, the input matrix is ​​processed by FFT, and then covariance estimation and QR decomposition are performed according to the method of this invention.

[0097] On a Xilinx Kintex-7 XC7K325T FPGA, implemented at a clock frequency of 200MHz: using the method of this invention, the average relative error of the diagonal of the QR decomposition matrix R does not exceed 0.1%, and the residual suppression ratio is greater than 65dB; the average single-frame processing latency is about 2.3μs, corresponding to a throughput of about 4.2Msps; the total power consumption is estimated to be about 1.2W by Xilinx XPower Analyzer.

[0098] A comparison was made using existing conventional QR decomposition methods. Based on traditional Givens rotation, this method requires calling division and square root units, and is processed serially. The comparison was conducted under the same dataset and hardware conditions. Results show: the average relative error on the diagonal is approximately 0.2%, the residual suppression ratio is less than 50dB; the average single-frame processing latency is approximately 6.5μs, corresponding to a throughput of approximately 1.5Msps; the total power consumption is estimated at approximately 1.8W.

[0099] As can be seen, compared with existing conventional QR decomposition methods, the method of the present invention not only has a lower relative error and a higher throughput, but also reduces the total power consumption.

Claims

1. An optimization method for QR decomposition matrix inversion, characterized in that, Includes the following steps: S1. Receive signal data, perform FFT preprocessing and data downsampling to obtain FFT signal data; S2. Double-buffered optimized parallel processing is adopted, based on a double-buffered FIFO data storage structure, to store FFT signal data; S3. Calculate the covariance matrix of the FFT signal data and use it as the Givens rotation input data; S4. The off-diagonal elements of the covariance matrix are gradually eliminated by Givens rotation, Givens rotation calculation is performed, and the Givens rotation factor is calculated using the optimized CORDIC calculation method. The Systolic array architecture is used to accelerate the calculation process of Givens rotation; S5. Based on the Givens rotation calculation results, output the final QR decomposition results; Step S1 includes the following steps: Receive signal data, perform FFT calculations to obtain frequency point data; Principal component analysis was performed on the frequency data to obtain the effective frequency data; The effective frequency data is downsampled to obtain FFT signal data; The principal component analysis specifically refers to: Calculate the signal power spectrum of frequency point data Use the following formula to calculate: ;in, The complex spectrum value of the frequency point data obtained by FFT calculation; The complex number is the conjugate of the complex spectral value of the frequency point data; Set the energy threshold again Select the frequency point where the cumulative energy percentage reaches the energy threshold, and express it using the following formula: ;in, The final number of effective frequency points selected; The downsampling ratio of the data is calculated using the following formula: .

2. The QR decomposition matrix inversion optimization method according to claim 1, characterized in that, When the source of the signal data to be processed is a GNSS signal or a Doppler frequency offset signal, step S1 is as follows: Receive signal data, perform FFT calculations to obtain frequency point data; Frequency-specific downsampling of frequency data is performed, and the non-uniform downsampling function within the frequency range is expressed by the following formula: ;in, The preset principal component threshold of the signal; Set low-frequency and high-frequency ranges. In the low-frequency range, use the first frequency point retention ratio to retain frequencies, and in the high-frequency range, use the second frequency point retention ratio to retain frequencies; the first frequency point retention ratio is greater than the frequency point retention ratio. Finally, the data downsampling ratio is calculated and the data is downsampled.

3. The QR decomposition matrix inversion optimization method according to claim 1, characterized in that, In step S2, the double-buffered FIFO data storage structure includes buffer A and buffer B; The dual-buffered optimized parallel processing is as follows: First, buffer A is used to store the currently obtained FFT signal data; when buffer A is insufficient, the process switches to buffer B to store the currently obtained FFT signal data, while simultaneously starting the calculation of step S3 on the FFT signal data already stored in buffer A; when buffer B is insufficient and the subsequent calculation of the FFT signal data already stored in buffer A is completed, buffer A is cleared, and the process switches back to buffer A to store the currently obtained FFT signal data, while simultaneously starting the calculation of step S3 on the FFT signal data already stored in buffer B; the two buffers are alternately switched for storage and reading calculation according to the above steps to achieve parallel processing.

4. The QR decomposition matrix inversion optimization method according to claim 1, characterized in that, In step S3, the covariance matrix of the FFT signal data is calculated using the following formula: ;in, This is FFT signal data; This is the conjugate transpose of the FFT signal data; The number of points in the FFT signal data; Let be the covariance matrix of the FFT signal data.

5. The QR decomposition matrix inversion optimization method according to claim 1, characterized in that, In step S4, Givens rotation is used to gradually eliminate the off-diagonal elements of the covariance matrix. The rotation matrix is ​​represented by the following formula: in, Let be the Givens rotation matrix for rotating between row i and row j, where i is the row index of the element to be eliminated and j is the column index of the element to be eliminated. is the rotation angle calculated from the element to be eliminated; c is the first rotation factor, and s is the second rotation factor; The CORDIC iterative optimization method is used to calculate the first twitch factor c and the second twitch factor s using a systolic array architecture, specifically as follows: systolic array uses The structure contains several Processing Element (PE) units that perform Givens rotation calculations in parallel. The systolic array calculation process is as follows: the first column PE calculates the Givens rotation of the first column data, eliminates the second diagonal elements of the first column, and flows the calculation result to the next column PE to continue eliminating the second diagonal elements of the second column, and so on, processing the subsequent columns in turn, and finally transforming the input matrix into an upper triangular matrix. The CORDIC iterative optimization method is used to calculate the first rotation factor c and the second rotation factor s based on the obtained upper triangular matrix.

6. The QR decomposition matrix inversion optimization method according to claim 5, characterized in that, The CORDIC iterative optimization method is specifically as follows: Pre-store the first rotation factor c and the second rotation factor s of several angles, and construct a lookup table; When the required accuracy of the calculation result is lower than the angle interval in the pre-stored lookup table, the first rotation factor c and the second rotation factor s are obtained directly by looking up the table. When the required accuracy of the calculation result is higher than the angle interval in the pre-stored lookup table, the error part is approximated by binary angles using the pre-stored lookup table to perform finite iterative calculations, and the first rotation factor c and the second rotation factor s are obtained.

7. The QR decomposition matrix inversion optimization method according to claim 6, characterized in that, The finite iterative calculation also includes an error control strategy. The error is calculated, and when the error is less than a set threshold, the iterative calculation stops, and the final values ​​of c and s are taken as the results. The error e is calculated using the following formula: .

8. A system for implementing the QR decomposition matrix inversion optimization method according to any one of claims 1 to 7, characterized in that, It includes a data receiving and preprocessing module, a double-buffered FIFO data storage module, a covariance matrix calculation module, a Givens rotation calculation module, and a QR decomposition result output module; The data receiving and preprocessing module receives signal data, performs FFT preprocessing and data downsampling to obtain FFT signal data, and uploads the data to the double-buffered FIFO data storage module; The double-buffered FIFO data storage module adopts double-buffered optimized parallel processing. Based on the double-buffered FIFO data storage structure, it stores the received FFT signal data and uploads the data to the covariance matrix calculation module. The covariance matrix calculation module calculates the covariance matrix of the FFT signal data based on the received data, uses it as the Givens rotation input data, and uploads the data to the Givens rotation calculation module. The Givens rotation calculation module uses Givens rotation to gradually eliminate the off-diagonal elements of the covariance matrix based on the received data, performs Givens rotation calculation using a Systolic array architecture, calculates the Givens rotation factor using an optimized CORDIC calculation method, and uploads the data to the QR decomposition result output module. The QR decomposition result output module outputs the final QR decomposition result based on the Givens rotation calculation result.

Images