# Recovery of Missing Samples with Sparse Approximations

Authors: Benjamin G. Salomon
Journal: ISRN Signal Processing
 (2013)

---

## 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.


---

## Body

## 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

Letx∈ℝN be a discrete signal with N samples. That is, x can be described by an N-dimensional real vector whose elements are denoted by x[0],x[1],…,x[N-1]. These elements correspond to the samples of the signal. Let Φ be an N×N orthonormal transform matrix; that is, ΦTΦ=ΦΦT=I, where I is the N×N identity matrix and ΦT is the transpose of Φ. The orthonormal matrix Φ defines an orthonormal transform of the vector x, denoted by x^∈ℝN, which is by definition the signal x^=Φx. The inverse transform of x^ is by definition the signal x=ΦTx^. In the complex case, the vectors belong to the space ℂN (space of N-dimensional complex vectors), and the conjugate transpose of the unitary matrix Φ has to be taken when computing the inverse transform.The columns ofΦT define an orthonormal basis of ℝN. That is, each signal x∈ℝN can be represented by a linear combination of the columns of ΦT, where the coefficients in this linear combination are given by the transform coefficients x^. Let ({ϕr}r=0,1,…,N-1) be the rows of ΦT (i.e., the columns of the matrix Φ ). It follows that for each x∈ℝN(1)x[k]=∑m=0N-1‍x^[m]ϕk[m],k=0,1,…,N-1.LetQ be a K-size proper subset of {0,1,…,N-1} and assume that only the K<N samples {x[k]}, k∈Q, are known. Consequently, the available K signal samples define  systems of K equations for the N transform coefficients:
(2)x[k]=∑m=0N-1‍x^[m]ϕk[m],k∈Q.
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 ΦT. That is, the signal is a linear combination of  certain subset of certain orthonormal basis vectors of ℝN.LetR be a P-size proper subset of {0,1,…,N-1}, where P≤K. A discrete bandlimited signal approximation to the signal x is obtained by
(3)xBL[k]=∑m∈R‍x^[m]ϕk[m],k=0,1,…,N-1.
That is, we assume that all transform coefficients with indices r∉R are equal to zero.For example, letΦ be the N×N discrete Fourier transform (DFT) matrix. That is, Φ is the N×N unitary matrix with elements (j is the imaginary unit)
(4)Φkl=1Ne-j(2π/N)kl.
A signal x is bandlimited if its DFT x^=Φx vanishes on some fixed nonempty proper subset S of {0,1,…,N-1}. The complement of S is the subset R. When the subset R is the set of 2M+1 elements (i.e., P=2M+1)
(5){0}∪{1,2,…,M}∪{N-M,…,N-1},
or equivalently (by periodicity)
(6){0,±1,±2,…,±M},
the signal is called a low-pass signal in DFT domain.Another example is the discrete cosine transform. LetΦ be the N×N discrete cosine transform (DCT) matrix. That is, Φ is the N×N orthonormal matrix with elements
(7)Φkl={1N,k=0,0≤l≤N-1,2Ncosπ(2l+1)k2N,1<k≤N-1,0≤l≤N-1.
A signal x is bandlimited if its DCT x^=Φx vanishes on some fixed nonempty proper subset S of {0,1,…,N-1}. When the complement of S (i.e., subset R) is the set of P elements
(8){0,1,2,…,P-1},
the signal is called a low-pass signal in DCT domain.Assuming thatx is a discrete bandlimited signal (i.e., x[k]=xBL[k], k=0,…,N-1), then, the available signal samples, x[k], k∈Q, can now be expressed as
(9)x[k]=∑m∈R‍x^[m]ϕk[m],k∈Q.
Equation (9) is a linear system consisting of K equations with P≤K unknowns. The existence and uniqueness of the solution depend, for any transform, on the subsets Q and R (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 K=P, we can determine the unique exact solution. When the system has a full column rank and K>P (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 fewunknown 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):
(10)x(k)=∑m=0N-1‍x^(m)ϕk(m),k∈Q,
where Q is the K-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. Letf∈ℝK be a signal of K samples that contains the given samples; that is, f[n]=x[kn], kn∈Q, n=0,1,…,K-1. Let φr, r=0,1,…,N-1, be the columns of the matrix ΦT. For each γ=0,1,…,N-1, let gγ=hγ/∥hγ∥ be a unit-norm vector of K samples, where hγ[n]=φγ[kn], kn∈Q, n=0,1,…,K-1, and the norm is the standard l2 norm in ℝK: ∥f∥2=∑n=0K-1‍f2[n]. The set of N unit-norm vectors gγ will be called a dictionary. The problem is to determine a sparse approximation to the vector f with the dictionary vectors gγ, γ=0,1,…,N-1.After determining a sparse solution usingL dictionary vectors:
(11)f≈∑i=1Lαγigγi=∑i=1Lβγihγi,
where
(12)βγi=αγi∥hγi∥,
the original signal x will approximated by
(13)x≈∑i=1L‍βγiφγi,
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ℝK of real-valued signals of size K. Let D={gγ}, γ=0,1,…,N-1, be a dictionary of N>K vectors, having a unit norm. Let Rnf be the residual of an n term approximation of a given signal f∈ℝK. MP subdecomposes the residue Rnf by projecting it on the dictionary vector that matches Rnf best. Starting from an initial approximation f0=0 and residue R0f=f, the MP algorithm builds up a sequence of sparse approximations stepwise.  MP begins by projecting f on a vector gγ0∈D and computing the residue R1f:
(14)f=〈f,gγ0〉gγ0+R1f,
where R1f is the residual vector after best approximating (l2 norm sense) R0f=f with gγ0. Since R1f is orthogonal to gγ0,
(15)∥f∥2=|〈f,gγ0〉|2+∥R1f∥2.
It follows that to minimize ∥R1f∥, we must choose gγ0∈D such that |〈f,gγ0〉| is maximum. The process is iterated on the residual. Suppose that Rnf for n≥0 is already computed. The next iteration chooses gγn such that
(16)gγn=argmingγ∈D∥Rn+1f∥2=argmingγ∈D∥Rnf-angγ∥2.
For unit-norm vectors, the last condition is equivalent to
(17)gγn=argmaxgγ∈D∥〈gγ,Rnf〉∥.
That is, the optimal vector is the one which best correlates with the residual. The MP projects Rnf on the chosen vector:
(18)Rnf=〈Rnf,gγn〉gγn+Rn+1f.
The orthogonality of Rn+1f and gγn implies
(19)∥Rnf∥2=|〈Rnf,gγn〉|2+∥Rn+1f∥2.
Summing (18) from n between 0 and L-1 (for any integer L>0) yields
(20)f=∑n=0L-1〈Rnf,gγn〉gγn+RLf.
Similarly, summation of (19) from n between 0 and L-1 gives
(21)∥f∥2=∑n=0L-1‍∥〈Rnf,gγn〉∥2+∥RLf∥2.
When dealing with finite-dimension signals (as in our case), ∥RLf∥ converges exponentially to 0 when L 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 gγn selected by the MP is a priori not orthogonal to the previously selected vectors gγp, p=0,1,…,n-1. After subtracting the projection of Rnf over gγn, new components are introduced in the direction of gγp, p=0,1,…,n-1. This is avoided by projecting the residues on an orthonormal set of vectors uγp, p=0,1,…,n-1, that span the same subspace that is spanned by gγp, p=0,1,…,n-1. That is, in each iteration, we determine the best approximation of the residual in the subspace spanned by uγp, p=0,1,…,n-1. When an approximation of sufficient accuracy is obtained, to expand f over the original dictionary vectors gγp, p=0,1,…,n-1, we perform a change of basis by expanding um in gγp, p=0,1,…,n-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:x=-0.5465φ7-0.8468φ19. 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 l2 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 x are given, the missing 30 samples are completely recovered by the OMP procedure.(a) Original signal (dash-dot: missing samples), (b) noisy signal (dash-dot: missing samples), (c) the vectorf containing the noisy given samples, and (d) original missing samples and the restored samples (dash-dot).
(a)
(b)
(c)
(d)


## 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.



---

*Source: 830723-2013-10-07.xml*