FPGA (Field Programmable Gate Array) implementation method of multiple signal classification algorithm based on real number processing
Through real-numeric processing, the complex covariance matrix is converted into the real covariance matrix, which solves the problem of large amount of calculation of the eigenvalue decomposition of the complex covariance matrix, and realizes the efficiency improvement of the multiple signal classification algorithm.
Patent Information
- Application Number
- CN202411889173.3
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2024-12-13
- Publication Date
- 2025-04-18
AI Technical Summary
In the existing array signal processing technology, the eigenvalue decomposition calculation of the complex covariance matrix is large, which limits the efficiency of the multiple signal classification algorithm.
Through real-numeric processing, the complex covariance matrix is converted into a real covariance matrix, and the real-numeric vector is used instead of the complex space search vector, reducing the calculation amount.
The calculation amount of multiple signal classification algorithm is reduced, hardware resource consumption is reduced, and computational convergence speed is improved.
Smart Images
Figure CN119848597A_ABST
Abstract
Description
Technical Field
[0001] The invention relates to the technical field of array signal processing, in particular to an FPGA implementation of a multiple signal classification algorithm based on real number processing. Background Art
[0002] As an important branch of signal processing, array signal processing has been widely used and developed rapidly in the fields of communication, radar, sonar, seismic exploration and radio astronomy. Array signal processing arranges a group of sensors in different positions in space in a certain way to form a sensor array. Using a sensor array to receive spatial signals is equivalent to sampling the spatially distributed field signals to obtain spatial discrete observation data of the signal source. The purpose of array signal processing is to enhance the required useful signals, suppress useless interference and noise, and extract useful signal features and information contained in the signals by processing the signals received by the array. The traditional direction finding algorithm that appeared in the early days has great limitations. Its resolution is limited by the Rayleigh limit. The so-called Rayleigh limit is the limitation of the antenna aperture size. Compared with the traditional direction finding algorithm, the multiple signal classification algorithm has many advantages, such as: higher resolution, accurately distinguishing multiple very close signal sources at a shorter array antenna spacing or a smaller number of array elements; stronger anti-noise performance, can be applied to complex working environments; has the ability to process multiple signal sources, and can process multiple signal sources arriving at the same time. In the subspace algorithm, the eigenvalue decomposition of the complex covariance matrix is a crucial step, but the computational complexity of this step is very high. Through real number processing, the original complex covariance matrix is converted into a real covariance matrix. At the same time, the complex space search vector in the traditional algorithm is also replaced by a real vector, which reduces the computational complexity of subsequent processing. Summary of the invention
[0003] The purpose of the present invention is to provide an FPGA implementation method of a multiple signal classification algorithm based on real number processing.
[0004] The technical solution for achieving the purpose of the present invention is: an FPGA implementation method of a multiple signal classification algorithm based on real number processing, comprising the following steps:
[0005] Step 1: The N-element linear array antenna receives the incident signal, which is composed of P mutually unrelated narrowband signals and noise in space. The N-element uniform linear array sampling data is output to the covariance matrix calculation module, and N and P are positive integers, P <N;
[0006] Step 2: The covariance matrix calculation module calculates the covariance matrix according to the input N-element uniform linear array sampling data and outputs the complex covariance matrix R.
[0007] Step 3: The covariance matrix realization module performs covariance matrix realization based on the input complex covariance matrix R, converts the complex matrix calculation to real matrix calculation, and outputs the real covariance matrix
[0008] Step 4: The information source number estimation module outputs the matrix R′ to the eigendecomposition module for eigendecomposition. The eigendecomposition module outputs the eigenvector U to the information source number estimation module for information source number estimation. The information source number estimation module outputs the calculated number of information sources.
[0009] Step 5: The eigendecomposition module is based on the input real covariance matrix Perform eigendecomposition and output real covariance matrix The feature vector of .
[0010] Step 6: The peak search module uses the input real covariance matrix The eigenvector and the number of signal sources are used to search for spectral peaks and output the spatial spectrum function.
[0011] Based on the above basic scheme, the array model of step one is:
[0012] A uniform linear array consisting of M array elements receives K far-field narrowband signals, assuming that the number of array elements is equal to the number of channels. According to the narrowband array model, the linear combination of these K incident signals and noise can express the data received by the array. The first array element of the array is selected as the reference element, and the data matrix of the array receiving signal is:
[0013] X(t)=AS(t)+N(t)
[0014] Where X(t) is the snapshot data; S(t) = [S 1 (t), S 2 (t), ..., S K (t)] T is the spatial signal; A is the steering vector matrix; N(t)=[N 1 (t), N 2 (t), ..., N M (t)] T is the noise vector matrix.
[0015] Its covariance matrix is:
[0016] R=E[XX H ]=AE[SS H ]A H +E[NN H ]
[0017] =AR s A H +R N=AR s A H +δ 2 I
[0018] Step 2 The working process of covariance matrix calculation is:
[0019]
[0020] Among them, x i is the sampling data of array element i, there are N array elements in total, and the superscript H is the conjugate transpose.
[0021]
[0022] where x 1 , x 2 ,...x N is the sampling data of array elements 1, 2, ...N, is the conjugate of the sampling data of array elements 1, 2, ...N.
[0023] Step 3 The working process of calculating the real number of the covariance matrix is:
[0024] Introduce a unitary matrix Q as:
[0025]
[0026] When the number of array elements N is an even number, I represents the N / 2-order unit matrix, and J represents the N / 2-order antisymmetric matrix.
[0027] For the covariance matrix calculated from a finite sample, the covariance matrix obtained by sampling is It satisfies Hermitian but not supersymmetric, and adopts a supersymmetric estimation method of sampling covariance matrix:
[0028]
[0029] in It satisfies the Hermitian and supersymmetric properties, that is, and
[0030]
[0031] In the sense of Euclidean distance, is the best Hermitian supersymmetry estimate,
[0032]
[0033] For any Hermitian supersymmetric matrix A in the above equation, the Euclidean norm is defined as ||B|| 2 =tr(B H B), since Q is a unitary matrix, Equivalent to Because QAQ H real number, if there exists an optimal Hermitian supersymmetric matrix A 0 , so that Minimum, must meet
[0034]
[0035] Re(·) means taking the real part, multiply the above formula by Q on the left H , multiply Q on the right, according to QQ H =J, yes
[0036]
[0037] From the above formula, we can conclude that: for It is the optimal Hermitian supersymmetric estimate in terms of Euclidean distance. We can obtain:
[0038]
[0039] So calculate Not necessary, just calculate Use its eigenvector to estimate the orientation.
[0040] Remember The matrix after unitary transformation is For the real covariance matrix after unitary transformation The eigenvalue of The eigenvalues are the same, and the eigenvector relationship between them satisfies the following relationship:
[0041]
[0042] in yes The characteristic vector of yes The feature vector of .
[0043] When the number of array elements N is an even number, I represents the N / 2-order unit matrix, J represents the N / 2-order antisymmetric matrix, and the covariance matrix R is unitarily transformed. The real part after the unitary transformation is:
[0044]
[0045] This is the real-numbered covariance matrix, which is the required covariance matrix for subsequent calculations. In the calculation of the covariance matrix, only matrix multiplication and addition operations are involved, and since the matrix elements of the unitary matrix are 0 and 1, only addition and subtraction of matrix elements are required in the actual calculation process.
[0046] The working process of step 4 source number estimation is as follows:
[0047] The real covariance matrix To split into blocks:
[0048]
[0049] Extract the elements of the first N-1 rows and columns of the covariance matrix to get R′, perform eigendecomposition on R′ to get the eigenvector matrix U, and then construct the transformation matrix T:
[0050]
[0051] After transforming the covariance matrix, it is:
[0052]
[0053] The judgment rule according to the Gale circle is:
[0054]
[0055] When GDE(n) is less than 0 for the first time, record the value of n this time, then the number of information sources is P=n-1.
[0056] The working process of feature decomposition in steps 4 and 5 is:
[0057] The eigendecomposition of the covariance matrix is realized by the parallel Jaocbi algorithm, and the covariance matrix is divided into diagonal units and non-diagonal units, where the diagonal unit is a matrix divided by the main diagonal elements of the covariance matrix, and the rest are non-diagonal units, each of which is a 2nd-order matrix. The rotation matrix is calculated through the diagonal units. First, the rotation angle θ needs to be calculated. For the following diagonal units:
[0058]
[0059] The calculation formula is:
[0060]
[0061] From this, we can get the rotation matrix elements sinθ and cosθ, and the rotation matrix is:
[0062]
[0063] After the rotation matrix is calculated by the diagonal unit, it is passed to the non-diagonal unit for rotation calculation, and the rotation calculation is also performed on itself. The formula for rotation calculation is as follows:
[0064]
[0065] After rotating the diagonal cells, a′ in the rotated matrix N，N-1 , a′ N-1，N All are 0. When the diagonal and non-diagonal units have completed the rotation calculation, it is necessary to perform row and column scheduling, that is, to exchange the elements in the matrix, in preparation for the next calculation of the rotation matrix and the rotation calculation. After completing the calculation of the eigenvalues and eigenvectors, it is also necessary to arrange the eigenvalues and corresponding eigenvectors in descending order.
[0066] 6. The FPGA implementation method of the multiple signal classification algorithm based on real number processing according to claim 1 is characterized in that the working process of the spectrum peak search in step 6 is:
[0067] The real-numbered covariance matrix of step three is subjected to eigendecomposition in step five to obtain the eigenvector matrix. Then, according to the number of signal sources obtained by the signal source number estimation module in step four, the noise subspace in the eigenvector matrix is extracted. The direction vector in the spatial spectrum function also needs to be real-numbered, pre-calculated in MALTAB, and then stored in ROM. The formula of the spatial spectrum function of the MUISC algorithm is as follows.
[0068]
[0069] Since spectrum peak search requires inverse calculation, which is not easy to do in FPGA, in order to reduce the complexity of calculation, P can be calculated. MUSIC The reciprocal of (θ).
[0070] The present invention greatly reduces the amount of calculation of the multiple signal classification algorithm through a real number processing, converts the original complex number operation into a real number operation, and achieves the desired estimation effect. Compared with the conventional implementation method, the real number processing can reduce the consumption of hardware resources and improve the calculation convergence speed. BRIEF DESCRIPTION OF THE DRAWINGS
[0071] Figure 1 A signal processing flow chart of an embodiment of the present invention
[0072] Figure 2 Structural diagram for covariance matrix processing
[0073] Figure 3 The structural diagram for the real number processing of complex matrix
[0074] Figure 4 The structural diagram of the source number estimation
[0075] Figure 5 Flowchart for matrix eigendecomposition
[0076] Figure 6 The structure diagram of the characteristic decomposition row and column scheduling
[0077] Figure 7 Structure diagram for peak search
[0078] Figure 8 Simulation comparison of the classic multiple signal classification algorithm and the real-number processed MUSIC algorithm DETAILED DESCRIPTION
[0079] The method of the present invention is further described below in conjunction with the accompanying drawings and embodiments.
[0080] The basic scheme of the present invention includes the following six steps, and the specific signal processing flow chart is as follows: Figure 1 shown.
[0081] Step 1: An N-element linear array antenna receives an incident signal, which is composed of P mutually unrelated narrowband signals and noise in space. The N-element uniform linear array sampling data is output to a covariance matrix calculation module, where N and P are positive integers and P < N.
[0082] Step 2: The covariance matrix calculation module calculates the covariance matrix according to the input N-element uniform linear array sampling data and outputs the complex covariance matrix R.
[0083] Step 3: The covariance matrix realization module performs covariance matrix realization according to the input complex covariance matrix R, converts the complex matrix calculation to real matrix calculation, and outputs the real covariance matrix
[0084] Step 4: The information source number estimation module outputs the matrix R′ to the eigendecomposition module for eigendecomposition. The eigendecomposition module outputs the eigenvector U to the information source number estimation module for information source number estimation. The information source number estimation module outputs the calculated number of information sources.
[0085] Step 5: The eigendecomposition module is based on the input real covariance matrix Perform eigendecomposition and output real covariance matrix The feature vector of .
[0086] Step 6: The peak search module uses the input real covariance matrix The eigenvector and the number of signal sources are used to search for spectral peaks and output the spatial spectrum function.
[0087] Establish the signal model. A uniform linear array composed of M array elements receives K far-field narrowband signals. Assume that the number of array elements is equal to the number of channels. According to the narrowband array model, the linear combination of these K incident signals and noise can express the data received by the array. Select the first array element as the reference array element, then the data matrix of the array receiving signal is:
[0088] X(t)=AS(t)+N(t)
[0089] according to Figure 2 Compute the covariance matrix:
[0090]
[0091] Among them, x i is the sampling data of array element i, there are N array elements in total, and the superscript H is the conjugate transpose.
[0092]
[0093] where x 1 , x 2 ,...x N is the sampling data of array elements 1, 2, ..., N, is the conjugate of the sampling data of array elements 1, 2, ..., N. Construct N complex multipliers for calculation. During the complex multiplication operation, the multipliers used for each row of the matrix are time-division multiplexed. This method can save multiplier resources and accelerate timing convergence. Since the array element sampling data is also input by row, this method will not block the pipeline. The calculated covariance matrix is cumulatively summed up in snapshots, and after setting the number of snapshots, the average is calculated and then output. Since the covariance matrix is input in parallel per row, each row is corresponded to a submodule, and a total of N identical submodules are instantiated for calculation.
[0094] according to Figure 3 To process the real numbers:
[0095] When the number of array elements N is an even number, I represents the N / 2-order unit matrix, J represents the N / 2-order antisymmetric matrix, and the covariance matrix R is unitarily transformed. The real part after the unitary transformation is:
[0096]
[0097] This is the real-numbered covariance matrix, which is the required covariance matrix for subsequent calculations. In the calculation of the covariance matrix, only matrix multiplication and addition operations are involved, and since the matrix elements of the unitary matrix are 0 and 1, only addition and subtraction of matrix elements are required in the actual calculation process.
[0098] according to Figure 4 Estimation of the number of sources:
[0099] First, select the first 7 rows and 7 columns of the covariance matrix, then fill the elements on the diagonal with 1 and the remaining positions with 0 to turn it into an 8-order matrix. Perform eigendecomposition on the matrix. In addition to the eigenvalue 1 and the corresponding eigenvector, the rest is the eigenvector matrix required for source number estimation. Then, reorganize the obtained eigenvector matrix to obtain the unitary transformation matrix T, perform unitary transformation on the covariance matrix, and design a division operation in obtaining the radius of the Gale circle. Since this design uses 8 array elements, according to the following formula:
[0100]
[0101] The divisor is 7, so the divider method can be used. This method consumes DSP. The load can also be set to 7*x. The value needs to be between 0 and 1, and x should be a negative power of 2. In this way, you only need to shift the required data to get the result. This design uses the shift method. Since the real covariance matrix is obtained first and then the eigendecomposition is performed, the input covariance matrix is first stored in RAM. After the eigendecomposition is completed to obtain the unitary matrix, subsequent calculations are performed. Since the unitary transformation requires two unitary transformation matrices T, they are stored in RAM when input. When the Gale circle estimation is performed, when the value of GDE is obtained, it is only necessary to sum the sign bits of all its numbers to obtain the source number.
[0102] according to Figure 5 , 6 Perform matrix eigendecomposition:
[0103] The eigendecomposition of the covariance matrix is realized by the parallel Jaocbi algorithm, and the covariance matrix is divided into diagonal units and non-diagonal units, where the diagonal unit is a matrix divided by the main diagonal elements of the covariance matrix, and the rest are non-diagonal units, each of which is a 2nd-order matrix. The rotation matrix is calculated through the diagonal units. First, the rotation angle θ needs to be calculated. For the following diagonal units:
[0104]
[0105] The calculation formula is:
[0106]
[0107] From this, we can get the rotation matrix elements sinθ and cosθ, and the rotation matrix is:
[0108]
[0109] After the rotation matrix is calculated by the diagonal unit, it is passed to the non-diagonal unit for rotation calculation, and the rotation calculation is also performed on itself. The formula for rotation calculation is as follows:
[0110]
[0111] After rotating the diagonal cells, a′ in the rotated matrix N,N-1 , a′ N-1,N All are 0. When the rotation calculation of diagonal units and non-diagonal units is completed, it is necessary to perform row and column scheduling, that is, to exchange elements in the matrix, in preparation for the next calculation of the rotation matrix and rotation calculation. Dividing the diagonal and non-diagonal units, calculating the rotation matrix, performing rotation calculation, and performing row and column scheduling is an iterative process of the parallel Jacobi algorithm. The non-diagonal elements in the matrix are continuously eliminated by rotation, so the eigenvalue matrix is calculated. FPGA will produce truncation errors during fixed-point calculations. Therefore, the error caused by rotation iteration for smaller numbers is relatively large. It is necessary to select a suitable threshold to limit the rotation process. The upper triangular matrix cannot be all zero as a judgment condition. When the iteration reaches a certain number of times, the iteration can be stopped, otherwise the error will not decrease but increase. After the above iterative calculation of the covariance matrix, the eigenvalue matrix can be obtained. For the calculation of the eigenvector matrix, the same iterative process needs to be performed on the unit matrix. After completing the calculation of the eigenvalue and eigenvector, the eigenvalue and the corresponding eigenvector need to be arranged in order from large to small.
[0112] according to Figure 7 To perform a peak search:
[0113] The real covariance matrix is subjected to eigendecomposition in step five to obtain the eigenvector matrix. Then, based on the number of signal sources obtained by the signal source number estimation module, the noise subspace in the eigenvector matrix is extracted. The direction vector in the spatial spectrum function also needs to be realized and pre-calculated in MALTAB and then stored in ROM. The formula of the spatial spectrum function of the MUISC algorithm is as follows.
[0114]
[0115] Since spectrum peak search requires inverse calculation, which is not easy to do in FPGA, in order to reduce the complexity of calculation, P can be calculated. MUSIC The reciprocal of (θ).
[0116] Figure 8 This is a simulation comparison chart of the classic multiple signal classification algorithm and the algorithm after real number processing. It can be seen that after real number processing, the expected estimation effect can also be achieved.
Claims
1. An FPGA implementation method of a multiple signal classification algorithm based on real number processing comprises the following steps: Step 1: The N-element linear array antenna receives the incident signal, which is composed of P uncorrelated narrowband signals and noise in space. The N-element uniform linear array sampling data is output to the covariance matrix calculation module, and N and P are positive integers, P <N; Step 2: The covariance matrix calculation module calculates the covariance matrix according to the input N-element uniform linear array sampling data and outputs the complex covariance matrix R. Step 3: The covariance matrix realization module performs covariance matrix realization based on the input complex covariance matrix R, converts the complex matrix calculation to real matrix calculation, and outputs the real covariance matrix Step 4: The information source number estimation module outputs the matrix R′ to the eigendecomposition module for eigendecomposition. The eigendecomposition module outputs the eigenvector U to the information source number estimation module for information source number estimation. The information source number estimation module outputs the calculated number of information sources. Step 5: The eigendecomposition module is based on the input real covariance matrix Perform eigendecomposition and output real covariance matrix The feature vector of . Step 6: The peak search module uses the input real covariance matrix The eigenvector and the number of signal sources are used to search for spectral peaks and output the spatial spectrum function.
2. The FPGA implementation method of the multiple signal classification algorithm based on real number processing according to claim 1 is characterized in that: The working process of the covariance matrix calculation module in step 2 is as follows: First, receive N-element uniform linear array data and obtain the autocorrelation matrix: Among them, x i is the sampling data of array element i, there are N array elements in total, and the superscript H is the conjugate transpose. where x1, x2, ...x N is the sampling data of array elements 1, 2, ..., N, is the conjugate of the sampling data of array elements 1, 2, …, N; finally, the complex covariance matrix R is calculated.
3. The FPGA implementation method of the multiple signal classification algorithm based on real number processing according to claim 1 is characterized in that: The working process of the calculation of the covariance matrix real number module in step 3 is: First, receive the complex covariance matrix R from the covariance matrix calculation module and introduce a unitary matrix Q as: When the number of array elements N is an even number, I represents the N / 2-order unit matrix, and J represents the N / 2-order antisymmetric matrix; then the complex covariance matrix R is transformed as follows: Finally, the real covariance matrix is calculated 4. The FPGA implementation method of the multiple signal classification algorithm based on real number processing according to claim 1 is characterized in that: The working process of the information source number estimation module in step 4 is as follows: First accept the real covariance matrix To split into blocks: Secondly, for the real covariance matrix The elements of the first N-1 rows and columns of are extracted to obtain the matrix R′, and the matrix R′ is output to the eigendecomposition module for eigendecomposition to obtain the eigenvector matrix U, and then the transformation matrix T is constructed: Then the real covariance matrix After transformation: The judgment rule according to the Gale circle is: Finally, when GDE(n) is less than 0 for the first time, the value of n is recorded and the number of information sources P=n-1 is calculated.
5. The FPGA implementation method of the multiple signal classification algorithm based on real number processing according to claim 1 is characterized in that: The working process of the feature decomposition module in step 5 is as follows: First accept the real covariance matrix It is divided into diagonal units and non-diagonal units, where the diagonal units are matrices divided by the main diagonal elements of the covariance matrix, and the rest are non-diagonal units. Each unit is a 2nd-order matrix as shown below: Next, calculate the rotation angle θ. For the diagonal unit in the above formula, the calculation formula is: According to the rotation angle θ, the rotation matrix elements sinθ and cosθ can be calculated, and the obtained rotation matrix is: Again, after the rotation matrix is calculated by the diagonal unit, it is passed to the non-diagonal unit for rotation calculation, and the rotation calculation is also performed on itself. The rotation calculation formula is as follows: After rotating the diagonal cells, a′ in the rotated matrix N,N-1 , a′ N-1，N are all 0; then when the diagonal units and non-diagonal units have completed the rotation calculation, it is necessary to perform row and column scheduling, that is, to exchange the elements in the matrix to prepare for the next calculation of the rotation matrix and the rotation calculation; finally, after completing the calculation of the eigenvalues and eigenvectors, arrange the eigenvalues and the corresponding eigenvectors in order from large to small, and output the real covariance matrix The feature vector of .
6. The FPGA implementation method of the multiple signal classification algorithm based on real number processing according to claim 1 is characterized in that: The working process of the spectrum peak search module in step 6 is: First, receive the source number and the real covariance matrix The eigenvector of is calculated according to the following formula: Finally, the calculated spatial spectrum function is output.