Real-time FPGA Design for OMP Targeting 8K Image Reconstruction

Assuming the length of original image signal $X$ is $N*1$. Using a transform matrix marked as $\psi$ with the size of $N*N$, transforming the original signal to the transform domain to gain the sparse signal called $\theta$ with the size of $N*1$. Sparse means, in the signal, there are only a few elements are non-zero elements. The sparsity level is marked as $k$. Projecting the original signal on the measurement matrix $\phi$, then the compressed signal $Y$ called measurements is obtained. The size of measurement matrix is $M*N$, the size of signal $Y$ is $M*1$. The relationship between $k,M,N$ is $K<<M<N$. The $M/N$ represents the sampling rate. Multiplying the measurement matrix with the transform matrix is the sensing matrix $A$. The theory of CS is as Equation.1 shows.

Compressed Sensing has two main principles. One is the sparsity of the original signals. This means the signal is compressible. The other one is incoherence. The measurement matrix with the transform matrix should have incoherence. The incoherence ensures the possibility of reconstructing the original. In other words, the sensing matrix, which is gained by multiplying the measurement matrix with the transform matrix, determines the quality of the reconstruction results. Thus, choosing proper sensing matrix is important.

Orthogonal Matching Pursuit (OMP) is a reconstruction algorithm that belongs to the greedy algorithm. To reconstruct the original image signal, OMP first uses Dot Product to find the index of the column vector of the sensing matrix with the highest contribution, then puts this index into the index set. After that, the column vectors in the index set are used to approximate the original signal and update the residual. Both the procedure of approximating the original signal and the procedure of updating the residual include the Least Square Problem. The procedure of OMP reconstruction is as follows:

a. Initialize the residual $r_{0}=Y$, the iteration counter $t=1$, the index set $S_{0}=\emptyset$.

b. Find the index of the highest contribution column vector of the sensing matrix. $Index_{t}$ = $argmax|A^{T}*r_{t-1}|$.

c. Update the index set $S_{t}=S_{t-1}\cup Index_{t}$.

d. Use Least Square Problem to approximate the original signal. $\hat{\theta}=(A_{S}^{T}*A_{S})^{-1}A_{S}^{T}*Y$

e. Judge if satisfy the quit condition or not. If yes, go to g. If not, go to f.

f. Update residual. $r_{t}=Y-A_{S}*\hat{\theta}$, then go back to b.

g. Inverse the approximate signal to obtain the reconstructed image signal. $\hat{X}=\psi_{N*N}*\hat{\theta}$.

OMP is an iteration-based reconstruction algorithm. If the reconstruction procedure satisfies the quit condition, the approximation signal will multiply with the transform matrix to obtain the reconstructed image signal. The iteration time is usually regarded as the sparsity level of the signal in the transform domain. But the sparsity is unknown for practical image reconstruction. If the iteration time is too large, the reconstruction time will be considerable. If the iteration time is too less, the quality of the result will be worse. Thus, the strategy of the quit condition is crucial in the design.

In this work, the Hadamard matrix is adopted as the measurement matrix, the Fast Walsh-Hadamard transform matrix is chosen as the transform matrix. The sensing matrix $A$, which is used in the reconstruction procedure, is as the Equation.2 and Fig.1 shows.

In the Equation.2, the value of $V$ changes with the different choices of M and N. Especially, the range of $i$ is $0<i\leq M$, the range of $j$ is $0<j\leq N$. It is easy to find that when $j>M$, $A(i,j)$ is zero. That is, if the vector index is bigger than $M$, the column vectors are all zero vectors. Besides, the non-zero column vectors of the chosen sensing matrix are sparse vectors, which contain zero elements mainly. Based on these features, the Dot Product and Least Square Problem are optimized further to reduce the time cost to achieve the aim of real-time reconstruction.

This subsection introduces the optimized Dot Product implementation. The Dot Product is used to find the index of the highest contribution column vector of the sensing matrix. The original calculation is illustrated in Equation.3.

Because the column vectors with index bigger than $M$ are all zero vectors, the vector-based multiplication results of these vectors are zero. Hence, the Dot Product in Equation.3 is optimized as Equation.4 shows.

In addition, the column vectors from 1 to $M$ are sparse vectors, the calculation in Equation.3 can be revised further. Deleting the calculations of zero elements, the Dot Product calculation is as Equation.5 shows.

After obtaining the result of the Dot Product, this result will be sort to find the index of the highest contribution column vector of the sensing matrix. This sort procedure is implemented by 2-to-1 comparison calculation. As Equation.5 presents, there is a common divisor in the result of Dot Product. The multiplication with this common divisor, which is represented as $V$ in Equation.5, can be deleted to save both time and hardware cost.

Because in the first iteration, the residual is equals to the measurements, and the values of measurements are positive numbers, the $result_{DP}(1)$ must be the biggest value. That is, the $Index$ of the first iteration must be 1. In OMP, all of the column vectors only can be chosen once. Based on these conditions, in the first iteration, the procedure in Equation.3 is omitted, the result of $Index$ in the first iteration is set as 1 and the $result_{DP}(1)$ in other iterations is set as 0. The implemented procedure of finding the index of the highest contribution column vector is as Equation.6, Equation.7 present.

Fig.2 shows the diagram of the implemented procedure when the iteration index is bigger than 1. First, the first element of the residual is set as 0, other elements are set as its absolute value. Then using the 2 to 1 comparison to find the biggest value. Last, outputting the index of the biggest value.

This subsection introduces the optimized Least Square Problem (LSP). The LSP is the key point at both residual update in the iteration and the original signal approximation after quitting the iteration. The calculation of LSP is as Equation.8 shows.

As mentioned in Section 3.1, the $A_{S}$ contains sparse column vectors. Hence, in the vector-based multiplication of $A_{S}^{T}*A_{S}$, only several situations, which are illustrated in Equation.9, get non-zero results.

Marking the $A_{S}^{T}*A_{S}$ as $Result_{matrix}$, the Equation.10 shows the equation of $Result_{matrix}$.

The range of $i$ is $1\leq i\leq index_{iteration}$, the range of $j$ is $1\leq j\leq index_{iteration}$. The $index_{iteration}$ represents the iteration index. The value of $P$ and $Q$ depends on the iteration index.

Because the value of matrix $Result_{matrix}$ is certain in each iteration, the value of $Result_{matrix}^{-1}$ is certain too. Hence, the complicated inverse calculation is substituted as using a Look Up Table (LUT) to find the value of $Result_{matrix}^{-1}$ of each iteration. But the inverse value leads to fractional numbers which is hard to implement in hardware, shifting the input is indispensable.

The $A_{S}^{T}*Y$ in Equation.8 is still a time cost procedure. As mentioned in Section 2.2, the index set $S$ unions the max index of the current iteration with the indexes stored in the previous iterations. The $A_{S}^{T}*Y$ also can be replaced by combining the result of previous iteration with the result of the current iteration. So only the new result of $A_{S}^{T}*Y$ is calculated in each iteration. Besides, the column vectors of $A_{S}^{T}$ are sparsity, so the vector-based multiplication can also be substituted by the simple shift and add calculations.

The implementation details of the Least Square Problem is as Fig.3 shows. First, multiply measurements with the transpose of sensing matrix. Then multiply it with the inverse matrix, which is stored as LUTs, to obtain the approximate signals in transform domain. Next, multiply the approximate signals in transform domain with the sensing matrix. After that, use the shift measurements to update the residual. Judge if meets the quit condition or not. If yes, output the shift value of the approximate signals in transform domain. If not, output the shift values of updated residual.

Fig.4 shows the overall architecture of this work. As illustrated in Section 3.2, the first iteration of the reconstruction omits the Dot Product procedure. The result of the Dot Product in the first iteration is defaulted as 1. So the reconstruction procedure of this work starts from LSP to update the residual. In the LSP, if the quit condition is satisfied, LSP will output a flag named $reconstruction\_finish$ and the shift of the approximate signals in transform domain. If the $reconstruction\_finish$ is not equals to 1, the output residual of LSP will put into the Dot Product to find the index of the highest contribution column vector of the sensing matrix. As mentioned in Section 3.3, there are calculations duplicate in every iteration. Thus, the max value of residual and max index of column vector are stored in every iteration. These parameters are inputs of LSP. This loop will be quit until the $reconstruction\_finish$ is 1. Then multiplying the output of the LSP with the transform matrix to obtain the reconstructed result.