KR102026958B1 - reduced computational complexity orthogonal matching pursuit using partitioned inversion technique for compressive sensing - Google Patents

Abstract

Disclosed is an OMP method for processing a signal received by a radar by using compressive sensing, and performing an OMP algorithm for the compressive sensing, more specifically, an OMP method for compressive sensing for a radar, which uses partitioned inversion in inverse matrix calculation for application of the least-squares method. The partitioned inversion includes a process of taking an i^th partial matrix when i is an integer to perform i^th inverse matrix calculation in a previous step, and using the result of the i^th inverse matrix calculation obtained in the previous step while taking an (i+1)^th partial matrix including the i^th partial matrix taken in the previous step and performing an (i+1)^th inverse matrix calculation in a next step of the (i+1)^th. Thus, the computational complexity can be reduced.

Description

Reduced computational complexity orthogonal matching pursuit using partitioned inversion technique for compressive sensing

The present invention relates to an OMP (OMP) method, and more particularly to an OMP method that can reduce the computational complexity by using a split inverse transform technique for compression sensing used in radar signal processing.

Radar is a device that emits a radio wave signal of a certain bandwidth, detects radio waves reflected from an object, obtains information from the reflected wave signal, and checks an object based on the information or identifies a feature of the object. In this case, a wideband wave must be used to accurately determine an object, and a short period radio wave must be used to detect a small object over a long distance.

In order to represent an object by converting an analog signal related to the object into a digital signal and restoring it into an analog form, a classic Nyquist Shannon sample theorem according to the classical Nyquist Shannon sample theorem acquires a discrete signal at a sampling rate of more than twice the maximum frequency or band. Sampling can restore the signal without distortion or loss. To accurately reconstruct the object-related signal detected by the radar, an analog-to-digital converter with a high sample rate of up to twice the frequency should be used to process the detected signal.

However, radar uses wideband signals, including very high frequency bands, so reconstructing the signal using transducers with adequate accuracy in real time requires transducers and computational resources with good capacity and computational power, but at the moment It is difficult to equip with a converter and a computer.

One way to solve this problem is to divide the signal bandwidth and process the signal using a corresponding converter and computer for each bandwidth, and combine multiple converters and computers used for the full-band signal to achieve overall signal recovery. However, even in this case, a very large number of converters and computation resources and time are required. For example, when a transducer and a computer are provided for the radar signal processing of a fighter aircraft, it requires too much weight and space, and thus practical application is difficult.

By the way, the Nyquist Shannon sample theorem has been understood to be accurate for a very long time, but because it is a very general theory of arbitrary baseband signals, in recent years, lossless restoration with less discrete signals is required when restoring object information for special types of signals. It is understood that it is possible.

An important signal type is the sparse signal. In the case of using wavelet or transform code for image compression, even if the actual image has an even distribution in all areas, only a few coefficients in the wavelet transform have very large values. It is known that images can be stored or transmitted at high compression ratios. For example, if one million pixels constituting one image have a significant value only in 25000 pixels, the image obtained by replacing the signal values of the remaining 25,000 pixels with 0 and replacing the signal values of the remaining pixels with 0 and then inverting again Looks pretty much the same as the video.

 At this time, the wavelet signal in which most of the pixels are replaced with 0 except for a part is an interspersed signal. Most images are approximated with scattered signals through proper conversion. This allows for a very effective representation of the video signal through scattered signals. The simplest example of an interspersed signal can also be seen through the Fourier transform. In the time domain, the constant signal is not an interspersed signal because it has a nonzero value for the entire time, but in the frequency domain it has a nonzero coefficient only if the frequency is zero, and the coefficients of all other frequencies are all zero. This is the scattering signal at.

In this example, it can be seen that an arbitrary signal can be a scattered signal through some transformation. Of course, not all signals have this feature, but it is well known that many signals (images, etc.) observed in nature can be scattered through some transformation. In general, interspersed signals are not observed or acquired as is. In the case of the above-mentioned images, the video signals we see are not scattered signals. Therefore, the problem of finding such scattered signals is important when a given signal is assumed to be scattered through some transformation. That is, the problem of finding the scattered signal from the observed signal is important, and recent research has shown that through the use of compressive sensing (CS), the scattered signal can be found by using only a small number of samples smaller than the length of the observed signal.

The CS can compress / restore the signal with very few samples when applied directly to an image signal. Therefore, we can obtain a new method that is very efficient from the traditional method of video signal compression and reconstruction, which obtains digital signal according to Nyquist Shannon sampling rate, converts wavelet, etc., and collects and stores and compresses only the largest coefficient. have. For example, when considering a one-dimensional signal, if you sample a 1 second signal with a bandwidth of 1 KHz at the Nyquist rate, you need 2000 samples, and only 20 coefficients to represent the sample signal through some transformation. Assuming a very large value, and unlike the traditional method, if we can get 20 coefficients from the scattered signal transformed into just 40 samples in one second of signal, then the latter method will be an efficient one. The theoretical basis of this new approach is CS.

Therefore, when the CS technique is used for radar signal processing, only the converter and the calculation equipment which are inferior to the conventional one can exhibit the same radar performance, and the same equipment can have much higher radar performance.

On the other hand, in order to understand CS, we first need to define a scattering signal, that is, a sparse signal. Assume an arbitrary N-dimensional linear space. When the scattering signal is expressed as a vector of this space, if the value of the nonzero coefficient of this vector is K (<< N) which is much smaller than N, then the signal or vector is a K-sparse signal or K- It is called an interstitial signal.

For example, if N = 100 and K = 1 0, then the K-pulse signal is a number where 10 coefficients are non-zero, and the remaining 90 coefficients are all zero. related to basis. The signal may or may not be scattered depending on the basis or transform used to represent the signal in some vector space. In the time domain, for example, constant signals have nonzero values for all times. These signals are not interspersed. But if you look at this signal in the frequency domain through the Fourier transform, it is an interspersed signal in the frequency domain where the frequency has a non-zero coefficient value for zero, and the coefficients for all other frequencies are all zeros. The same phenomenon can be seen with the (cosine) waveform.

To see this in more detail, define the following signal:

Where χ is the scattering signal vector of length N × 1. A is the M × N sensing matrix. Where M <N. In this case, the problem of determining χ by y is an underdetermined system where the number of equations is less than the number of variables, and the number of solutions is infinitely large. As a simple example, in the case of “10 = χ ₁ + χ ₂ ”, the solutions of both variables are infinitely large, but if χ is a K-pulse signal, this problem can be seen differently. Suppose we simply know where K = M and the nonzero values in χ. If you write the system again using a matrix that removes the row vector of A matrix corresponding to 0 (marked horizontal line on A and writes 'A'), it is as follows.

Where χ bar is a χ vector of only nonzero coefficients. In this case, the size of the A bar is M × M, so it is easy to determine the χ bar. For example, M = 3, N = 4, and K = 2.

The first system is a deficient system with three variables in four equations, but two of the variables (that is, if χ ₂ and χ ₃ are zero) are redundant or overdetermined systems with two variables in three equations.

Through this simple example, the general signal as in Equation 1 can not be solved, but it can be seen that the scattered signal is different. However, in order to move on to Equation 2, there must be an answer about how to obtain the A bar by finding a position where the coefficient of the χ vector is not zero. Before answering this question, let's first look at how Equation 1 relates to a realistic problem.

In Equation 1, χ is an interspersed signal. In general, a given signal cannot be expressed directly as an interspersed signal. It is assumed that a given signal is equal to the following equation.

In other words, the signal f we can observe is related to the signal χ vector scattered through a known matrix (or transformation) called Ψ. If signal f is a signal in the time domain, then Ψ is the inverse Fourier transform, and χ vector is the coefficient vector in the frequency domain.

Suppose that the length of this signal is N. The interspersed signal is a K-pulse signal. It is assumed that the signal can be obtained from the signal as follows rather than directly obtaining the signal f .

Where Φ is the M × N measurement matrix. In other words, it is assumed that a signal having a shorter scene length than that of the signal f can be obtained. For audio or video signals, Φ can be thought of as a device that acts as an analog-to-digital converter (ADC). However, in the case where M is much smaller than N, it can be regarded as taking a very small number of samples (i.e., M) than the number of samples required (i.e., N) in a general manner. Now, the relationship between the acquired signal y and the scattered signal χ is

Where A = ΦΨ. That is, it can be seen that Equation 4 becomes as Equation 1. Now we can see that Equation 1 is very important. The nature of the magnitude M × N sensing matrix A plays an important role in estimating the scattered signal χ from the acquired signal y . In this regard, to see the properties of the sensing matrix, constants such as Equation 5 may be defined.

Where δ _k is the smallest constant that satisfies the above inequality for all K-pulse signals χ. These constants are called isometry constants.

When δ _k approaches zero, the left and right terms of the two inequalities approach each other. In this case, the sensing matrix A is a very good sensing matrix with little loss of information about the signal χ. In particular, if M is small, this property can be considered as a very good sensing matrix. If matrix A satisfies Equation 5 for some positive integer K, it is said to satisfy the K-order-limited Isometry Property (RIP). For the sensing matrix A can be obtained as follows.

In summary, in an optimization problem

If A satisfies RIP with δ _2k <√2-1, then this solution satisfies the equation

Where χ _k is a vector of all zero K max values of χ.

According to the above theorem, if χ is K-pulse, the optimization of Equation 6 gives the correct solution. This result is an important result of CS. This theorem shows that the way to solve Equation 1 is the optimization of Equation 6, and the conditions for obtaining the solution can be seen together.

In practice, it is important to obtain a sensing matrix A having good RIP properties. This problem is equivalent to finding a good survey matrix Φ. It is necessary to determine such a matrix well. If A is random, it is known that it can have good RIP properties, and it is important to design the survey matrix Φ to have a random sensing matrix A.

For example, the signal can be sampled at random times. It can be seen that the signal can be reconstructed even with a much lower ratio than the actual Nyquist ratio.

In general, when χ is estimated with respect to the random sensing matrix A as shown in Equation 6, it is known that a scattered signal can be obtained with a very high probability if the condition of Equation 7 is satisfied.

It can be seen that the sample number, M, of a signal satisfying Equation 7 is not large if K is small enough even if N is very large. In addition, it can be seen that the scattered signal does not have much meaning when determining the sample ratio.

Intuitively, we need to look at the l p-norm, which determines the magnitude of the vector, to see how Equation 6 can estimate the scattered signal. The l p-norm of the vector χ of arbitrary length is defined by the following equation.

To write as a real norm, p≥1 must be expanded to p≥0. The shape of the circle according to each norm is shown in FIG. 1 in two dimensions.

 That is, when the value of the above formula is 1, it has a different shape as shown in FIG.

As shown in Equation 6, the shape of the cost function of l 1-norm is as shown in Fig. 1 (a). Commonly used l 2-norm is a circle. In general, when p≤1, the shape of the circle has a more pointed shape on the axis. Now consider the problem of finding a solution that minimizes norm while satisfying a linear constraint such as Equation 6. If the linear constraint equation is shown as a straight line, it can be seen as shown in FIG. 2 (a) shows p = 1, FIG. 2 (b) shows p = 2, FIG. 2 (c) shows p = ∞ and FIG. 2 (d) shows p = 1/2.

In each figure, χ ^ represents the solution for a given norm. For Figures 3 (a) and (d), we get a scattered solution (χ ^ = 0), whereas Figures 2 (b) and (c) have solutions. This is not the case with scattered signals. This example shows that if the scattered vector or signal is a solution that satisfies the linear constraint, the value of p must be 1 or less.

Therefore, l 1-norm optimization plays a very important role in solving the scattered signal in CS. We need an algorithm that efficiently solves the problem of equation (6): Basis Pursuit, Homotopy algorithm, Least Absolute Shrinkage and Selection Operator (LASSO) regression, minimum angle Various methods such as (Least Angle Regression: LARS) regression model have been proposed, and Orthogonal Matching Pursuit (OMP) is one of them.

Looking at the basic radar operation, the pulse-radar system uses the same antenna for both time division and transmission. The transmission signal is transmitted in a specific direction as a short pulse modulated at a specific frequency. In this case, it is assumed that the antenna is very directional. After transmission, the antenna switches to receive mode and receives reflected waves, if any. If there is no reflected wave, it can be considered that there is no target object in this particular direction. If there is a reflected wave, this is regarded as the signal returned by the transmitted pulse signal reflected on the target. Of course, it can be a reflection file reflected by the surrounding feature, not the target, but it can be sufficiently distinguished from the reflected wave by the target since it is a scattered signal that can be recognized by the previous radar operation.

If we can get the time that the transmitted signal has been reflected on the target, we can estimate the distance between the target and the antenna. If the time delay between a transmitting wave and a receiving wave is T, the distance r is as follows.

Where c is the propagation speed. And dividing by 2 is because T is the round trip time. This process proceeds sequentially in the direction in which the target may be, and keeps looking over and over again. Therefore, radar performance is greatly affected by the transmission pulse signal and time delay estimation. In brief, the application of CS to such a radar system, a transmission pulse signal of the radar is called sT (t), when the signal is reflected back to the object, the received signal can be written as follows.

Here u (t) can be thought of as a unique function of an object and can be thought of as an impulse response of an object. In this case, the transmission signal is an input signal and the reflected signal is an output signal. Therefore, if u (t) is obtained, it is possible to know the characteristics of the object, that is, its position, velocity, and shape of the object, so it is important to estimate this function.

In order to estimate a function u (t) having an inherent characteristic of the object, a transmission signal may be designed to have a characteristic as shown in Equation (8).

Signals having this property can be made into pseudo random or pseudo-noise (PN) signals. Given a random pseudo noise sequence, the pseudo noise correlation function is similar to the Dirac-delta function, and assuming that the transmitted signal of this property is matched with the received signal, , The output signal is expressed by the following equation (9).

The output signal of the matched filter is then an approximation of the intrinsic function of the object. Analog matched filters are difficult to implement if the transmit pulse signal is complex. In fact, in the past, this problem was limited by the transmission pulse signal. The development of the digital signal processor can solve the problem by replacing the analog match filter with the digital match filter. That is, it is now possible to use a transmission pulse signal having excellent performance. In this case, u (t) sampled with a digital matching filter can be obtained after first sampling the signal with a high-speed analog-to-digital converter (ADC).

Ultimately, it is important to design a good receiving system using radar signal processing to obtain u (t), which can be used to implement CS with lower complexity. For this purpose u (t) must of course be an interspersed signal. This signal need not be scattered in the time domain. This signal needs to be an interspersed signal through some specific transformation.

However, assuming that it is scattered in the time domain for convenience, u (t) can be written as Equation 10 below.

Where a _k is a coefficient of an arbitrary function of the delay time τ _k , and this function is a K-pulse signal if we sample at a high sampling rate. In this case, the received signal is represented by Equation 11 below.

Suppose that s _T (t) is created using a PN sequence called p. For convenience, the signal may be expressed as shown in Equation 12 below.

Where is a PN signal sampled with a p _n n delay time. A is a K-pulse signal, and the τ _k th term has a nonzero value, that is, a _k and all other values are 0. here

[p0…] By writing pn-1] = Ψ, we can see that the received signal s _R is equal to f. Since M is made of PN signals, if M samples are randomly obtained from s _R (t), the obtained signal can be written as shown in Equation (13).

Where Ψ- is a matrix of row vectors corresponding to M randomly sampled times in Ψ. Such a matrix may have random characteristics according to the PN sequence p and may have good RIP properties. Therefore, the desired K-pulse signal a can be estimated with only a few samples.

In addition, CS can be used in other ways in radar systems. For example, distributed transmit antennas transmit different narrowband signals to detect different targets. If the number of targets is small enough, the received signal is scattered in angle-range-Doppler space. Therefore, even with a small number of samples, the targets can be represented as scattered signals in the angular range Doppler space using CS and estimated.

Seeing targets as signals interspersed in two-dimensional time-frequency space, and looking at and looking at the target in two-dimensional space, like a high-resolution image, for a case where the number of targets is smaller than the number of points in discrete two-dimensional time-frequency space The CS method has already been proposed. These studies show that CS has similar characteristics to those applied for image reconstruction. The other, however, is that the multidimensional signal space has been replaced by the time-frequency or angular range Doppler space, resulting in information needed for the radar, such as the position and velocity of the target.

As mentioned briefly above, among the various algorithms proposed to recover the original signal x from the compressed sensed signal, OMP is one of the representative means. One of the most important problems of OMP is the inverse of the matrix in the least squares process. The problem is that the inverse matrix repeatedly introduces a high degree of computational complexity.

To solve this problem, QR decomposition or Cholesky-basis factorization has been proposed to improve the computational efficiency of the inverse matrix. However, in spite of these methods, the calculation of inverse matrix still has some inefficiencies due to complexity.

On the other hand, OMP using partitioned inverse transform or partitioned inversion has not been proposed because it is recognized that general partitioned inverse results in low computational efficiency.

Republic of Korea Patent Application No. 10-2014-0111141 Republic of Korea Patent Registration No. 10-1703773 Republic of Korea Patent Registration No. 10-1921095

1. Candes, E .; Wakin, M. An introduction to compressive sampling. IEEE Signal Process. Mag. 2008, 25, 21-30. 2. Pati, Y .; Rezaiifar, R .; Krishnaprasad, P. Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. In Proceedings of the IEEE Record of The Twenty-Seventh Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, USA, 1-3 November 1993; pp. 40-44. 3. Tropp, J .; Gilbert, A. Signal recovery from random measurements via orthogonal matching pursuit.IEEE Trans. Inf. Theory 2007, 53, 4655-4666. 4. Chen, S .; Donoho, D .; Saunders, M. Atomic decomposition by basis pursuit. SIAM Rev. 2001, 43, 129-159. 5. Khan, I .; Singh, M .; Singh, D. Compressive sensing-based sparsity adaptive channel estimation for 5 Gmassive MIMO systems. Appl. Sci. 2018, 8, 2076 ~ 3417. 6. Sahoo, S .; Makur, A. Signal recovery from random measurements via extended orthogonal matching pursuit.IEEE Trans. Signal Process. 2015, 63, 2572-2581. 7. Rabah, H .; Amira, H .; Mohanty, B .; Almaadeed, S .; Meher, P. FPGA implementation of orthogonal matching pursuit for compressive sensing reconstruction. IEEE Trans. Very Large Scale Integr. Syst. 2015, 23, 2209-2220. 8. Bai, L .; Maechler, P .; Muehlberghuber, M .; Kaeslin, H. High-speed compressed sensing reconstruction on FPGA using OMP and AMP. In Proceedings of the 19th IEEE International Conference on Electronics, Circuits, and Systems (ICECS), Seville, Spain, 9-12 December 2012; pp. 53-56. 9. Blache, P .; Rabah, H .; Amira, A. High level prototyping and FPGA implementation of the orthogonalmatching pursuit algorithm. In Proceedings of the 11th IEEE International Conference on Information Science, Signal Processing and their Applications (ISSPA), Montreal, QC, Canada, 2-5 July 2012; pp. 1336-1340 10. Tropp, J .; Gilbert, A .; Strauss, M. Algorithms for simultaneous sparse approximation. Part 1: Greedy pursuit.Signal Process. 2006, 86, 572-588. 11. Determe, J .; Louveaux, J .; Jacques, L .; Horlin, F. On the exact recovery condition of simultaneous orthogonal matching pursuit. IEEE Signal Process. Lett. 2016, 23, 164-168. 12. Quan, Y .; Li, Y .; Gao, X .; Xing, M. FPGA implementation of realtime compressive sensing with partial Fourierdictionary. Int. J. Antennas Propag. 2016,2016, 1671--1687. 13. Septimus, A .; Steinberg, R. Compressive sampling hardware reconstruction. In Proceedings of the IEEE International Symposium on Circuits and Systems (ISCAS), Paris, France, 30 May-2 June 2010; pp. 3316 ~ 3319 14. Compressive Sensing Radar Research Trends: Jin Ho Choi, The Magazine of the IEIE, 2014.6, pp 498 ~ 506 15. Compression Sensing and Broadband Radar Signal Acquisition Techniques: Park, Jae-Hyuk Jang, Heung-No Lee, Electronics Journal of the Institute of Electronics Engineers 2016.10, pp 809 ~ 814

The present invention uses the OMP method, which is one of various algorithms proposed to recover the original signal x from the compressed sensing signal as described above, and is a least squares process that is recognized as one of the most important problems of OMP. It aims to provide a new way to reduce the computational complexity in obtaining inverse matrix in.

It is an object of the present invention to provide an OMP method that can increase the efficiency of compression sensing, which can increase the detection performance and resolution compared to hardware performance, especially when operating a radar, which is a representative anti-vehicle device.

An object of the present invention is to provide an OMP method that uses a matrix partitioning inverse different from the known partitioning while reducing the complexity of the calculation when performing the inverse matrix calculation of the least-squares process, particularly in using the OPM method.

In the OMP method for the compression sensing for radar of the present invention for achieving the above object, in the OMP method for processing the received signal of the radar by compression sensing, performing an OMP algorithm for compression sensing, the least square method (least square) The partitioned inversion method (partitioned inversion) is used to calculate the inverse matrix for application.

In the present invention, the partitioned inversion method performs partitioning by taking a partial matrix of an input matrix with respect to the inverse matrix, and in the previous step, when i is an integer, the i-th inverse matrix is calculated by taking the i-th order submatrix, and then i + In the next step of the first step, the process of using the result of the i-th order inverse matrix obtained in the previous step is performed by taking an i + 1 th order matrix including the i th order partial matrix taken in the previous step and performing the i + 1 order inverse matrix. It is made including.

This process can be done repeatedly in ascending order until a final inverse of all the inverses is obtained.

For example, the partitioning inverse matrix method takes a partial matrix from the upper left corner of the matrix rather than splitting it into regions that do not overlap each other in partitioning the input matrix, and then adds rows and columns to the first partial matrix of the input matrix. The addition can be done in ways that take incremental and overlapping partial matrices.

In this case, in calculating the next inverse matrix, it is possible to save the calculation time for obtaining the inverse matrix by using the one obtained in the previous step as it is as the first inverse matrix.

In the present invention, three properties of the input matrix with respect to the inverse matrix, that is, conjugated symmetry, positive definiteness, and overlapped regions are used in the process of deriving the inverse matrix by iterative calculation. On the premise that

According to the present invention, the radar resolution and accuracy are improved even when using the existing radar equipment by using a new OMP method that can reduce the computational complexity in obtaining the inverse matrix in the least squares process. This can increase radar performance.

According to the present invention, it is possible to increase radar performance without the cost and time-consuming elements such as the introduction and replacement of additional expensive equipment through the application of the operating program in the radar system for operation in the existing radar equipment operation.

1 is a view showing shapes of a cost function of l 1-norm that satisfies Equation 6 associated with a linear constraint equation of the prior art;
Fig. 2 is graph shapes when the linear limiting equation mentioned in the description of Fig. 1 is a straight line, in which Fig. 2 (a) is p = 1, Fig. 2 (b) is p = 2, and Fig. 2 (c) is p =. 2 (d) is a graph showing the case where p = 1/2,
FIG. 3 is a table showing, as members, corresponding region overlap states of the partial matrixes of the input matrix for executing the partitioned inverse method in the present invention; FIG.
4 is a table showing comparisons between execution steps of the conventional partitioned inverse method and the partitioned inverse method of the present invention;
5 is a conceptual diagram showing a concept of simplifying inverse matrix calculation using an inverse matrix obtained from a previous iteration step in the present invention;
6 is a graph showing a comparison of the sum of the number of operation steps in the partitioned inverse method of the present invention and the conventional methods as m of the m-th symmetric matrix is increased;
FIG. 7 is a table showing the computational complexity in the partitioned inverse method of the present invention and the conventional methods for the m-th order symmetric matrix U and the single supporter system. FIG.
8 is a configuration and flow conceptual diagram illustrating an example of an SOMP structure for performing a partitioned inverse method of the present invention;
9 is a table showing a comparison of hardware utilization according to the conventional method and the method of the present invention in the experiment of the present invention,
10 is a graph showing a comparison between the restored state of the original input signal to MALTAB and HLS,
11 is a performance comparison table for the 38 * 128 measurement matrix size.

Hereinafter, the present invention will be described in more detail with reference to the accompanying drawings.

The present invention uses three properties of the input matrix for the inverse of each OMP iteration: conjugated symmetry, positive definiteness, and overlapped regions. These properties have been found to reduce computational complexity by reusing the computational results obtained from previous OMP iterations. In terms of comparison of computational complexity, the new partitioned inverse method can improve the complexity over the conventional Cholesky-based inverse method and the conventional partitioned inverse method based on the derived equations.

In addition, multiple measurement vectors (MMV) are applied to the OMP algorithm to improve sparse signal recovery compared to the single measurement vector (SMV) method, which is called simultaneous OMP (SOMP). The inventors implemented SOMP and measured execution time in the field-programmable gate array (FPGA) using the matrix inverse transformation scheme (divisional inverse matrix method) proposed in the present invention. The experimental results show that the total hardware usage is significantly reduced compared to the conventional partitioned inverse method, and the restoration or reconstruction time was measured to 27 μs.

The following describes the condition of the input matrix in the least-squares (LS) problem, together with an overview of the SOMP algorithm, and describes the new partitioned inverse matrix method and experimental results in more detail.

If we first describe the SOMP algorithm,

In SOMP, a linear equation is defined as in Equation 14 below.

[Equation 14]

y = Fx,

Where y ∈ R ^{M × L} , F ∈ R ^{M × N} , and x ∈ R ^{N × L.}

The SOMP process given in Algorithm 1 consists of optimization problems (1) and (2) and least squares (LS) problems (3) and (4). In this process, the residue r ⁰ (i = 1) is initially set to y. During the i th iteration, the optimization problem selects one of the columns of F, which is strongly correlated to the residue of y and retrieves the position k of this column. F _k is a partial matrix including columns according to k, and F ⁱ is updated by summing F _k with the previous sub matrix. The least-squares problem removes the contributed column for the new estimate xtild and computes the new r remaining. Finally, when m iteration is achieved, the final estimate of the original signal is calculated.

In this case, Algorithm 1 inputs the measurement matrix F ∈ R ^{M × N} as the input value and y ∈ R ^{M × L} as the multiple measurement vector, and obtains the output value x tilde ∈ R ^{N × L} as an estimate of the original signal. The variable is the sparseity level m of the original signal x and the remaining r ∈ R ^{M × N} , with initial values r ⁰ = y and xtilde = 0. For the i th iteration, equation (1) k is the index number and j is in the range of 1 or more and N or less k ⁱ = argmaxj∥ <r ^i-1 , F _j > ∥ ₂ , (2) F = [Fk ^{i -1} Fk ⁱ ], (3) x tilde ⁱ = argmim∥ <yF ⁱ x> ∥ ₂ , (4) r ⁱ = yF ⁱ x tilde ⁱ We then repeat the process until i = m to derive the final estimate of x using.

Next, describe the least-squares problem.

In the SOMP algorithm, the recovered signal x ~ (xtilt) gradually becomes similar to the input signal x, and the equation for the new residual in the LS problem is given by the following equation (15).

[Equation 15]

r _i = x ⁱ yF = yy tilde tilde ^{i = yF i (U i)} - 1 (F i) H y

In Equation 15, a reduction in the computational complexity of the inverse of the input matrix U becomes a problem of the OMP algorithm.

Next, look at the conditions of the input matrix in the least-squares problem.

The input matrix U suggests various conditions as follows. U is always symmetrical, as shown in Equation 16 below, and may be expressed as in Equation 17 by having positive limitations.

[Equation 16]

U ^T = (F ^T F) ^T = F ^T (F ^T ) ^T = F ^T F, where T represents the transpose matrix.

[Equation 17]

x ^T (F ^T F) x = (Fx) ^T (Fx) = (Fx) ² > 0, Fx is not 0 unless x is zero.

Equation 17 shows that the eigenvalue Fx is always greater than zero.

That is, U = F ^T F has (1) an inverse of U, (2) the diagonal element C _tt of the matrix is greater than 0 for all t, and (3) the main submatrix also has positive determinism. Satisfies the positive confinement condition by following three characteristics.

In Algorithm 1 (2), F ⁱ is added one column per loop (regression) and the next input matrix U is determined. As the loop proceeds, the next U ⁱ contains the U ^{i- 1} of the previous loop as shown in FIG. For example, U ² is the extension matrix of U ¹ by the previous iteration, so U ⁱ can be divided into a previously obtained portion and a newly added portion.

Hereinafter, a look at the partitioned inverse method proposed by the present invention using the characteristics of the input matrix described above.

First, the conventional inverse matrix method is a method of obtaining an inverse matrix by dividing an input matrix into four parts. Because of the low number of conditions, this inverse method is noise-resistant and works by using the lower part of the input matrix in symmetric and positive limiting (SPD) characteristics. The steps for the inverse of the conventional method are shown in the left part of FIG. However, as the number of steps increases and the matrix size increases, the number of calculations increases significantly.

On the other hand, in the partitioned inverse matrix method of the present invention, in order to reduce the complexity of calculating the inverse matrix, a relational expression as shown in Equation 18 is obtained by using the above-described input matrix property.

Equation 18

A = A ^T → A ^-1 = (A ^{- 1} ) ^T , D = D ^T → D- ¹ = (D ^{- 1} ) ^T , B = C ^T

U as the input matrix is a symmetric matrix and determines whether or not the constituent matrix (A, B, C, D) of the input matrix satisfies Equation 18 to be used in the next optimization process.

Step 1 in the conventional partitioned inverse method of Figure 4 can be eliminated. That is, in the i th OMP iteration, U ⁱ is divided into four partial matrices containing U ^{i −1} as shown in FIG. 5. Unlike the conventional partitioned inverse method, the steps of the method include the inverse U ^{i -1} that has already been used in the previous OMP iteration. Therefore, the conventional first step is unnecessary because it uses (U ^i-1 ) ⁻¹ shown in gray in FIG.

Conventional steps 2 and 3 may be reduced to one step based on Equation 18 as follows.

[Equation 19]

(A ^-1 B) ^T = (A ^{- 1} C ^T ) ^T = C (A ^-1 ) ^T = CA ^-1

Conventional steps 7 and 8 may also be simplified by Equation 18 as follows.

[Equation 20]

F = A ^-1 B (D-CA ^-1 B) ^-1

[Equation 21]

G = (- (D-CA -1 B) - 1 CA - 1) T = (A - 1) T C T ((D-CA -1 B) -1) = -A -1 B (D-CA ^-1 B) ^-1

Equation 18-As shown in Fig. 2 from Equation 21, the number of steps of the inverse matrix method decreases from 10 to 7.

Looking at the computational complexity for the partitioned inverse method proposed in the present invention,

6 shows a comparison of the complexity between the three inverse matrix methods. The choresky-based method and the conventional partitioned inverse method continuously calculate the inverse of the matrix with increased size per iteration, whereas the inverse matrix method of the present invention calculates the remaining partial matrix except for the iterative calculation of the previously obtained part. . Therefore, the proposed method improves the computational complexity problems of multiplication, addition / subtraction, and division according to the input matrix size. In this case, the matrix size may increase from 1 × 1 to m × m per loop. To visualize the improved complexity, MATLAB (2013b, The MathWorks, Natick, MA, 2007), to calculate and compare the total number of steps according to the matrix size of the choresky-based inverse matrix method, the existing partitioned inverse method, and the proposed partitioned inverse method, USA) simulation was performed to obtain the same results as the graph of FIG. In FIG. 6, m denotes the size (m × m) of the matrix.

Looking at the partitioning inverse for a multi-support system, the computational complexity features as shown in the table of FIG. 7 apply only to a single supporter system. However, the inverse matrix method of the present invention can be used in a multi-support system, and the complexity equation of FIG. 7 changes from the single system to the multiple system, and thus the following Equations 22 to 8 are applied to the multiplication operation, the addition subtraction operation, and the division operation. Changes to 24

[Equation 22]

(2m ² -1) n-2mn ² + n ³

[Equation 23]

(2m ² -2m + 1) n-2mn ² + n ³

[Equation 24]

n

In this case, n is the number of the supporters, and as n increases, the computational complexity increases, but the inverse matrix method proposed in the present invention can be applied to a multi-support system.

Referring to the SOMP structure involving the partitioned inverse method, the SOMP structure with the partitioned inverse method is configured as shown in FIG. Referring to Figure 8, the correlation matching block (CRMB) is looking at the position k of the associated column or columns of the measurement matrix F, F ⁱ is the summation of in the previous iteration stored by the memory block F ^{i -1} and F ⁱ Fk Is updated by U block (GUB) generation separates the rows of the measurement matrix corresponding to k and generates the matrix U by the inner product. The inverse of matrix U is generated by a U block inverse transform (IUB) stored in a U- ¹ memory block, which provides the inverse of U calculated in the previous step. Thus, IUBs can reduce computational complexity by using previously computed inverse matrices. The least-squares block (LSB) calculates x-tilde (an estimate of the original signal x), and the residual r ⁱ is generated by the residual block (RB). Finally, the final estimate signal is obtained by the system at m iterations.

Looking at the experiment for a while, we used a field programmable gate array (FPGA) implementation approach in the effect experiment of the present invention.

As a test equipment, we verified the operation of the proposed inverse matrix and the entire SOMP algorithm using the Xilinx Kintex UltraScale board and the XCKU115-FLVA2104-2-I chipset. The board has a large number of digital signal processing (DSP) slices, which is advantageous for the design of large algorithms.

Xilinx's High Level Synthesis (HLS) tools for FPGA design allow (1) data types, math, IP, linear algebra and other library processing with arbitrary precision, and (2) register transfer levels due to the co-simulation provided. (RTL) automatic extraction and easy identification of resistor transistor logic circuits (RTL), and (3) C synthesis provides design resources such as clock, region, and I / O port descriptions. Can be easily changed, and (4) facilitates the use of IP modules in other design tools such as system generators, enabling high-performance hardware design with little hardware knowledge.

Taking advantage of these advantages, the C / C ++ code is converted to a hardware description language (HDL) such as Verilog, and the device is based on 16-bit and 32-bit fixed-point complex operations. Most of the calculations were done with 16 bits, but the split reversal process used 32 bits. This is because the partitioning operation requires a more accurate calculation.

In terms of the parameters used in the experiment, the number of physical channels was four, and the measurement matrices (M × N) were 32 and 128, respectively. At channel reuse factor 8, sampling frequency of 250 MHz, Nyquist bandwidth is 4 gigahertz, MMV extension 16, F _P is 31.25 megahertz, and FFT point is 128.

For further optimization of FPGA implementations, SOMP algorithms require large matrix internals. You need three as loops to represent the inner product of your code. We apply the pipeline second as a loop to perform real-time signal reconstruction or reconstruction. Although the number of DSPs increases in proportion to the innermost number as a loop, the latency can be greatly reduced because the calculations can be parallel.

Most of the operations that occur internally are complex types, which increases the complexity of the operations. Even if complex division is required in the implementation of the partitioned inverse method, the denominator is always a positive real number and can therefore be replaced by two actual partitions.

In terms of SOMP hardware utilization, various components are required for the hardware implementation of the SOMP algorithm. Block RAM effectively stores vectors and matrices. In this case, the measurement vector is y, the measurement matrix is F, and the rest are r. DSPs are multiply-accumulate (MAC) circuits prefabricated in FPGAs and are essential hardware because they quickly perform multiplication and add / subtract operations. Flip flop (FF) is a shift register used to synchronize logic. As can be seen from the table of Figure 9, the overall hardware utilization for the partitioned inverse method of the present invention is reduced compared to the conventional partitioned inverse method.

However, only the use of BRAM for the inverse matrix method of the present invention is increased compared to the conventional method. This is because the U matrix and the U matrix of the previous step must be stored at the same time.

Regarding signal recovery, perform MATLAB simulations and HLS to verify the recovered signal from the input pulse signal, and the data precision is converted from 16 (0.12) to 32 (0.24) bits, where precision (p) and partial bits (factional bits) (f) is used to define the data type p (.f). Fig. 10A shows an input pulse signal mixed with a random carrier signal. The signals are well reconstructed or reconstructed by simulation and HLS as shown in Figs. 10 (b) and 10 (c), respectively.

In the experiment, the peak signal-to-noise ratio (PSNR) was obtained for objective evaluation, and the PSNR can be represented by the following equations (25) and (26).

[Equation 25]

PSNR = 20log ₁₀ (MAX / √MSE)

[Equation 26]

MSE = 1 / N × Σ [x ⁱ -x ~ ⁱ ] ²

Where MAX is the maximum possible value of signal x, N is the total number of samples, x to ⁱ (x til ⁱ ) are the sample values of the reconstructed i-th signal, and x ⁱ is the sample value of the i-th original signal.

The PSNR of the SOMP structure was determined to be 30.26 dB for 16-bit and 32-bit data precision.

The table in FIG. 11 compares the performance of the conventional pulse signal recovery algorithm apparatus. The table focuses on the OMP algorithm rather than the inverse method of other studies, but an objective comparison is needed to emphasize the inverse method of the present invention. For this reason, papers of the same measurement matrix size (32 × 128) were selected while applying different inversion methods (ie, QR decomposition based and Coreski based inverse method). Although higher clock frequencies are used than other studies in the table, the restoration time by the method of the present invention is fast enough when it has greater scattering than other literature data (Non-Patent Documents 7, 12, 13) showing the results of conventional experiments. You can see that.

In conclusion, the SOMP algorithm generally results in high computational complexity for the LS problem because the internal problem is very large. The SMV model also increases this complexity. In order to reduce the complexity, we propose a new SOMP algorithm for the MMV model, but apply a new partitioned inverse method that eliminates the recalculation of the inverse of the input matrix in the previous loop. The inverse matrix method of the present invention restores a signal with a small amount of hardware in comparison with the conventional inverse partition scheme. Using a Xilinx Kintex Ultrascale board and HLS, the signal recovery was verified with high PSNR and the original signal recovery time was reduced. Thus, compared to the conventional partitioning inverse and Coreski-based partitioning, it can be used more appropriately for applications requiring high level restoration such as radar with less hard air usage.

The present invention has been described above by way of limited embodiments, which are only illustratively described to help understanding of the present invention, and the present invention is not limited to these specific embodiments.

Therefore, one of ordinary skill in the art to which the present invention pertains may make various modifications or applications based on the present invention, and such modifications and applications belong to the appended claims.

Claims (3)

delete In the OMP method of processing the received signal of the radar by compression sensing, performing an OMP algorithm for compression sensing,
Use the partitioned inversion method to calculate the inverse for least squares,
In the partitioned inversion method, when performing a partition taking a partial matrix of an input matrix with respect to an inverse matrix, i + 1th order is calculated by taking an i-th order sub-matrix when i is an integer and performing i-th order inverse matrix. The next step of the step includes taking the i + first order matrix containing the i order sub-matrix taken in the previous step and using the result of the i-order inverse matrix obtained in the previous step while calculating the i + 1 order inverse matrix ,
Repeating the process while increasing the order to obtain the final inverse matrix consisting of all inverse matrix parts OMP method for radar compression sensing. The method of claim 2,
The OMP method, when the basic linear equation is given by y = Fx, y ∈ R ^{M × L} , F ∈ R ^{M × N} , x ∈ R ^{N × L} ,
Input the measurement matrix F ∈ R ^{M × N} as the input value and y ∈ R ^{M × L} as the multiple measurement vector, and obtain the output x x til ∈ R ^{N × L} as an estimate of the original signal.
The variable is the spacing level m of the original signal x and the remaining r ∈ R ^{M × N} , initially set to r ⁰ = y, xtilde = 0,
For the i-th iteration, (1) k is the index number and j is in the range of 1 or more and N or less. k ⁱ = argmaxj∥ <r ^i-1 , F _j > ∥ ₂ , (2) F = [Fk ^i-1 Fk ⁱ ], (3) x tilde ⁱ = argmim∥ <yF ⁱ x> ∥ ₂ , (4) r ⁱ = yF ^{i Use} x tilde ⁱ to proceed the process until i = m to derive the final estimate of x. Repeating
OMP method for compression sensing for radar comprising a SOMP process.

Images