Wavelet Compressive Sampling Signal Reconstruction Using Upside-Down Tree Structure

Abstract

This paper suggests an upside-down tree-based orthogonal matching pursuit (UDT-OMP) compressive sampling signal reconstruction method in wavelet domain. An upside-down tree for the wavelet coefficients of signal is constructed, and an improved version of orthogonal matching pursuit is presented. The proposed algorithm reconstructs compressive sampling signal by exploiting the upside-down tree structure of the wavelet coefficients of signal besides its sparsity in wavelet basis. Compared with conventional greedy pursuit algorithms: orthogonal matching pursuit (OMP) and tree-based orthogonal matching pursuit (TOMP), signal-to-noise ratio (SNR) using UDT-OMP is significantly improved.

1. Introduction

Compressive sampling (CS) [1, 2], as a new emerging signal processing theory, which has received considerable attention in signal processing. CS exploits the sparse structure in signal, and it enables signal reconstruction from a small number of random samples. A variety of signal recovery algorithms have been proposed to reconstruct the sparse signal. Generally, there are two classes of methods: convex optimization and greedy pursuit algorithms. Although the convex optimization method, such as basis pursuit (BP) [3], is powerful for sparse signal reconstruction, it may be computationally burdensome. Considering the computational complexity and difficulty of realization, greedy pursuit algorithms, especially matching pursuit (MP) [4] and orthogonal matching pursuit (OMP) [5], are attractive for engineering problems.

MP algorithm is computationally efficient and often features good performance; however, when the basis Ψ is not given by an orthogonal basis, MP algorithm cannot find the best approximation for original signal [6]. As an alternative, OMP orthogonalizes each selected column vector associated with maximum projection in basis Ψ, and OMP does not suffer the aforementioned flaw. In CS field, although there are many works about theory analysis and practice implementation of OMP, these recovery algorithms are generic in the sense that they do not exploit any particular structure in the signal besides its sparsity in some basis. However, for some signals there is additional a priori information that we can exploit for improving recovery performance. For example, the piecewise smooth signals, which are widely used in practice, are not only sparse in wavelet domain, but also form a connected subtree.

In this paper, we present an improved OMP signal recovery algorithm by employing an upside-down tree structure of signal in wavelet domain (we refer to this tree-based algorithm as UDT-OMP). Proposed algorithm is evaluated by signal-to-noise ratio (SNR) as a measure of quality of reconstructed signal. We have compared the performance of UDT-OMP with OMP and tree-based orthogonal matching pursuit (TOMP) with SNR as a function of number of measurements.

2. Compressive Sampling (CS) Background

CS is a novel sampling paradigm that goes against the common wisdom in data acquisition. CS states that a sparse or compressible signal can be recovered from a small salient set of random projections. To make it possible, there are two fundamental premises [7]: sparsity, which pertains to the signals of interest, and incoherence, which pertains to the sensing modality. Sparsity expresses the idea that the “information rate [8]” of a high-dimensional signal may be much lower than the maximum frequency presented in the signal. Consider an 𝑁-dimensional signal 𝑥, which can be sparsely represented in an appropriate basis Ψ, and the transfer coefficients of 𝑥 are given by 𝛼_𝑖=<𝑥,𝜓_𝑖>(𝑖=1,2,…,𝑁). There are no more than 𝐾 nonzero entries in 𝜶=[𝛼₁,𝛼₂,…,𝛼_𝑁]^𝑇(𝐾≪𝑁), and the signal 𝑥 is called K-sparse signal. Just the opposite, incoherence means that the measurement matrix Φ has dense representation in the basis Ψ, and Φ is independent of Ψ.

CS also extends to that so-called compressible signals that are not exactly sparse but can be closely approximated as such (i.e., wavelet coefficients of signal and image). Sparse signals have coefficients 𝜶 that, when sorted, decay according to a power law: |𝛼_𝑖|<𝐶⋅𝑖^−1/𝑝 for some positive constant 𝐶 and 𝑝⩽1; the smaller the decay exponent 𝑝, the faster the decay and the better the recovery performance we can expect from CS. In practice, most of man-made and natural signals are sparse or compressible in the sense that they have concise representations when expressed in an appropriate basis, such as Fourier basis and wavelet basis. For a K-sparse signal 𝐱, we can find its 𝑀(𝑀=𝑂(𝐾log(𝑁/𝐾)) [9]) linear measurements 𝐲=Φ𝐱(𝐱=Ψ𝜶), and Φ is the measurement matrix of size 𝑀×𝑁(𝑀<𝑁), and signal 𝐱 can be reconstructed by solving the following inverse problem:min‖𝜶‖₀s.t.𝐲=𝚽𝐱=𝚽𝚿𝜶, (2.1) where the 𝑙₀ norm used here simply counts the number of nonzero entries in 𝜶.

Since 𝑀<𝑁, (2.1) is an ill-posed problem, and there are many possible solutions of (2.1). The original signal 𝐱 can be reconstructed from 𝐲 by exploring its sparse expression, that is, among all possible 𝜶 that satisfy 𝐲=ΦΨ𝜶, seek the sparsest. It is known to be NP-hard to solve problem (2.1) and different suboptimal strategies are adopted in practice such as BP and OMP.

3. Orthogonal Matching Pursuit (OMP) Algorithm

MP approximation is improved by orthogonalizing the directions of projection with a Gram-Schmidt procedure; it is known as orthogonal matching pursuit (OMP). This orthogonalization was introduced in [6], and it was intensively studied in [5]. The improved precision of this orthogonalization is penalized by the higher computation of the Gram-Schmidt. In CS, OMP selects column vectors from dictionary 𝐃(=ΦΨ) that minimizes the difference, also called the residual 𝐫_𝑘, between measurement 𝑦 and approximation. Specially, starting with 𝐫₀ = y, the OMP algorithm searches for the kth atom with maximum projection as𝜆_𝑘=argmax_𝑖𝐫_𝑘−1,𝐝_𝑖, (3.1) and it updates the residual as𝐫_𝑘=𝐲−𝑃_{span{𝐝𝜆1,𝐝𝜆2,…,𝐝𝜆𝑘}}𝐲. (3.2) Here, 𝐝_𝑖 is the ith column vector of D, and 𝑃_𝑉 denotes the orthogonal projection onto the subspace 𝑉, and 𝑃_𝑉=𝑉(𝑉^𝑇∗𝑉)−1∗𝑉^𝑇. OMP algorithm in parallel applies the Gram-Schmidt orthogonalization upon chosen atoms for efficient computation of projections. It can be formulated as a subset selection problem where a minimum subset of columns of D matrix is chosen to approximate the observation vector 𝐲 in the least square sense. The OMP algorithm successively chooses an additional column of D matrix to reduce the approximation error. Equivalent, the OMP method begins with a tentative solution of 𝜶 with a single nonzero entry, and gradually adds nonzero entries one by one until the approximation error of 𝐲 meets a pre-determined criterion. The accuracy of approximation increases with the number of iterations; however, the iterative number must not be bigger than the number of random measurements and we cannot expect to recover signal from higher dimensions than that of measurements.

4. UDT-OMP Algorithm

4.1. Tree Structure in Wavelet Domain

In CS, both BP and OMP reconstruct signal based on sparsity without considering any other particular structure that may exist in the signal. The wavelet transform of a piecewise smooth signal (many punctured real-world phenomena give rise to such signal), an important subclass of sparse signals, yields a sparse, structured representation of signals in this class: the significant coefficients tend to form a connected subtree of the wavelet coefficient tree. Figure 1 shows an example of wavelet representation of such signal. Only a few wavelet coefficients are significant, which form sparse subtrees when represented in multiscale wavelet transform. The number of such subtrees equals to that of the discontinuities in signals.

(a)

(b)

(c)

(d)

Figure 1 

Example of sparse representations in the wavelet domain: (a) piecewise smooth signal; (b) wavelet transform of piecewise smooth signal, only few coefficients are significant and they are compressible or near sparse; (c) the significant coefficients are well organized in tree structure across the scales; (d) binary tree for piecewise smooth signal, the significant wavelet coefficients arise from the discontinuities in the signal (the black circles denote the large wavelet coefficients).

In this work, we only focus on 1D signals, and similar arguments apply for 2D and multidimensional signals. Consider a signal 𝐱 of length 𝑁=2^𝐽, given a bandpass wavelet function 𝜓(𝑡) and a lowpass scaling function 𝜙(𝑡), the discrete wavelet transform (DWT) coefficients 𝜶 of size 𝑁 in terms of shifted versions of 𝜙(𝑡) and shifted and dilated versions of 𝜓(𝑡). The wavelet representation of 𝐱 is given [10]𝐱=^𝑁𝐽−1_𝑖=0𝑣_𝐽,𝑖𝜙_𝐽,𝑖+^𝐽𝑁_𝑗=1^𝐽−1_𝑖=0𝑤_𝑗,𝑖𝜓_𝑗,𝑖, (4.1) where 𝐽 denotes the scale of analysis and scale 𝐽 indicates the coarsest scale or lowest resolution of analysis. 𝑁_𝑗=𝑁/2^𝑗 is the number of coefficients at scale 𝑗∈{1,…,𝐽}, and 𝑖 is the position, 0≤𝑖≤2^𝑗−1. In terms of matrix notation, 𝐱 has the representation 𝐱=Ψ𝜶, where Ψ is a matrix containing the scaling and wavelet functions as columns and the scaling and wavelet coefficients are as follows:𝜶=𝑣_𝐽,𝑖𝑤_0,0𝑤_1,0𝑤_2,0…^𝑇. (4.2)

According to the statistical analysis in wavelet domain, wavelet coefficients have the following two properties [11].Compression. The wavelet transforms of real-world signals decay exponentially as the scale become finer, and they tend to be sparse (as depicted in Figure 1(b)).Tree structure. Those significant wavelet coefficients propagate across scale in wavelet tree, and they are well organized in a tree structure (as depicted in Figure 1(c)).

This tree structure was exploited by previous CS reconstruction algorithms known as iterative reweighted 𝑙₁-norm minimization with the wavelet Hidden Markov tree model (HMT-IRWL1) [12] and TOMP [13]. HMT-IRWL1 integrates the HMT model to enforce the wavelet coefficient structure during IRWL1, which updates the weight values with state probability of HMT model and highly depends upon accuracy of model. If accurate Markov model is available, then the HMT-IRWL1 could be powerful to recover the sparse signal. Although, in practice CS applications, only a small set of random measurements are available, one cannot get an accurate Markov model from these measurements. TOMP evaluates the sums of projections along each wavelet coefficients connected by a subtree, which is an improved OMP recovery algorithm. Both HMT-IRWL1 and TOMP are based on upright subtree, and one father node connects with two children nodes, as depicted in Figure 1(d).

4.2. UDT-OMP Algorithm

In previous tree-based CS reconstruction algorithms, trees are assumed to be upright, which means that those significant coefficients propagate from coarser scale to finer scale. Examining the tree structure in Figure 1(c), we noticed that not only do they organize in subtrees, but also those subtrees are upside down. From finer scale to coarser scale, wavelet coefficients become larger, and there are more significant coefficients in coarser scales. In this paper, we motivate an improved version OMP by employing weighting an upside-down tree (UDT-OMP).

The input of UDT-OMP is a dictionary D with size of 𝑀×𝑁(𝑀≪𝑁), a measurement vector 𝐲 with length of 𝑀, upward extending coefficient 𝑙, weighting value 𝑐, and iterative number 𝐾. UDT-OMP returns a reconstruction sparse vector 𝜶 of length 𝑁 that is subject to 𝐲=ΦΨ𝜶.

UDT-OMP evaluates the projection of each single column vector and searches for the maximum projection. According to the columns associated with the maximum projections, UDT-OMP constructs an upside-down subtree upward coarser scale with depth of 𝑙, which will be weighted in the next iteration of searching for the maximum projection.

In (4.2), we noticed that the wavelet transform of signal 𝜶 can be divided into two parts: scaling coefficients and wavelet coefficients [13]𝐯𝐰𝜶=, (4.3) where 𝐯 contains all the scaling coefficients and 𝐰 contains all the wavelet coefficients of signal 𝐱. Since the scaling coefficients are significant, we aim to recover all the scaling coefficients 𝐯 and those significant wavelet coefficients in 𝐰. In particular, we recover 𝐯 and 𝐰 separately and revise (2.1) as follows:min‖𝐰‖₁𝐯𝐰s.t.𝐲=𝚽𝐱=𝚽𝚿. (4.4)

UDT-OMP algorithm consists of two steps. We first limit our search space in the columns associated with the scaling coefficients and then in the columns associated with the wavelet coefficients. Let Λ_𝑘 be the selected set of columns in dictionary D from the beginning 𝑘 iterations. Let 𝐔_𝑘 be the Gram-Schmidt version of Λ_𝑘. According to (3.1) and (3.2), OMP aims to search for the columns with the maximum projections, which associate with the significant coefficients. Firstly, we set the selected sets Λ₀ and U₀ to be empty. In step 1, we select all columns of D corresponding to the scaling coefficients since these coefficients are significant. According to the construction of wavelet basis Ψ, the wavelet transform level 𝐽 (also called the scale of analysis) is deterministic and so as the positions of scaling coefficients (as depicted in (4.2)):𝚲₀=𝐃allcolumnsassociatedwiththescalingcoefficients. (4.5)

All vectors in Λ₀ are sequentially orthogonalized using Gram-Schmidt and stored in U₀. At end of step 1, the residual is updated by𝐫₀=𝐲−𝑃_span{𝚲0}𝐲. (4.6) In step 2, we focus on recovering the wavelet coefficients which have sparse tree structure. Step 2 is the repetition of iterations. Let {𝐒_𝑘}^𝐾_𝑘=0 be the candidate set of vectors, which will be weighted in the next searching of maximum projection. We initialize all the elements of 𝐒₀ as 1.

In the first iteration, we initialize the counter 𝑘=1. Then UDT-OMP investigates the weighted projections of the current residual on all the column vectors of D𝜽^𝑘_𝑖=𝐒_𝑘−1(𝑖)⋅𝐫_𝑘−1,𝐝_𝑖,𝜆^𝑘=argmax_{1≤𝑖≤𝑁}𝜽^𝑘_𝑖. (4.7)

According to the maximum projection position 𝑖, we construct the kth upside-down subtree 𝑇_𝑘 with depth of 𝑙 (𝑙 is an integer, 𝑙 should be greater than 1 and not greater than the level of wavelet transform). At each finer scale, every four consecutive nodes connect two nodes of nearest coarser scale, which forms an upside-down subtree. For example, if 𝑤_3,2 (as depicted in Figure 1(d)) is current maximum projection, 𝑤_2,0 and 𝑤_2,1 (both connect with four consecutive nodes 𝑤_3,0, 𝑤_3,1, 𝑤_3,2, and 𝑤_3,3) would be children nodes in the nearest coarser scale, and 𝑤_1,0 and 𝑤_1,1 are children nodes in the next coarser scale:𝐒_𝑘(𝑖)=𝑐∶𝑖∈𝑇_𝑘. (4.8) Since we want to enforce the selected columns, the weighting value 𝑐 should be greater than 1. We add the newly chosen node into our selected set Λ_𝑘:𝚲_𝑘=𝚲_𝑘−1∪𝜆^k,𝐔_𝑘=𝑃_GS𝚲_𝑘, (4.9) where 𝑃_GS(·) denotes orthogonalization using Gram-Schmidt. The residual is updated as follows:𝐫_𝑘=𝐲−𝑃_{span{𝚲𝑘}}𝐲. (4.10) After 𝐾 iterations of (4.7)–(4.10), the significant coefficients are determined by the columns in Λ_𝑘, and they are represented as follows:𝜶=𝐔_𝑘𝐲, (4.11) and the approximation of the original signal is𝐱=𝚿𝜶. (4.12)

The pseudocodes of UDT-OMP are described as in Algorithm 1.

Inputs: CS matrix D, measurement 𝐲, upward extending coefficient l, weighting value c, and iterative number K, wavelet transform level J, 𝐒0 = {1} Outputs: approximated signal 𝑥 Initialize: 𝜶0 = {0}, 𝐫0 = {0}, k = 0, Λ0 = {1,2,…, 𝑁𝐽 (= N/2J)} 𝑟0=𝑦−𝑃span{Λ0}𝑦 While halting criterion false do    (1) k = k +1    (2) 𝜽𝑘𝑖=𝐒𝑘−1(𝑖)⋅|⟨𝐫𝑘−1,𝐝𝑖⟩|    (3) 𝜆𝑘=argmax1≤𝑖≤𝑁{𝜽𝑘𝑖}    (4) construct upside-down tree 𝑇𝑘      𝐒𝑘(i) = {𝑐∶𝑖∈𝑇𝑘}    (5) Λ𝑘=Λ𝑘−1∪𝜆𝑘    (6) 𝐔𝑘=𝑃GS(Λ𝑘)    (7) 𝐫𝑘=𝐲−𝑃span{Λk}𝐲 end while 𝜶=𝐔𝑘𝐲 return 𝐱=ΦΨ𝜶

Algorithm 1 

Upside-down tree orthogonal matching pursuit.

For simplicity, sparsity level 𝐾 (assumed known) can be used as halting criterion here. If 𝐾 is unknown, we can modify the iteration in the above by letting 𝑘 run from 1 to 𝑀 but adding a threshold 𝜀 for |𝐫_𝑘|_𝟏 below which the iteration is terminated.

5. Experiment Result

To demonstrate the advantage of upside-down tree structure, we evaluated UDT-OMP algorithm by comparing the performance of OMP, TOMP, and UDT-OMP. We used piecewise smooth signal of length 𝑁=1024. 4-level Daubechies 8 wavelets were applied to sparsify the test signal, upward extending coefficient was set 𝑙=2, and weighting value 𝑐=𝑒 (natural number). Samples were obtained using a measurement matrix with i.i.d. Gaussian entries.

In the first experiment, we reconstructed the signal from 300 random measurements by OMP, TOMP, and UDT-OMP, respectively. Figure 2 depicts the reconstructions, the SNR of the OMP reconstructed signal is 27.63 dB, the SNR of the TOMP reconstructed signal is 34.57 dB, and the UDT-OMP reconstruction achieves SNR of 39.62 dB. In this experiment, UDT-OMP gains more than about 12 dB and 5 dB over OMP and TOMP, respectively. From experimental result, the advantage of the UDT-OMP method is clearly demonstrated.

(a)

(b)

(c)

(d)

Figure 2 

An example piecewise smooth signal of length 1024 and its reconstructions from 300 random measurements using OMP, TOMP, UDT-OMP. (a) Original signal. (b) OMP reconstructed signal, SNR = 27.63 dB. (c) TOMP reconstructed signal, SNR = 34.57 dB. (d) UDT-OMP reconstructed signal, SNR = 39.62 dB.

In the second experiment, we reconstructed the test signal from different numbers of measurements using OMP, TOMP, and UDT-OMP. Numbers of measurements over the range of 100 to 400 in increment of 30 are tested. 100 trials are repeated for each specific number, and the averaged reconstruction SNR is plotted in Figure 3. Through exploring the upside-down sparse tree structure, the UDT-OMP method outperforms OMP and TOMP in piecewise smooth signal reconstruction.

Figure 3 

Comparison between the performance of OMP, TOMP, and UDT-OMP.

6. Conclusions

This paper introduced an upside-down tree structure weighting scheme for OMP algorithm in wavelet domain CS signal reconstruction. The UDT-OMP weights the nodes that connect in the subtree with significant values. Different from tree structure presented in the previous CS recovery algorithms, UDT-OMP constructs the tree using upside-down structure rather than upright structure. UDT-OMP weights the projections that should have larger coefficients. The experimental results show that our method outperforms OMP and TOMP, and it can achieve more accurate approximation in piecewise smooth signal reconstruction. In this paper, we only considered the constant weight value. We can also adopt different weight values at different scales.

Acknowledgments

The work is partially supported by National Natural Science Foundation of China (no. 60827001) and China Scholarship Council (no. 2009607046).