Recovery of Missing Samples with Sparse Approximations

Abstract

In most missing samples problems, the signals are assumed to be bandlimited. That is, the signals are assumed to be sparsely approximated by a known subset of the discrete Fourier transform basis vectors. We discuss the recovery of missing samples when the signals can be sparsely approximated by an unknown subset of certain unitary basis vectors. We propose the use of the orthogonal matching pursuit to recover missing samples by sparse approximations.

1. Introduction

Discrete signals are usually represented by their samples taken on a uniform sampling grid. However, in many applications, it may happen that some samples are lost or unavailable. In such cases, it is required to convert the irregularly sampled signal to a regularly sampled one, that is, to restore the missing samples. One approach for the recovery of missing data in discrete signals is based on the assumption that the underlying continuous-time signals are bandlimited. The celebrated sampling theorem by Whittaker [1], Kotel'nikov [2], and Shannon [3] implies that any continuous-time bandlimited signal can be reconstructed by its regularly spaced samples if the sampling frequency is higher than two times the maximum frequency component of the signal. The solution of the nonuniform sampling problem poses more difficulties and there exists a vast literature dealing with necessary and sufficient conditions for unique reconstruction and methods for reconstructing a function from its samples, for example, [4–6]. The numerical reconstruction methods, however, have to operate in a finite-dimensional model, whereas the theoretical results are usually derived for the continuous-time bandlimited functions (an infinite-dimensional subspace). The use of a truncated sampling series results with a finite-dimensional model but may lead, however, to severely ill-posed numerical problems [7].

Another approach is to address this problem in a finite-dimensional model of discrete bandlimited signals. A discrete bandlimited (DBL) signal has a sparse representation in terms of a certain unitary basis vector (e.g., discrete Fourier transform). That is, the signal can be represented by only a known subset of the unitary basis vectors. The recovery of missing samples is, in this case, equivalent to solving a linear system of equations [8]. In this paper, we focus on the recovery problem, when the a priori knowledge is that the signal is sparsely represented by an unknown subset of certain unitary basis vectors. This problem is much harder to solve, and there is an infinite number of possible solutions. Our approach is to choose the sparsest solution. That is, we are interested in approximating the known samples with a minimum number of basis vectors. Standard methods for solving linear systems of equations cannot provide the sparsest solution. We suggest the use of the orthogonal matching pursuit algorithm [9–11] for determining the sparsest approximation to the given samples from which the unknown samples can be determined.

2. Recovery of Missing Samples of Discrete Bandlimited Signals

Let 𝑥∈ℝ^𝑁 be a discrete signal with 𝑁 samples. That is, 𝑥 can be described by an 𝑁-dimensional real vector whose elements are denoted by 𝑥[0],𝑥[1],…,𝑥[𝑁−1]. These elements correspond to the samples of the signal. Let Φ be an 𝑁×𝑁 orthonormal transform matrix; that is, Φ^𝑇Φ=ΦΦ^𝑇=𝐼, where 𝐼 is the 𝑁×𝑁 identity matrix and Φ^𝑇 is the transpose of Φ. The orthonormal matrix Φ defines an orthonormal transform of the vector 𝑥, denoted by 𝑥∈ℝ^𝑁, which is by definition the signal 𝑥=Φ𝑥. The inverse transform of 𝑥 is by definition the signal 𝑥=Φ^𝑇𝑥. In the complex case, the vectors belong to the space ℂ^𝑁 (space of 𝑁-dimensional complex vectors), and the conjugate transpose of the unitary matrix Φ has to be taken when computing the inverse transform.

The columns of Φ^𝑇 define an orthonormal basis of ℝ^𝑁. That is, each signal 𝑥∈ℝ^𝑁 can be represented by a linear combination of the columns of Φ^𝑇, where the coefficients in this linear combination are given by the transform coefficients 𝑥. Let ({𝜙_𝑟}_{𝑟=0,1,…,𝑁−1}) be the rows of Φ^𝑇 (i.e., the columns of the matrix Φ ). It follows that for each 𝑥∈ℝ^𝑁𝑥[𝑘]=^𝑁−1_𝑚=0𝑥[𝑚]𝜙_𝑘[𝑚],𝑘=0,1,…,𝑁−1. (1)

Let 𝑄 be a 𝐾-size proper subset of {0,1,…,𝑁−1} and assume that only the 𝐾<𝑁 samples {𝑥[𝑘]}, 𝑘∈𝑄, are known. Consequently, the available 𝐾 signal samples define systems of 𝐾 equations for the 𝑁 transform coefficients: 𝑥[𝑘]=^𝑁−1_𝑚=0𝑥[𝑚]𝜙_𝑘[𝑚],𝑘∈𝑄. (2) The recovery of missing samples from an arbitrary signal is of course impossible: there is an infinite number of solutions to the underdetermined system of equations (2). The signals that occur in applications, however, often satisfy known constraints. The constraint used in this paper is called discrete band-limitedness: the signal can be represented by only a known subset of the columns of Φ^𝑇. That is, the signal is a linear combination of certain subset of certain orthonormal basis vectors of ℝ^𝑁.

Let 𝑅 be a 𝑃-size proper subset of {0,1,…,𝑁−1}, where 𝑃≤𝐾. A discrete bandlimited signal approximation to the signal 𝑥 is obtained by 𝑥^BL[𝑘]=_𝑚∈𝑅𝑥[𝑚]𝜙_𝑘[𝑚],𝑘=0,1,…,𝑁−1. (3) That is, we assume that all transform coefficients with indices 𝑟∉𝑅 are equal to zero.

For example, let Φ be the 𝑁×𝑁 discrete Fourier transform (DFT) matrix. That is, Φ is the 𝑁×𝑁 unitary matrix with elements (𝑗 is the imaginary unit) Φ_𝑘𝑙=1√𝑁𝑒^{−𝑗(2𝜋/𝑁)𝑘𝑙}. (4) A signal 𝑥 is bandlimited if its DFT 𝑥=Φ𝑥 vanishes on some fixed nonempty proper subset 𝑆 of {0,1,…,𝑁−1}. The complement of 𝑆 is the subset 𝑅. When the subset 𝑅 is the set of 2𝑀+1 elements (i.e., 𝑃=2𝑀+1) {0}∪{1,2,…,𝑀}∪{𝑁−𝑀,…,𝑁−1}, (5) or equivalently (by periodicity) {0,±1,±2,…,±𝑀}, (6) the signal is called a low-pass signal in DFT domain.

Another example is the discrete cosine transform. Let Φ be the 𝑁×𝑁 discrete cosine transform (DCT) matrix. That is, Φ is the 𝑁×𝑁 orthonormal matrix with elements Φ_𝑘𝑙=1√𝑁,𝑘=0,0≤𝑙≤𝑁−1,2𝑁𝜋()𝑘cos2𝑙+12𝑁,1<𝑘≤𝑁−1,0≤𝑙≤𝑁−1. (7) A signal 𝑥 is bandlimited if its DCT 𝑥=Φ𝑥 vanishes on some fixed nonempty proper subset 𝑆 of {0,1,…,𝑁−1}. When the complement of 𝑆 (i.e., subset 𝑅) is the set of 𝑃 elements {0,1,2,…,𝑃−1}, (8) the signal is called a low-pass signal in DCT domain.

Assuming that 𝑥 is a discrete bandlimited signal (i.e., 𝑥[𝑘]=𝑥^BL[𝑘], 𝑘=0,…,𝑁−1), then, the available signal samples, 𝑥[𝑘], 𝑘∈𝑄, can now be expressed as 𝑥[𝑘]=_𝑚∈𝑅𝑥[𝑚]𝜙_𝑘[𝑚],𝑘∈𝑄. (9) Equation (9) is a linear system consisting of 𝐾 equations with 𝑃≤𝐾 unknowns. The existence and uniqueness of the solution depend, for any transform, on the subsets 𝑄 and 𝑅 (see, e.g., [8] for a discussion on the recovery of samples when the signals are bandlimited in the discrete Fourier transform domain). When the system of equations (9) has a full column rank and 𝐾=𝑃, we can determine the unique exact solution. When the system has a full column rank and 𝐾>𝑃 (overdetermined system of equations), we will be interested in the least squares solution.

3. Sparse Approximations

In the previous section, we discussed the recovery of missing samples, when the band-limitedness was known a priori. That is, we had the a priori knowledge of which transform coefficients are equal to zero. The recovery of the missing samples is, in this case, equivalent to solving a linear system of equations for which many algorithms can be used. If we know that the signal is a low-pass bandlimited signal (e.g., in DFT or DCT domain), but the bandwidth is not given, we may start with a low bandwidth model and solve the resulting linear system of equations for increasing bandwidths until a sufficient accuracy in approximating the given samples is obtained. If the only a priori knowledge is that the signal is sparse in terms of the basis vectors, that is, the signal can be represented (or accurately approximated) with a few unknown basis vectors, the situation is very different. In this case, we have also to determine the basis vectors that sparsely approximate the signal.

We have to solve the underdetermined system of equations (2): 𝑥(𝑘)=^𝑁−1_𝑚=0𝑥(𝑚)𝜙_𝑘(𝑚),𝑘∈𝑄, (10) where 𝑄 is the 𝐾-size subset of indices of known samples. The underdetermined system of equations (2) has an infinite number of solutions. Since we know that the signal is sparsely approximated by the basis vectors, one reasonable approach is to choose the sparsest solution among the infinite number of solutions. That is, the optimal approximation is defined as either the sparsest approximation (i.e,. with the fewest basis vectors) that yields an approximation error that is smaller than a prespecified threshold or the approximation using a fixed number of basis vectors that minimizes the approximation error. Finding these approximations is an NP hard problem [12, 13].

We will use the following notations. Let 𝑓∈ℝ^𝐾 be a signal of 𝐾 samples that contains the given samples; that is, 𝑓[𝑛]=𝑥[𝑘_𝑛], 𝑘_𝑛∈𝑄, 𝑛=0,1,…,𝐾−1. Let 𝜑_𝑟, 𝑟=0,1,…,𝑁−1, be the columns of the matrix Φ^𝑇. For each 𝛾=0,1,…,𝑁−1, let 𝑔_𝛾=ℎ_𝛾/‖ℎ_𝛾‖ be a unit-norm vector of 𝐾 samples, where ℎ_𝛾[𝑛]=𝜑_𝛾[𝑘_𝑛], 𝑘_𝑛∈𝑄, 𝑛=0,1,…,𝐾−1, and the norm is the standard 𝑙₂ norm in ℝ^𝐾: ‖𝑓‖²=∑^𝐾−1_𝑛=0𝑓²[𝑛]. The set of 𝑁 unit-norm vectors 𝑔_𝛾 will be called a dictionary. The problem is to determine a sparse approximation to the vector 𝑓 with the dictionary vectors 𝑔_𝛾, 𝛾=0,1,…,𝑁−1.

After determining a sparse solution using 𝐿 dictionary vectors: 𝑓≈^𝐿_𝑖=1𝛼_𝛾𝑖𝑔_𝛾𝑖=^𝐿_𝑖=1𝛽_𝛾𝑖ℎ_𝛾𝑖, (11) where 𝛽_𝛾𝑖=𝛼_𝛾𝑖ℎ_𝛾𝑖, (12) the original signal 𝑥 will approximated by 𝑥≈^𝐿_𝑖=1𝛽_𝛾𝑖𝜑_𝛾𝑖, (13) from which the missing samples can be determined.

Several algorithms have been proposed for reducing the computational complexity by searching for efficient but nonoptimal approximations. The matching pursuit (MP) [10, 11] is a popular iterative greedy algorithm for approximate decomposition that addresses the sparsity issue directly. Vectors are selected one by one from the dictionary, picking at each iteration the vector that best correlates with the present residual, and thus optimizing the signal approximation (in terms of energy) at each iteration. An intrinsic feature of this algorithm is that when stopped after a few steps, it yields an approximation using only a few dictionary vectors.

We consider the space ℝ^𝐾 of real-valued signals of size 𝐾. Let 𝐷={𝑔_𝛾}, 𝛾=0,1,…,𝑁−1, be a dictionary of 𝑁>𝐾 vectors, having a unit norm. Let 𝑅_𝑛𝑓 be the residual of an 𝑛 term approximation of a given signal 𝑓∈ℝ^𝐾. MP subdecomposes the residue 𝑅_𝑛𝑓 by projecting it on the dictionary vector that matches 𝑅_𝑛𝑓 best. Starting from an initial approximation 𝑓₀=0 and residue 𝑅₀𝑓=𝑓, the MP algorithm builds up a sequence of sparse approximations stepwise. MP begins by projecting 𝑓 on a vector 𝑔_𝛾0∈𝐷 and computing the residue 𝑅₁𝑓: 𝑓=𝑓,𝑔_𝛾0𝑔_𝛾0+𝑅₁𝑓, (14) where 𝑅₁𝑓 is the residual vector after best approximating (𝑙₂ norm sense) 𝑅₀𝑓=𝑓 with 𝑔_𝛾0. Since 𝑅₁𝑓 is orthogonal to 𝑔_𝛾0, 𝑓²=𝑓,𝑔_𝛾0²+𝑅₁𝑓². (15) It follows that to minimize ‖𝑅₁𝑓‖, we must choose 𝑔_𝛾0∈𝐷 such that |⟨𝑓,𝑔_𝛾0⟩| is maximum. The process is iterated on the residual. Suppose that 𝑅_𝑛𝑓 for 𝑛≥0 is already computed. The next iteration chooses 𝑔_𝛾𝑛 such that 𝑔_𝛾𝑛=argmin_{𝑔𝛾∈𝐷}𝑅_𝑛+1𝑓²=argmin_{𝑔𝛾∈𝐷}𝑅_𝑛𝑓−𝑎_𝑛𝑔_𝛾². (16) For unit-norm vectors, the last condition is equivalent to 𝑔_𝛾𝑛=argmax_{𝑔𝛾∈𝐷}𝑔_𝛾,𝑅_𝑛𝑓. (17) That is, the optimal vector is the one which best correlates with the residual. The MP projects 𝑅_𝑛𝑓 on the chosen vector: 𝑅_𝑛𝑓=𝑅_𝑛𝑓,𝑔_𝛾𝑛𝑔_𝛾𝑛+𝑅_𝑛+1𝑓. (18) The orthogonality of 𝑅_𝑛+1𝑓 and 𝑔_𝛾𝑛 implies 𝑅_𝑛𝑓²=𝑅_𝑛𝑓,𝑔_𝛾𝑛²+𝑅_𝑛+1𝑓². (19) Summing (18) from 𝑛 between 0 and 𝐿−1 (for any integer 𝐿>0) yields 𝑓=^𝐿−1_𝑛=0𝑅_𝑛𝑓,𝑔_𝛾𝑛𝑔_𝛾𝑛+𝑅_𝐿𝑓. (20) Similarly, summation of (19) from 𝑛 between 0 and 𝐿−1 gives 𝑓²=^𝐿−1_𝑛=0𝑅_𝑛𝑓,𝑔_𝛾𝑛²+𝑅_𝐿𝑓². (21) When dealing with finite-dimension signals (as in our case), ‖𝑅_𝐿𝑓‖ converges exponentially to 0 when 𝐿 tends to infinity [10].

The approximations of MP are improved by orthogonalizing the directions of projection with a Gram-Schmidt procedure [9]. The resulting orthogonal MP (OMP) converges with a finite number of steps, which is not the case for a standard MP. The vector 𝑔_𝛾𝑛 selected by the MP is a priori not orthogonal to the previously selected vectors 𝑔_𝛾𝑝, 𝑝=0,1,…,𝑛−1. After subtracting the projection of 𝑅_𝑛𝑓 over 𝑔_𝛾𝑛, new components are introduced in the direction of 𝑔_𝛾𝑝, 𝑝=0,1,…,𝑛−1. This is avoided by projecting the residues on an orthonormal set of vectors 𝑢_𝛾𝑝, 𝑝=0,1,…,𝑛−1, that span the same subspace that is spanned by 𝑔_𝛾𝑝, 𝑝=0,1,…,𝑛−1. That is, in each iteration, we determine the best approximation of the residual in the subspace spanned by 𝑢_𝛾𝑝, 𝑝=0,1,…,𝑛−1. When an approximation of sufficient accuracy is obtained, to expand 𝑓 over the original dictionary vectors 𝑔_𝛾𝑝, 𝑝=0,1,…,𝑛−1, we perform a change of basis by expanding 𝑢_𝑚 in 𝑔_𝛾𝑝, 𝑝=0,1,…,𝑛−1.

We demonstrate the applicability of the approach by example. The original 128 samples length signal is the linear combination of two basis vectors of the DCT with randomly chosen coefficients: 𝑥=−0.5465𝜑₇−0.8468𝜑₁₉. The given signal is a noisy version of the original signal with additive white Gaussian noise with standard deviation equal to 0.05. We assume that 30 contiguous samples (9 to 38) of the noisy signal are missing. The results of the OMP are depicted in Figure 1, and we can see that the missing samples were approximately restored. The 𝑙₂ norm of the missing samples is 0.4339, while the norm of the error is 0.0315. If the samples of the original exact DBL signal 𝑥 are given, the missing 30 samples are completely recovered by the OMP procedure.

(a)

(b)

(c)

(d)

Figure 1 

(a) Original signal (dash-dot: missing samples), (b) noisy signal (dash-dot: missing samples), (c) the vector 𝑓 containing the noisy given samples, and (d) original missing samples and the restored samples (dash-dot).

4. Conclusions

In this paper, we addressed the recovery of missing samples of discrete bandlimited signals. We have suggested the use of the orthogonal matching pursuit to recover missing samples by sparse approximations, when the signal can be sparsely represented by an unknown subset of basis vectors of a certain orthonormal transform.