A Method for the Recovery of Gaps in General Analytic Signals

Abstract

We propose a generalized Papoulis-Gerchberg algorithm for the recovery of gaps in general analytic functions. The continuous-time algorithm is based on signal expansions in terms of the Chebyshev polynomials of the first kind, and the discrete-time implementation is based on a suitable nonuniform sampling scheme and the discrete cosine transform.

1. Introduction

Gapped (or missing) data is encountered when it is hard to obtain contiguous measurements of a function for a long time. The recovery of the missing function values is possible, at least in theory, given a certain a priori knowledge about the function, for example, analyticity on a certain interval. A real-valued function 𝑓, defined on a closed interval 𝐼 on the real line, is called analytic on 𝐼, if there exists an analytic extension of 𝑓 onto some open set 𝐺 of the complex plane ℂ that contains 𝐼 [1]. That is, there is a unique single-valued analytic function, defined on 𝐺, that coincides with 𝑓 on 𝐼. Let 𝑓 be analytic on 𝐼. When the values of 𝑓 on a subinterval of 𝐼 are given, 𝑓 can be determined on 𝐼, by means of analytic continuation.

In many missing data problems, the functions are assumed to be bandlimited. Let 𝐿²(ℝ) denote the Hilbert space of all square integrable functions defined on the real line ℝ. For 𝜏>0, and a function 𝑓∈𝐿²(ℝ), we define the operators 𝑃_𝜏 and 𝑄_𝜏 as follows:𝑃_𝜏𝑓𝑄=𝑓0,(𝑡),||𝑡𝑡||≤>𝜏𝜏,,_𝜏𝑓=𝑓−𝑃_𝜏𝑓.(1) Let ℱ𝑓=𝑓(𝜔)=∫^∞−∞𝑓(𝑡)𝑒^{−𝑗𝜔𝑡}𝑑𝑡, and, ℱ⁻¹𝑓 be the Fourier transform and inverse Fourier transform of 𝑓, respectively. For 𝜎>0, and a function 𝑓∈𝐿²(ℝ), we define the bandlimiting operator 𝑃_𝜎 as follows:𝑃_𝜎𝑓=ℱ⁻¹𝐿_𝜎ℱ𝑓=^∞−∞𝑓(𝑡−𝜏)sin𝜎𝜏𝜋𝜏𝑑𝜏,(2) where𝐿_𝜎(𝜔)=10,,||𝜔𝜔||≤>𝜎𝜎,.(3) The bandlimited signal extrapolation is defined as follows. Find 𝑓 having 𝑔=𝑃_𝜏𝑓 and knowing that 𝑓=𝑃_𝜎𝑓. Various methods were proposed for solving the bandlimited extrapolation problem; [2] contains a detailed comparison of some of the methods. The iterative bandlimited signal extrapolation algorithm of Papoulis and Gerchberg [3, 4] is attractive due to its relative simplicity, requiring only Fourier transform and inverse Fourier transform in each iteration. Based on the identity𝑔+𝑄_𝜏𝑃_𝜎𝑓=𝑃_𝜏𝑓+𝑄_𝜏𝑓=𝑓,(4) Papoulis [3] and Gerchberg [4] proposed the following algorithm:𝑚=I1n,it2i,a…liza:tio𝑓n_𝑚step:𝑓⁰=𝑔,=𝑔+𝑄_𝜏𝑃_𝜎𝑓^𝑚−1.(5) A proof for the convergence of the algorithm was given by Papoulis [3] using signal expansions in terms of the prolate spheroidal wave functions. Gerchberg [4] actually solved the dual problem of extrapolating the spectrum of a finite object beyond its diffraction limits, leading to superresolution and image enhancement.

A bandlimited function 𝑓(𝑡), defined on the real-line ℝ, is an entire function of exponential type when 𝑡 is extended to ℂ [5]. Hence, bandlimited signals are a very special case of analytic functions, and other analytic signals may be encountered in practice. In this paper, we propose a method for the recovery of gaps in general analytic signals using a generalization of the Papoulis-Gerchberg algorithm, based on signal expansions in terms of the Chebyshev polynomials of the first kind. In Section 2, we briefly review the properties of Chebyshev polynomials. We generalize in Section 3 the Papoulis-Gerchberg algorithm from continuous-time bandlimited signals to continuous-time signals in polynomials spaces and general analytic signals. We then focus, in Section 4, on the discrete implementation of the continuous-time algorithm, based on a suitable nonuniform sampling scheme and the discrete cosine transform. The performance of the recovery algorithm is demonstrated, in Section 6, by a numerical example.

2. Chebyshev Polynomials

The Chebyshev polynomial 𝑇_𝑛(𝑥) of the first kind is a polynomial in 𝑥 of degree 𝑛, defined by the relation𝑇_𝑛(𝑥)=cos𝑛𝜃when𝑥=cos𝜃,(6) for 𝑥∈[−1,1] [6]. The range of the corresponding variable 𝜃 is [0,𝜋] (these ranges are traversed in opposite directions). Intervals [𝑎,𝑏] other than [−1,1] are easily handled by the change of variables𝑠=2𝑥−(𝑎+𝑏)𝑏−𝑎.(7) Without loss of generality, we focus hereafter on the interval [−1,1].

Let 𝐿^2𝑤([−1,1]) be the Hilbert space of all real-valued square integrable functions on [−1,1] with respect to the nonnegative weight function𝑤(𝑥)=1√1−𝑥².(8) This is the space of functions 𝑓 such that the norm‖𝑓‖=¹⁻¹𝑤(𝑥)||𝑓(𝑥)||²𝑑𝑥^1/2(9) is finite. The associated inner product is⟨𝑓,𝑔⟩=¹⁻¹𝑤(𝑥)𝑓(𝑥)𝑔(𝑥)𝑑𝑥.(10) The set {𝑇_𝑖(𝑥),𝑖=0,1,…} forms a complete orthogonal polynomial system in 𝐿^2𝑤([−1,1]). The orthogonality relation is given by [6]⟨𝑇_𝑚(𝑥),𝑇_𝑛(𝑥)⟩==¹⁻¹𝑑0𝑤,(𝑥)𝑇𝑚_𝑚≠(𝑥𝑛),_𝑚𝑇_𝑛(𝑥)𝑑𝑥𝜋2,𝑚=𝑛,(11) where𝑑_𝑚=21,,𝑚𝑚=≥01,.(12) The 𝑁 zeroes of the Chebyshev polynomial 𝑇_𝑁(𝑥) are𝑥_𝑘=cos(2𝑘+1)𝜋2𝑁,𝑘=0,1,…,𝑁−1.(13) These zeroes are listed in decreasing order in 𝑥. The (𝑁−1)th degree polynomial 𝑝_𝑁−1(𝑥), interpolating a function 𝑓(𝑥) in the zeroes of 𝑇_𝑁, can be written as a sum of Chebyshev polynomials in the form𝑝_𝑁−1(𝑥)=^𝑁−1_𝑛=0𝑐_𝑛𝑇_𝑛(𝑥),(14) where the coefficients 𝑐_𝑛 in (14) are given by the explicit formula𝑐_𝑛=𝛾_𝑛𝑁−1_𝑘=0𝑓𝑥_𝑘𝑇_𝑛𝑥_𝑘,(15) where𝛾_𝑛=⎧⎪⎨⎪⎩𝑁12,𝑛=0,𝑁,0<𝑛≤𝑁−1.(16) It follows [7–9] that𝑐_𝑛=𝛾_𝑛𝑁−1_𝑘=0𝑓𝑥_𝑘cos𝜋(2𝑘+1)𝑛2𝑁,(17) which is the well-known discrete cosine transform (type II) [10] applied on the samples of the function 𝑓(𝑥) taken at a nonuniform sampling grid corresponding to the zeroes of 𝑇_𝑁. Efficient algorithms exist for computing the discrete cosine transform (see, e.g., [11] and references therein).

3. Gap Recovery in General Analytic Signals Based on A Chebyshev Polynomial Series Expansion

Let ℙ_𝑁−1 denote the subspace of all polynomials of degree <N (𝑁>1) in 𝐿^2𝑤([−1,1]). Let 𝑃_𝑁−1 denote the orthogonal projection operator onto ℙ_𝑁−1𝑃_𝑁−1𝑓=^𝑁−1_𝑘=0𝑓_𝑘𝑇_𝑘(𝑥),(18) where𝑓_𝑘=2𝜋𝑑_𝑘⟨𝑓,𝑇_𝑘⟩,(19)𝑑_𝑘 is defined in formula (12) and 𝑓∈𝐿^2𝑤([−1,1]).

Let 𝐶_𝑆 be the subspace of all the functions in 𝐿^2𝑤([−1,1]) which vanish outside a collection of disjoint closed intervals in [−1,1] and let 𝐶^+𝑆 be the orthogonal complement of 𝐶_𝑆 in 𝐿^2𝑤([−1,1]). We denote by 𝑃_𝑆 and 𝑄_𝑆 the orthogonal projection operators onto subspace 𝐶_𝑆 and 𝐶^+𝑆, respectively. It follows that 𝑄_𝑆=𝐼𝑑−𝑃_𝑆, where 𝐼𝑑 is the identity operator in 𝐿^2𝑤([−1,1]). If we assume that 𝑓=𝑃_𝑁−1𝑓, then we have:𝑓=𝑃_𝑆𝑓+𝑄_𝑆𝑓=𝑃_𝑆𝑓+𝑄_𝑆𝑃_𝑁−1𝑓(20) The continuous-time recovery problem can be stated as follows. Find 𝑓 having 𝑔=𝑃_𝑆𝑓 and knowing that 𝑓=𝑃_𝑁−1𝑓. We propose the following Papoulis-Gerchberg type algorithm for recovering gaps in polynomial signals in 𝐿^2𝑤([−1,1]) (PGP algorithm).

The PGP Algorithm
𝑚I=nit1i,a2li,z…atio:n𝑓st_𝑚ep:𝑔=𝑓⁰=𝑃_𝑆𝑓,(21)=𝑔+𝑄_𝑆𝑃_𝑁−1𝑓^𝑚−1.(22) The iteration process is stopped when ‖𝑓^𝑚−𝑓^𝑚−1‖≤𝜖, where 𝜖>0 is a user specified threshold. 𝑓^𝑚 or 𝑃_𝑁−1𝑓^𝑚 is taken as the result of the recovery algorithm after 𝑚 iterations.
If 𝑓∈ℙ_𝑁−1, then 𝑓=𝑃_𝑁−1𝑓, and it follows that 𝑔=𝑓−𝑄_𝑆𝑃_𝑁−1𝑓.(23) We know from (22) and (23) that 𝑓^𝑚−𝑓=𝑄_𝑆𝑃_𝑁−1𝑓^𝑚−1−𝑓.(24) Therefore, ‖𝑓^𝑚−𝑓‖≤‖‖𝑄_𝑆‖‖‖‖𝑃_𝑁−1‖‖‖‖𝑓^𝑚−1−𝑓‖‖≤‖‖𝑓^𝑚−1−𝑓‖‖,(25) which proves that the PGP algorithm has the property that it reduces the error energy during the iteration process.

Lemma 1. A signal 𝑓∈ℙ_𝑁−1 can be uniquely determined from 𝑃_𝑆𝑓 by using the PGP algorithm (21), and (22).

Proof. Iterating (22), we obtain 𝑓^𝑚==𝑔𝑔++𝑄𝑄_𝑆𝑆𝑃_𝑁−1𝑓^𝑚−1=𝑔+𝑄_𝑆𝑃_𝑁−1𝑔+𝑄_𝑆𝑃_𝑁−1𝑓^𝑚−2=_𝑚𝑃_𝑁−1𝑔+𝑄_𝑆𝑃_𝑁−1𝑄_𝑆𝑃_𝑁−1𝑔+𝑄_𝑆𝑃_𝑁−1𝑓^𝑚−3_𝑛=0𝑄_𝑆𝑃_𝑁−1^𝑛𝑔.(26) Substituting (23) in (26), we obtain 𝑓^𝑚=𝑓−𝑄_𝑆𝑃_𝑁−1^𝑚𝑓.(27) According to von Neumann's theorem on alternating projections [12], for every 𝑓∈ℙ_𝑁−1, lim_𝑚→∞𝑄_𝑆𝑃_𝑁−1^𝑚𝑓=𝑓_𝑐,(28) where 𝑓_𝑐 is the projection of 𝑓 onto the closed subspace 𝐶^+𝑆∩ℙ_𝑁−1. Any 𝑓∈ℙ_𝑁−1 is analytic on [−1,1] and can vanish in an interval if and only if it is the zero function. Hence, the subspace 𝐶^+𝑆∩ℙ_𝑁−1 contains only the zero function. We obtain that 𝑓_𝑐=0 and that the PGP converges in norm to the required function 𝑓.

𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1 is a bounded positive self-adjoint compact operator. Hence, there exists an orthonormal basis of ℙ_𝑁−1 consisting of eigenfunctions 𝜙_𝑛, 𝑛=0,1,…,𝑁−1, of the operator 𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1, with corresponding real positive eigenvalues 𝜆_𝑛, such that 1>𝜆₀>𝜆₁>⋯>𝜆_𝑁−1>0 [13].

Lemma 2. Let 𝑓=𝜙_𝑘 be an eigenfunction of 𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1 with corresponding eigenvalue 𝜆_𝑘, and let 𝑔=𝑃_𝑆𝑓 be the given values of 𝑓. Then, 𝑃_𝑁−1𝑓^𝑚=𝐴_𝑚𝜙_𝑘where𝐴_𝑚=1−1−𝜆_𝑘^𝑚+1.(29)

Proof. We will prove (29) by induction. It is true for 𝑚=0 that 𝑃_𝑁−1𝑓⁰=𝑃_𝑁−1𝑔=𝑃_𝑁−1𝑃_𝑆𝜙_𝑘=𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1𝜙_𝑘=𝜆_𝑘𝜙_𝑘=𝐴₀𝜙_𝑘,(30) where we used the identity 𝑃_𝑁−1𝜙_𝑘=𝜙_𝑘 (since 𝜙_𝑘∈ℙ_𝑁−1).
Suppose that it is true for some 𝑚≥1. By using (22), we obtain 𝑃_𝑁−1𝑓^𝑚+1=𝑃_𝑁−1𝑔+𝑃_𝑁−1𝑄_𝑆𝑃_𝑁−1𝑓^𝑚=𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1𝜙_𝑘+𝑃_𝑁−1𝐼_𝐿^2𝑤_([−1,1])−𝑃_𝑆𝑃_𝑁−1𝑓^𝑚=𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1𝜙_𝑘+𝑃_𝑁−1𝑓^𝑚−𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1𝑓^𝑚.(31) Substituting 𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1𝜙_𝑘=𝜆_𝑘𝜙_𝑘 and using the induction assumption, we obtain 𝑃_𝑁−1𝑓^𝑚+1=𝜆_𝑘𝜙_𝑘+𝐴_𝑚𝜙_𝑘−𝜆_𝑘𝐴_𝑚𝜙_𝑘=𝐴_𝑚+1−𝐴_𝑚𝜆_𝑘𝜙_𝑘.(32) In (31) and (32), we used the idempotent property of the orthogonal projection operator 𝑃_𝑁−1 (i.e., 𝑃_𝑁−1=𝑃_𝑁−1𝑃_𝑁−1). It is easy to see that 𝐴_𝑚+1−𝐴_𝑚𝜆_𝑘=1−1−𝜆_𝑘^𝑚+1+1−𝜆_𝑘^𝑚+1𝜆_𝑘=1−1−𝜆_𝑘^𝑚+2.(33) Hence, 𝑃_𝑁−1𝑓^𝑚+1=𝐴_𝑚+1𝜙_𝑘,(34) which completes the proof.

Any 𝑓∈ℙ_𝑁−1 can be represented as𝑓=^𝑁−1_𝑘=0𝛼_𝑘𝜙_𝑘,(35) where𝛼_𝑘=¹⁻¹𝑤(𝑥)𝑓(𝑥)𝜙_𝑘(𝑥)𝑑𝑥.(36) Using (29) and (35), we obtain that after 𝑚 iterations of the PGP algorithm𝑃_𝑁−1𝑓^𝑚=^𝑁−1_𝑘=0𝛼_𝑘1−1−𝜆^𝑘^𝑚+1𝜙^𝑘.(37) The energy of error is given by‖‖𝑓−𝑃_𝑁−1𝑓^𝑚‖‖²=^𝑁−1_𝑘=0𝛼^2𝑘1−𝜆_𝑘^2(𝑚+1).(38) It follows that the rate of convergence of the PGP algorithm is dependent on the spectral representation of the signal in terms of the eigenfunctions of the operator 𝑃_𝑁−1𝑃_𝑆𝑃_𝑁−1. If the signal contains only significant components which correspond to relatively high eigenvalues, the convergence is fast. If this is not the case, the PGP will be slowly convergent.

The PGP can also be viewed as an example of signal restoration by the method of projections onto convex sets (POCS) [14, 15]. The two convex sets, or constraints, are the affine set 𝑉₁, which contains all the functions 𝑓∈𝐿^2𝑤([−1,1]) with 𝑔=𝑃_𝑆𝑓 and the subspace ℙ_𝑁−1. If 𝑓∈𝑉₁ but 𝑓∉ℙ_𝑁−1, we have two inconsistent convex constraints. The convergence of a POCS algorithm with two inconsistent convex constraints is a well-studied problem [16, 17]. In this case, the PGP algorithm iterates 𝑃_𝑁−1𝑓^𝑚 converge to the signal 𝑞∈ℙ_𝑁−1 closest to 𝑉₁, that is, with the minimum distance 𝑑, where 𝑑 is given by𝑑𝑞,𝑉₁=inf_𝑦∈𝑉1‖𝑞−𝑦‖.(39) Thus, the PGP algorithm is exact for all 𝑓∈P_𝑁−1, and a certain “best polynomial approximation” in the sense of (39) is obtained when 𝑓∉ℙ_𝑁−1.

The expansion of analytic function in terms of the Chebyshev polynomials converges rapidly (exponentially). It can be proved that if the function 𝑓(𝑥) can be extended to a function 𝑓(𝑧) analytic on the ellipse {𝑧∶|𝑧+√𝑧²−1|=𝑟}, where 𝑟>1, then |𝑓−𝑃_𝑁−1𝑓|=𝑂(𝑟^{−(𝑁−1)}) for all 𝑥 in [−1,1] [6]. This property of rapid convergence combined with the approximation property (39) makes the proposed algorithm a practical gap recovery tool for general analytic functions.

4. Discrete Implementation of the Continuous-Time-Recovery Algorithm

A Chebyshev polynomial-based series expansion allows a convenient and accurate discrete implementation by a nonuniform sampling scheme, where the samples are taken at time locations that are Chebyshev polynomial zeroes and a discrete cosine transform (DCT) applied on the nonuniform samples [8, 9]. Equation (15) is a Gauss-Chebyshev quadrature for the original integral (19), which is known to be a very accurate numerical method [6]. In addition, if the maximum absolute error is used as the optimality criterion, then the interpolating polynomial given by (14) provides optimal approximation among polynomials of degree = 𝑁−1 to the function 𝑓 [18].

Let ℝ^𝑁 denote the real 𝑁-dimensional space. The 𝑙₂ inner product is given by⟨𝑦,ℎ⟩=^𝑁−1_𝑚₌₀𝑦[𝑚]ℎ[𝑚],(40) and the induced norm is given by‖𝑥‖²=⟨𝑥,𝑥⟩.(41) Let 𝐶 be the 𝑁×𝑁 unitary DCT matrix defined as𝐶_𝑚,𝑛=⎧⎪⎪⎨⎪⎪⎩1√𝑁,𝑚=0,0≤𝑛≤𝑁−1,2𝑁cos𝜋(2𝑛2𝑁+1)𝑚,10<≤𝑚𝑛≤≤𝑁𝑁−−11.,(42) The DCT of the vector 𝑦∈ℝ^𝑁 (arranged as a column vector) is given by 𝑦̂=𝐶𝑦. The inverse DCT of the vector 𝑦̂ is given by 𝑦=𝐶^𝑇𝑦̂, where 𝐶^𝑇 is the transpose of 𝐶. A discrete signal 𝑦 is bandlimited if its DCT 𝑦̂ vanishes on some fixed set, that is, if𝑦̂[𝑚]=0,𝑚∈𝑆,(43) where 𝑆 is a fixed nonempty proper subset of {0,1,…,𝑁−1}. The set of signals bandlimited to a specific set 𝑆 is a linear subspace of ℝ^𝑁, of dimension equal to the cardinal of the complement of 𝑆.

Let 𝑦[𝑘]=𝑓(𝑥_𝑘), 𝑘=0,…,𝑁−1 be a vector in ℝ^𝑁, whose elements are the samples of the analytic function 𝑓(𝑥) taken in the zeroes of 𝑇_𝑁 (13). Since only 𝑃_𝑆𝑓 is available, only certain 𝑄<𝑁 samples are known. Assuming that 𝑓∈ℙ_𝑀−1 and 𝑀≤𝑄, the discrete recovery problem is to determine the DCT coefficients 𝑦[̂𝑚] for 𝑚=0,1,…,𝑀−1 from the 𝑄 known samples. Using the DCT coefficients (with an appropriate simple scaling) in (14), we obtain the (𝑀−1)th degree polynomial that approximates 𝑓(𝑥).

Let 𝑆₁ be a subset of size 𝑄 in [0,1,2,…,𝑁−1] containing the indices of the known samples, let 𝑆₂ denote the subset of indices [0,1,…,𝑀−1], and denote 𝐷=𝐶⁻¹=𝐶^𝑇. The required DCT coefficients can be determined by solving the least squares problem_𝑖_∈𝑆2𝐷_𝑘𝑖𝑦̂[𝑖]≃𝑦[𝑘],𝑘∈𝑆₁.(44) It can be proved [19] that any square submatrix (𝑀=𝑄) of the matrix 𝐷 in (44) is the product of a Vandermonde and upper triangular matrix, and therefore, it is invertible. It follows that the submatrix 𝐷_𝑘𝑖, 𝑖∈𝑆₂, 𝑘∈𝑆₁, is a full column rank matrix when 𝑀≤𝑄. The resulting least squares problem (44) can be solved by many methods (e.g., singular value decomposition, QR factorization, conjugate-gradient method for solving the normal equations, etc.) [20, 21]. A simple popular option (not necessarily the best in terms of computational complexity) is a version of the discrete Papoulis-Gerchberg algorithm to be described below.

Let the operator 𝑃_𝐴 map a signal 𝑦∈ℝ^𝑁 into another, 𝑃_𝐴𝑦=𝐴𝑦, by multiplying 𝑦 with a diagonal matrix 𝐴 containing only zeroes or ones. 𝐴_𝑚𝑚=1, if 𝑚∈𝑆₁, and 𝐴_𝑚𝑚=0, otherwise. It is easy to see that this operation is idempotent and self-adjoint. Hence, 𝑃_𝐴 is the orthogonal projector onto the subspace of ℝ^𝑁 signals with zero-valued elements at the indices not belonging to 𝑆₁.

Let the operator 𝑃_𝐵 map a signal 𝑦 into another 𝑃_𝐵𝑦=𝐶^𝑇𝐵𝐶𝑦, where 𝐵 a diagonal matrix containing only zeroes or ones. 𝐵_𝑚𝑚=1, if 𝑚∈𝑆₂, and 𝐵_𝑚𝑚=0, otherwise. Using the identities 𝐵²=𝐵 and 𝐵^𝑇=𝐵, it follows that 𝑃_𝐵 is idempotent𝑃_𝐵𝑃_𝐵𝑦=𝐶^𝑇𝐵𝐶𝐶^𝑇𝐵𝐶𝑦=𝐶^𝑇𝐵𝐶𝑦=𝑃_𝐵𝑦,(45) and self-adjoint ⟨𝑃_𝐵𝑦,ℎ⟩==𝐶𝑦,^𝑇𝐶𝐵_𝑇𝐶𝑦,ℎ=⟨𝐵𝐶𝑦,𝐶ℎ⟩=⟨𝐶𝑦,𝐵𝐶ℎ⟩𝐵𝐶ℎ=⟨𝑦,𝑃_𝐵ℎ⟩,(46) for every 𝑦,ℎ∈ℝ^𝑁. Hence, 𝑃_𝐵 is the orthogonal projector onto the subspace of bandlimited (in terms of the DCT and 𝑆₂) signals in ℝ^𝑁.

A discrete version of the PGP algorithm is obtained by working in the finite-dimensional space ℝ^𝑁 and using the operators 𝑃_𝐴 instead of 𝑃_𝑆 and 𝑃_𝐵 instead of 𝑃_𝑁−1. The discrete version of the PGP is also a projection onto convex sets algorithm, and it follows that the iterates 𝑃_𝐵𝑦^𝑚 converge to the required least squares solution (see the discussion before (39)).

5. Discussion

Our model is based on two assumptions. First, the original analytic signal can be accurately approximated by a linear combination of a few Chebyshev polynomials (exponential convergence), and, secondly, that the Gauss-Chebyshev quadrature (15) provides a good approximation to the original integral (19). When these assumptions are not (approximately) satisfied (e.g., due to noise), the algorithm may converge (it always converges) to a crude approximation of the original signal. In such cases, a higher-order polynomial approximation and a denser sampling grid would provide an improvement. If the model assumptions are (approximately) satisfied, the resulting approximation is very accurate and minimizes a certain objective function in the sense of (39).

A Chebyshev polynomial series expansion is a Fourier cosine series with a change of variable 𝑥=cos𝜃 for 𝑓(cos𝜃) [6]. While 𝑓(𝑥) itself is not periodic in 𝑥, the function 𝑓(cos𝜃) is periodic in 𝜃. The smoother the function is, the more rapidly its Fourier coefficients decrease. Since 𝑓(cos𝜃) is periodic and analytic, its Fourier series must have exponential convergence. For a Fourier cosine series expansion for the nonperiodic function 𝑓(𝑥), the first derivative of the function is discontinuous at the interval's borders and for a Fourier series expansion, the function itself is discontinuous at the interval’s borders. It follows that a Chebyshev polynomial series expansion of an analytic function provides a better approximation than a Fourier series expansion (using a discrete Fourier transform in discrete implementations [22]) with the same number of terms, which makes it attractive to recovery problems. The improved approximation comes, however, at the price of a nonuniform sampling grid. The sampling points are clustered near the boundaries of the interval [−1,1]. It follows that the proposed algorithm is more suited to the recovery of gaps that occur in the middle of the signal or nearby. In such cases, we benefit from both the excellent approximation properties of the Chebyshev polynomials and the fact that less missing samples have to be recovered.

6. Numerical Example

We demonstrate the performance of the discrete implementation of the recovery algorithm by a gap recovery example. The function 𝑓(𝑥) defined on the interval [−1,1]𝑓(𝑥)=11+25𝑥².(47) It is analytic on [−1,1] but has a pole at 𝑥=±0.2𝑗, where 𝑗=√−1. We assume that the piece [−0.4,0.4] is missing and that signal is well approximated by 64 Chebyshev polynomials. The samples are taken at the 128 zeroes of the polynomial 𝑇₁₂₈, and it follows that 34 samples are missing. The proposed algorithm is used to recover the missing samples: the 𝑙₂ norm of the error is 0.16, while the norm of the missing samples is 3.5244. We see that the recovery algorithm succeeds in approximating the missing data. Finally, the function and the resulting polynomial approximation are given at a uniform grid of 128 points in the interval [−1,1] in Figure 1.

Figure 1 

Original function: solid, approximation by a Chebyshev polynomial expansion: dashdot.

7. Summary and Conclusions

The recovery of gaps in general analytic functions was obtained by using signal expansions in terms of the Chebyshev polynomials of the first kind. The recovery algorithm is exact for polynomial signals, and certain best polynomial approximations are obtained for general analytic signals. For the continuous-time case, we have extended the well-known Papoulis-Gerchberg algorithm for bandlimited signals to signals in polynomial subspaces. The discrete implementation is based on a specific nonuniform sampling grid (at the zeroes of a Chebyshev polynomial) and the discrete cosine transform and results in solving a linear least squares problem.