Laurent Gosse
email: l.gosse@area.ba.cnr.it

Practical band-limited extrapolation relying on Slepian series and compressive sampling

Keywords: Band-limited extrapolation, Prolate spheroidal wave functions, Slepian series, analytic continuation, finite Fourier transform, compressive sampling, sparse and compressible signals recovery. 1991 MSC: 42A16, 41A58, 65D15, 68P30

We consider a rather simple algorithm to address the fascinating field of numerical extrapolation of (analytic) band-limited functions. It relies on two main elements: namely, the lower frequencies are treated by projecting the known part of the signal to be extended onto the space generated by "Prolate Spheroidal Wave Functions" (PSWF, as originally proposed by Slepian), whereas the higher ones can be handled by the recent so-called "Compressive Sampling" (CS, proposed by Candès) algorithms which are independent of the largeness of the bandwidth. Slepian functions are recalled and their numerical computation is explained in full detail whereas Compressive Sampling techniques are summarized together with a recent iterative algorithm which has been proved to work efficiently on so-called "compressible signals" which appear to match rather well the class of smooth bandlimited functions. Numerical results are displayed for both numerical techniques and the accuracy of the process consisting in putting them altogether is studied for several test-signals.

Introduction

Band-limited functions are lying at the very core of various questions arising in many areas of application, especially signal processing, information and approximation theories. Let us begin by recalling that any function f belonging to L1 (R) or L 2 (R) is said to be band-limited if there exists a positive number σ such that:

∀t ∈ R, f (t) = 1 2π σ -σ
f (ξ) exp(itξ).dξ, the Fourier transform f being defined according to the functional space 1 of f . The number σ ∈ R + is called the bandwidth of f ; it is not an obvious task to decide whether or not real-life signals are endowed with a finite bandwidth, see [START_REF] Slepian | On bandwidth[END_REF]. The classical Paley-Wiener theorem states that any such function belonging to L 2 (R) extends to the whole complex plane as an entire function of exponential type; more precisely, there holds (see e.g. [START_REF] Zayed | On the notion of bandlimitedness and its generalizations[END_REF])

∀z ∈ C, |f (z)| ≤ sup t∈R |f (t)| exp(σℑ(z)).
As a consequence of analytic continuation theory for functions of complex variable, the knowledge of such a function restricted to any arbitrary interval of R allows to deduce all the remaining values corresponding to any z ∈ C.

Despite this very neat theory, numerical extrapolation of band-limited functions remains a delicate practical problem having quite a long history; several survey papers are available, like for instance [START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Moore | Prolate Spheroidal wave functions, an introduction to the Slepian series and its properties[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF][START_REF] Zayed | On the notion of bandlimitedness and its generalizations[END_REF]. Two main approaches coexist relying on different algorithmic methodologies: iterative methods generally follow the classical Gerchberg-Papoulis algorithm [START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF][START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF], but recently, a renewed interest grew around non-iterative methods relying on the so-called Slepian series [START_REF] Slepian | Prolate spheroidal wave functions, Fourier analysis and uncertainty. I[END_REF][START_REF] Slepian | Prolate spheroidal wave functions, Fourier analysis and uncertainty, IV[END_REF]. Iterative routines for scale-limited extrapolation have been developed around similar ideas and led to a so-called "generalized Gerchberg-Papoulis algorithm": see [START_REF] Lin | On theory and regularization of scale-limited extrapolation[END_REF][START_REF] Lin | On the convergence of wavelet-based iterative signal extrapolation algorithms[END_REF][START_REF] Xia | Signal extrapolation in wavelet subspaces[END_REF]. In the present text, we shall not deal with iterative extrapolation routines for band-limited functions; instead we refer to existing papers [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF][START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF][START_REF] Youla | Image Restoration by the Method of Convex Projections: Part 1 Theory[END_REF] and references therein. It turns out that, for a long time, iterative methods have gained from the fact that extrapolation algorithms based on Slepian's Prolate functions (outlined in the paper [START_REF] Slepian | Some comments on Fourier analysis, uncertainty and modeling[END_REF], page 388) were far beyond practical computing possibilities. However, the situation has changed dramatically nowadays and several papers have been recently published to present feasible algorithms for the computation of Slepian series: see for instance [START_REF] Boyd | Prolate spheroidal wavefunctions as an alternative to Chebyshev and Legendre polynomials for spectral element and pseudospectral algorithms[END_REF][START_REF] Gubner | A Simple Method for Computing Projections onto Subspaces of Prolate Spheroidal Wave Functions[END_REF][START_REF] Khare | Sampling theory approach to prolate spheroidal wave functions[END_REF][START_REF] Kirby | Calculation of spheroidal wave functions[END_REF][START_REF] Karoui | New efficient methods of computing the prolate spheroidal wave functions and their corresponding eigenvalues[END_REF][START_REF] Walter | A new friendly method of computing prolate spheroidal wave functions and wavelets[END_REF][START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF]. We aim at taking advantage of them in order to propose a realizable and efficient numerical scheme to extrapolate accurately smooth band-limited functions.

Compressive Sensing differs greatly from the aforementioned algorithms in the sense that it doesn't rely on the straightforward projection of a recorded signal onto the specific basis of some Hilbert functional space; its way of proceeding goes more like correlating it with a collection of randomly chosen but perfectly known test-signals. Hence Compressed Sensing is less naturally an extrapolation method because it needs unknown values to be located on the complementary of a random set of points. In most cases, unknown values are located in a connected domain (or two of them, each one around the edges of some observations interval). At this level, what saves the day is that CS can reconstruct signals which are strongly under-sampled: hence in many cases one can manage a random set of observations not seeing the whole interval containing the unknown values. For instance, in case there is a collection of 100 measures and it is to deduce the 10 ones located on the right, one needs for CS to work one particular random sequence of measures taken inside the whole set of 110 values but seeing only the 100 first ones. This can easily happen when one needs to measure only a very small subset of values, the smallness being directly related to the sparse character of (for instance) the Fourier transform of the signal to extrapolate. This is in this sense that we propose here to exploit CS algorithms to handle efficiently extrapolation problems involving high frequencies which typically aren't easy to process by the standard Gerchberg-Papoulis method; to the best of our knowledge, this idea seems to be new.

The paper is organized as follows. In §2.1, we recall useful theoretical considerations concerning the derivation of Prolate functions and present an extrapolation algorithm following Slepian [START_REF] Slepian | Some comments on Fourier analysis, uncertainty and modeling[END_REF]. An illustration of these functions together with their L 2 norm is given in §2.

"Adiabatic decoupling" between lower and higher frequencies

Part of our methodology relies on the assumption that a quite general class of interesting signals allow for being split between a "trend" corresponding to lower frequencies and a "fluctuation" which accounts for higher frequencies. Moreover, we shall also assume in all what follows that higher frequencies don't have much effect on low ones during a limited amount of time. Put in another way, extrapolating on an interval of moderate length permits to uncouple between low and high frequencies inside the signal. This is quite appealing because modern numerical strategies to handle economically low and high frequencies are quite different and such a decoupling opens the way of taking the best from each of them to treat both parts of the observed signal. We used the terminology "adiabatic decoupling" in a loose sense to underline the fact that high frequencies are separated from the rest of the signal and that they are not supposed to have a sensible effect on the lower ones during a limited amount of time; this is strongly reminiscent of the adiabatic theory as encountered in quantum mechanics, for instance in the context of molecular dynamics where light particles are driven by the movement of heavier ones (consult e.g. [START_REF] Spohn | Adiabatic decoupling and time-dependent Born-Oppenheimer theory[END_REF] for details in this direction).

Moreover, this framework is also well suited for other multiscale phenomena being of a given size d, but which are dominated by only the c ≪ d trend's Fourier coefficients (see e.g. [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Zhou | Theoretical and experimental analysis of randomized algorithm for sparse Fourier transform analysis[END_REF]). In this paper, we shall limit ourselves to check whether or not such a frequency splitting is appropriate or not based on the elementary signals correlation. More precisely, assume that we know a signal t → f (t) ∈ R inside some observation interval [a, b] and there is a possible decomposition f = f lo + f hi such that the convex hull of the support of flo and the support of fhi are disjoint. Then this splitting will be admissible if moreover,

κ := b a f lo (t)f hi (t).dt f lo L 2 (a,b) f hi L 2 (a,b) (1)
is lower that some threshold value ε 0 ≪ 1 in modulus. By Cauchy-Schwarz inequality, we have that κ ∈ [-1, 1] with |κ| = 1 if and only if there exists α ∈ R such that f lo = αf hi .

2.1 The theory of Prolate Spheroidal Wave Functions (PSWF)

We begin by presenting a theoretical framework adapted to the processing of low frequencies.

Several equivalent definitions

We now briefly resume the derivation and the main properties of Slepian's Prolate functions, following various survey articles, see [START_REF] Karoui | New efficient methods of computing the prolate spheroidal wave functions and their corresponding eigenvalues[END_REF][START_REF] Larsson | On prolate spheroidal wave functions for signal processing[END_REF][START_REF] Moore | Prolate Spheroidal wave functions, an introduction to the Slepian series and its properties[END_REF][START_REF] Slepian | Some comments on Fourier analysis, uncertainty and modeling[END_REF][START_REF] Walter | A new friendly method of computing prolate spheroidal wave functions and wavelets[END_REF][START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF]. Originally, these functions were sought as possible real-valued solutions of the maximum concentration problem for bandlimited functions raised by Shannon himself (see [START_REF] Moore | Prolate Spheroidal wave functions, an introduction to the Slepian series and its properties[END_REF], first page): for any T ∈ R + , consider the highest value of

0 ≤ α(T ) 2 = T 2 -T 2 |f (t)| 2 dt R | f (ξ)| 2 dξ
.

Assuming f to be σ-bandlimited, this expression rewrites as: ( * denotes complex conjugate)

α(T ) 2 = T 2 -T 2 dt σ -σ f (ξ) exp(itξ)dξ σ -σ f (ξ ′ ) * exp(itξ ′ )dξ ′ σ -σ | f (ξ)| 2 dξ = σ -σ dξ σ -σ sin T (ξ -ξ ′ ) π(ξ -ξ ′ ) f (ξ) f (ξ ′ ) * dξ ′ σ -σ f (ξ) f (ξ) * dξ .
The last step is completed by integration in the t variable. A maximizing f has to satisfy the following Fredholm equation of the second kind:

∀ξ ∈ [-σ, σ], σ -σ sin T (ξ -ξ ′ ) π(ξ -ξ ′ ) f (ξ ′ )dξ ′ -α(T ) 2 f (ξ) = 0.
After easy changes of variables, this reduces to the well-known eigenvalue problem:

∀ξ ∈ [-1, 1], 1 -1 sin c(ξ -ξ ′ ) π(ξ -ξ ′ ) ψ(ξ ′ )dξ ′ = λψ(ξ), λ ≥ 0, c = T σ. (2)
This presentation clearly shows that Prolate functions are indeed objects defined in the Fourier domain, even if many authors use the generic variables x or t to write them down, possibly suggesting the opposite. As the "sinc" kernel is symmetric positive definite, the integral operator is compact and equation (2) has countably many solutions λ n , ψ n , n ∈ N, which satisfy furthermore:

• λ n ≥ 0, λ 0 ≥ λ 1 ≥ ... ≥ λ n ≥ ... and λ n → 0 as n → +∞,
• real-valued eigenfunctions ψ n constitute a complete orthogonal system of L 2 (-1, 1).

We stress that besides the index n, eigenvalues and eigenfunctions λ n , ψ n do heavily depend on the parameter c. At this point, it has been qualified as a "lucky accident" by Slepian himself that these Prolate functions ψ n also furnish solutions bounded on (-1, 1) of the eigenvalue problem for the following differential operator:

∀x ∈ [-1, 1], - d dx (1 -x 2 ) dψ dx + c 2 x 2 ψ = χψ, χ > 0. ( 3
)
This fact is of great utility when it comes to computing numerically the functions ψ n as discretizing directly the convolution equation (2) yields a matrix which coefficients decay very slowly to zero; hence standard algorithms won't be able to furnish accurately its eigenvalues and eigenvectors. If one plugs c = 0 inside (3), one obtains the differential equation leading to the orthogonal polynomial basis of Legendre on [-1, 1] with χ n = n(n + 1). For this reason, Prolate functions with c > 0 are sometimes seen as a "bandlimited generalization of Legendre polynomials"; one has to be careful with this terminology as the Legendre basis is not primarily defined in the Fourier domain.

2.1.2 Orthogonal bases of both L 2 (-1, 1) and L 2 (R)

Another amazing property is that, contrary to Legendre polynomials, Prolate functions can be easily extended so as to furnish an orthogonal system in L 2 (R). Following Slepian, we define for any n ∈ N:

∀|ξ| > 1, ψ n (ξ) = 1 λ n 1 -1 sin c(ξ -ξ ′ ) π(ξ -ξ ′ ) ψ n (ξ ′ )dξ ′ .
They can be normalized so as to satisfy:

R ψ m (ξ)ψ n (ξ)dξ = δ m,n , 1 -1 ψ m (ξ)ψ n (ξ)dξ = λ n δ m,n , (4)
with δ m,n the Kronecker symbol. With such a normalization, eigenvalues can be interpreted as the fraction of the total mass which is concentrated in the interval [-1, 1]. Eigenfunctions constitute therefore an orthonormal basis of a subset of L 2 (R) called the Paley-Wiener space, constituted of c-bandlimited functions endowed with a finite L 2 (R) norm. It has been believed for long time that Slepian's functions were the unique ones to possess this surprising property: another example of such a system has been found only recently in [START_REF] Zayed | A generalization of the prolate spheroidal wave functions[END_REF].

A property which reveals itself of crucial importance for extrapolation goes as follows: any Prolate function ψ n also satisfies

∀t ∈ R, 1 -1 exp(ictξ)ψ n (ξ)dξ = µ n ψ n (t), µ n ∈ C. (5)
Besides emphasizing the bandlimited character of these functions, this feature expresses that, up to a change of scale, the function ψ n has the same shape as its Fourier transform. The values of µ n and λ n are related through:

λ n (c) = c 2π |µ n (c)| 2 . ( 6
)
Concerning extrapolation purposes, the expression (5) opens the way to completing the program outlined by Slepian in [START_REF] Slepian | Some comments on Fourier analysis, uncertainty and modeling[END_REF]. Namely, it allows for expressing the value of ψ n (t), for any t ∈ R, relying only on the knowledge of this function restricted in the interval [-1, 1], restriction that can be computed accurately thanks to the eigenvalue problem for the differential operator given in [START_REF] Boyd | Prolate spheroidal wavefunctions as an alternative to Chebyshev and Legendre polynomials for spectral element and pseudospectral algorithms[END_REF]. So, as a consequence of this, since we know that ψ n constitute an orthogonal system of both L 2 (-1, 1) and of the subset of L 2 (R) made of c-bandlimited functions, a practical extrapolation routine for a c-bandlimited function f known only in [-1, 1] will work as follows:

• compute as accurately as possible λ n and ψ n restricted to [-1, 1],

• express f in the orthogonal basis made of ψ n through the formula:

∀t ∈ [-1, 1], f (t) = n∈N f n ψ n (t), f n = 1 λ n 1 -1 f (τ )ψ n (τ ).dτ. (7)
• extend the definition of ψ n by means of (5) to any other interval of R,

• the extrapolation f of f is then given by:

∀t ∈] -∞, -1[∪]1, +∞[, f (t) = n∈N f n ψ n (t) = n∈N f n µ n 1 -1 exp(ictξ)ψ n (ξ)dξ. ( 8
)
The passage from [-1, 1] to R could be also achieved by means of the formula (2) which may look more interesting because in this last case, only the nonnegative λ n values are needed; however, we found that, relying on the algorithms proposed in [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF], the computing cost of the µ n values was equivalent to the one of the λ n . Therefore, the aforementioned extrapolation algorithm is completely operational as soon as one has reliable numerical algorithms to compute the Prolate functions ψ n , the complex values µ n and the scalar products f n for any c-bandlimited function f known on [-1, 1] only.

Numerical examples

We shall present some reliable algorithms in §3 of the present paper. Before that, we display some numerical simulations to illustrate the behavior of both Prolate functions and their corresponding eigenvalues. On the left of Fig. 1, we display the first even/odd Prolate functions on the interval [-1, 1] for the bandwidth parameter c = 7. For this moderate value, we observe that the first eigenfunctions, besides looking rather similar to ordinary Legendre polynomials, are quite spread across the interval [-1, 1]. Their absolute values grow strongly when approaching the edges ±1. The corresponding eigenvalues λ n are arranged like a step function (which explains why trying to solve directly a discretized version of (2) cannot be considered a good way of deriving Slepian functions) with a sharp transition zone between zero and one. It is a well-known fact that this transition occurs around the value n * such that λ n * ≃ c; accordingly, this value of c has been drawn with a vertical red line on both figures 1 and 2. Figure 2 deals with the higher value c = 21. In this case, the bandwidth is larger and consequently, Prolate functions are more peaked and more concentrated inside the interval [-1, 1]. For moderate values of n ∈ N, their moduli remain bounded in the vicinity of the edges ±1, which also differs from the preceding case. Eigenvalues λ n still display their "step function" profile, but the numerical transition zone is now located before the value c, which reflects some difficulties in the task of computing them accurately as c grows (see e.g. [START_REF] Karoui | New efficient methods of computing the prolate spheroidal wave functions and their corresponding eigenvalues[END_REF]).

Some important results of Compressed Sampling (CS) theory

Compressive Sampling (or Compressed Sensing) consists now in a set of firmly founded mathematical algorithms which allow to reconstruct a so-called sparse signal out of a surprisingly restricted number of its observations, regardless to its possibly high frequencies. Roughly speaking, the required number of observed values is proportional only to sparsity, that is to say the amount of nonzero coefficients one needs in order to describe the signal in a supposedly well-suited base; of course, a signal can be sparse in a given basis and not in another one. For an introductory text, see http://compressedsampling.com/.

The "Restricted Isometry Property" (RIP)

Here, we shall not attempt to give the more advanced theoretical results of Compressed Sampling; rather we shall try to convey the main ideas underlying the algorithms and state some of the results having far-reaching consequences for our purposes. Details are to be found in the original papers [START_REF] Candès | Compressive sampling[END_REF][START_REF] Candès | The restricted isometry property and its implications for compressed sensing C[END_REF][START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies[END_REF].

A signal f ∈ L 2 (R) (or belonging to any Hilbert functional space) is said to be s-sparse in the basis e i (t), i ∈ N if the synthesis of f (t) for any t ∈ R involves only a small number s ∈ N of elements in this basis. That is,

∀t ∈ R, f (t) = s i=0
α i e i (t), s small.

In the discrete case, given a suitable basis e i , i ∈ {1, 2, ..., d} of R d , we shall deal with s-sparse vectors v which satisfy v = s i=0 α i e i and s ≪ d. For any x ∈ R d , we also define its s-sparse approximation x s as the vector made of the s largest entries of x. Of course, any vector v ∈ R d being fixed, one can always consider an ad hoc basis such that its coordinates read (1, 0, ...0); here we work with the unitary Inverse Discrete Fourier Transform (IDFT) reading (exp(i(k -1)(j -1)/d)/ √ d) k=1,...,d, j=1,...,d or the canonical basis of R d .

This being said, it makes sense to state the following definition:

Definition 1 (Candès [START_REF] Candès | The restricted isometry property and its implications for compressed sensing C[END_REF]) For any integer d ≥ s ∈ N, the restricted isometry constant δ s of a (possibly rectangular) matrix Φ is defined as the smallest number such that it holds for every s-sparse vector x: (we recall that for 1 ≤ p < +∞,

x ℓ p = ( d i=0 |x i | p ) 1 p ) (1 -δ s ) x 2 ℓ 2 ≤ Φ x 2 ℓ 2 ≤ (1 + δ s ) x 2 ℓ 2 . (9)
A vector x ∈ R d is said to be s-sparse if it has at most s ≤ d non-zero entries.

A matrix endowed with this Restricted Isometry Property for a certain integer s is constituted such that every set of of its d columns with cardinality less than s approximately behaves like an orthonormal system. This means in particular that if δ s ≪ 1, then all the submatrices of Φ made of s columns are well conditioned; moreover, ( 9) is a frame inequality [START_REF] Christensen | An introduction to frames and Riesz bases[END_REF] restricted to s-sparse vectors for the frame defined as the set of d columns of Φ.

Recovery theorems: sparse and noisy cases

From this, Candès and his collaborators established the following recovery theorems:

Theorem 1 (Candès [START_REF] Candès | The restricted isometry property and its implications for compressed sensing C[END_REF]) Consider y ∈ R N , N ≤ d, a matrix Φ which satisfies δ 2s < √ 2-1, and x opt ∈ R d the solution of the following ℓ 1 minimization process:

inf x∈R d x ℓ 1 , Φ x = y ∈ R d . (10)
Then the following error estimates hold:

x opt -x ℓ 1 ≤ C x s -x ℓ 1 , x opt -x ℓ 2 ≤ C √ s x s -x ℓ 1 .
In particular, in case x ∈ R d is s-sparse, the recovery is exact.

A stability result is also available in case the observed vector y is corrupted by noise:

Theorem 2 (Candès [6]) Consider y ∈ R N , N ≤ d, a matrix Φ which satisfies δ 2s < √ 2-1, a noise vector z ∈ R N such that Φ x = y + z, z ℓ 2 ≤ ε and x opt ∈ R d the solution of the convex minimization process: inf x∈R d x ℓ 1 , y -Φ x ℓ 2 ≤ ε. ( 11
)
Then the following error estimate holds with constants C, C ′ being explicitly computable:

x opt -x ℓ 2 ≤ C √ s x s -x ℓ 1 + C ′ ε.
One issue remaining with both these results is that the user has to know in advance whether or not the observed vector y ∈ R N is noisy because depending on the answer to this question, one of the two programs [START_REF] Christensen | An introduction to frames and Riesz bases[END_REF] or (11) will be selected.

Remark 1

The RIP (9) as advocated in both Theorems 1 and 2 has a particular meaning when Φ is any N × d rectangular matrix extracted from the d × d unitary IDFT matrix. The requirement δ 2s < √ 2 -1 means that for any 2s-sparse vector of Fourier coefficients ẑ ∈ C d of a (generally non-sparse) vector z ∈ C d , there holds:

(2 - √ 2) z 2 ℓ 2 (C d ) < Φẑ 2 ℓ 2 (C N ) < √ 2 z 2 ℓ 2 (C d ) ,
where

ẑ ℓ 2 (C d ) = z ℓ 2 (C d
) by Parseval's equality and Φẑ = z with z ∈ C N is obtained by retaining only the coordinates of z corresponding to the N rows of Φ. So, this is a way of expressing that the original signal z cannot be strongly concentrated on the set of d -N remaining points; on the contrary, it has to be rather uniformly spread on the whole set of d points, which is a consequence of ẑ being sparse since the Fourier basis, having infinite support in the time variable, is strongly delocalized. This idea is even reinforced in case the set of N rows is "generic", in the sense that it is chosen equally probably as any other set of N rows among the available d ones: this is somewhat the reason behind exact recovery.

Numerical examples

The Restricted Isometry Property [START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies[END_REF], which is actually crucial when it comes to underdetermined recovery, is verified by few classes of matrices only. Well-known examples (see [START_REF] Candès | Compressive sampling[END_REF]) are as follows:

• Gaussian measurements: Entries of N × d matrix Φ are independently sampled from the normal distribution (0, 1 N ). Then if N ≥ Cs log(d/N ), Φ meets the requirements of Theorems 1 and 2 with probability 1 -C exp(-γd), γ > 0.

• Fourier measurements: Rows of N × d matrix Φ are randomly chosen among the ones of the unitary IDFT matrix. The resulting columns are normalized so as to have norm 1 in ℓ 2 (and the entries of the resulting vector x opt may have to be modified accordingly after the minimization process in case a FFT is necessary afterward). It is recalled in [START_REF] Needell | Uniform Uncertainty Principle and signal recovery via Regularized Orthogonal Matching Pursuit[END_REF] that for N ≥ Cs log(d) γ (γ = 6 or γ = 4), Theorem 1 holds with overwhelming probability; it is conjectured that γ = 1 suffices in [START_REF] Candès | Compressive sampling[END_REF].

As an illustration, we display on Fig. 3 the recovery of a random signal in R 256 involving 20 spikes (so s = 20 and d = 256 so this signal is sparse in the canonical basis of R 256 ). We selected a random Gaussian measurement matrix with N is taken as the integer part of s log(d): the reconstruction can be considered as perfect. Comparison with a standard ℓ 2 minimization is also shown for which the numerical outcome is not sparse at all.

3 Extrapolation of low frequencies with Prolate functions

Numerical computation of Prolate functions inside [-1, 1]

The numerical computation of Slepian functions relies on the solving of the differential Sturm-Liouville equation (3) which is a perturbation of the one admitting the Legendre polynomials P r as eigenfunctions. Relying on this last feature, it sounds natural to seek ψ n as a Legendre series of the form: The polynomials Pr are called normalized Legendre polynomials and satisfy 12) inside (3) and rearranging leads to an eigenvalue problem for a symmetric matrix A; more precisely, we have the following result.

ψ n (t) = r∈N β n,r Pr (t), Pr (t) = P r (t) r + 1/2. ( 12
)
1 -1 | Pr (t)| 2 dt = 1. Plugging (
Theorem 3 (Xiao et al. [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF]) For any value of c > 0, the numbers χ n and β n,. are the eigenvalues and eigenvectors of the endomorphism of ℓ 2 represented by the infinite matrix A whose entries are zero except:

A r,r = r(r + 1) + 2r(r + 1) -1 (2r + 3)(2r -1) c 2 , A r,r+2 = A r+2,r = (r + 2)(r + 1) (2r + 3) (2r -1)(2r + 5) c 2 .
Thanks to its particular structure, the matrix A can be split into two submatrices A even and A odd containing only even-numbered and odd-numbered rows and columns respectively. However, truncation of these infinite matrices isn't straightforward because their entries don't decay quickly with increasing rows and columns indexes. But thanks to the rapid decay in r of the scalar products | 1 -1 ψ n (t)P r (t).dt|, see Theorem 3.4 in [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF], the eigenvectors of interest can be obtained considering only the leading rows and columns of A even and A odd .

Hence a feasible algorithm to compute ψ n through the Legendre series representation [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF] goes as follows: solving the eigenvalue problems following Theorem 3:

(A even -χ 2n .Id) β 2n,. = 0, (A odd -χ 2n+1 .Id) β 2n+1,. = 0, n ∈ N,
and then generating efficiently the values of the normalized Legendre polynomials Pr (t) at particular abscissas. Especially, it is interesting to derive values of Pr at points t k , k = 1, 2, ..., K corresponding to Legendre-Gauss quadrature [START_REF] Laurie | Computation of Gauss-type quadrature formulas[END_REF]. For computing the nodes t k and weights ω k of Gaussian quadrature rules, the fundamental tool is the three-term recurrence relation satisfied by the set of orthogonal polynomials associated to the corresponding weight function. In case P r is the monic orthogonal polynomial of degree r, it is a well-known fact that such orthogonal polynomials are related through a three-terms recurrence relation ∀t ∈ [-1, 1], P r+1 (t) + a r P r (t) + b r P r-1 (t) = tP r (t), P 0 ≡ 1, P -1 ≡ 0.

One introduces therefore the so-called Jacobi K × K matrices reading,

J 1 =         a 0 1 0 • • • • • • 0 b 1 a 1 1 0 • • • 0 0 b 2 a 2 1 0 . . . 0 . . . . . . . . . . . . . . .         , J 2 = 0 √ b 2 a 2 √ b 3 0 . . . 0 . . . . . . . . . . . . . . .        
, and if P := (P 0 (t k ), P 1 (t k ), P 2 (t k ), ...P K-1 (t k )) T , then J 1 P = t k P . Following [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF] (see also [START_REF] Laurie | Computation of Gauss-type quadrature formulas[END_REF]), the nodes of a Gaussian quadrature can actually be computed as the eigenvalues of the tridiagonal symmetric matrix J 2 , and also the weights through the corresponding eigenvectors. This is the approach proposed by Gubner [START_REF] Gubner | A Simple Method for Computing Projections onto Subspaces of Prolate Spheroidal Wave Functions[END_REF] in order to compute efficiently any scalar product of the form which is exactly what is required for (7),

1 -1 f (t)ψ n (t).dt ≃ K k=1 ω k f (t k ) 2N,2N +1 r=0,1 β n,r Pr (t k ), N ∈ N, (13)
with f any smooth function defined on [-1, 1] and (ω k , t k ) k=1,...,K the weights and nodes of a Legendre-Gauss quadrature rule (see Fig. 4). The summation on index r is carried out according to the parity of n; Matlab codes are also furnished in the papers [START_REF] Gubner | A Simple Method for Computing Projections onto Subspaces of Prolate Spheroidal Wave Functions[END_REF] and [START_REF] Laurie | Computation of Gauss-type quadrature formulas[END_REF]. Because of [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF] and the convention which normalizes eigenvectors with unit ℓ 2 norm, both [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF] and [START_REF] Gubner | A Simple Method for Computing Projections onto Subspaces of Prolate Spheroidal Wave Functions[END_REF] present algorithms which involve Prolate functions normalized so as to be of unit norm in the interval [-1, 1], on the contrary of (4).

Eigenvalues: how to pass from [-1, 1] to any abscissa of R

We saw before in §2.1.2 that one of the key points to complete the extrapolation process is the efficient computation of eigenvalues appearing in either (2) or (5); following [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF] and taking advantage of our knowledge of the functions ψ n restricted to the interval [-1, 1], we present now a concrete way to complete this step. According to formula [START_REF] Candès | Compressive sampling[END_REF], it suffices to discretize the Fourier integral to produce the value of ψ n (x) for any x ∈ R. This can be done accurately exploiting the Gaussian quadrature rule derived in the preceding section; namely, we write that:

∀x ∈ R, ψ n (x) ≃ K k=1 ω k exp(ic x t k ) 2N,2N +1 r=0,1 β n,r Pr (t k ),
which is just a special case of [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF]. All the problem of setting up (8) reduces therefore in extracting accurately the eigenvalues µ n in [START_REF] Candès | Compressive sampling[END_REF]. Following [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF], we recall an useful result:

Theorem 4 (Xiao et al. [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF]) For any value of c > 0 and m, n ∈ N 2 such that m = n (mod 2), there holds:

|µ m | 2 |µ n | 2 = 1 -1 ψ ′ n (t)ψ m (t).dt 1 -1 ψ n (t)ψ ′ m (t).dt .
It becomes possible to compute the eigenvalues by induction: first, λ 0 can be derived through

µ 0 ψ 0 (t = 0) = 1 -1 ψ 0 (ξ)dξ.
Then, all the successive values of |µ n |, n ≥ 1 can be obtained relying on Theorem 4 by using again the Gaussian quadrature derived in §3.1. At last, we recall from [START_REF] Xiao | Prolate spheroidal wave functions, quadrature and interpolation[END_REF] that µ n = i n |µ n | and it remains to use [START_REF] Candès | The restricted isometry property and its implications for compressed sensing C[END_REF]. The figures 1 and 2 have been produced this way.

Remark 2 Another algorithm to derive eigenvalues µ n follows from the primary knowledge of Prolate functions ψ n obtained from solving the diagonalization problem raised by Theorem 3: from the numerical values ψ n (t k ) known on the Gaussian quadrature points, one can compute through [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF] the Fourier integral in ( 5) by the same quadrature rule. Thus it suffices to divide and compute an average to derive each corresponding eigenvalue. Numerical results are similar to those coming from the algorithm by Xiao et al.

Numerical results

A simple test-case

We aim at completing a 15% extrapolation of the following simple and C ∞ (R) (but not rigorously bandlimited) function,

f pro (t) = 1 2 exp(-2t 2 ) cos(5πt), t ∈ [-1, 1], (14)
that is to say, finding its numerical values on the interval ]1, 1.3] by means of its Prolate expansion following equations ( 7), [START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF] and the numerical procedures explained in the preceding sections. The result is shown in Fig. 5 with the choice c = 30 for the bandwidth. One can observe that for this value, the compression is very good since roughly 10 coefficients acting on even Prolate functions (f pro is clearly even) suffice to get a very good approximation: absolute errors inside the interval [-1, 1] are of the order of 10 -14 : this agrees with the values indicated in Table 1 and Theorem 3.11 of [START_REF] Shkolniskya | Approximation of bandlimited functions[END_REF]. Coefficients for odd Prolate functions are of the order of 10 -16 . When it comes to the extrapolation, errors grow no more than values around 10 -3 , which is still satisfying. Gaussian quadrature was used with 256 nodes.

First composite problem

We are now interesting in extrapolating more involved signals, for which we shall split between lower and higher frequencies, each one being treated by means of a specific algorithm. The examination of absolute errors made on the global problem will be deferred to §5. First, we consider the following smooth function for t ∈ [-1, 1]:

f (t) = exp t 2 sin(πt) f lo (t)
+ 0.15[sin(13πt) -0.5 cos(23πt) + 0.3 cos(7πt)]

f hi (t) . ( 15
)
In order to check whether or not the splitting between low and high frequencies is admissible, we computed the empirical correlation between these two parts following [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] where integrals are obtained with the Gaussian quadrature like (13): we found that for a = -1, b = 1, κ ≃ 0.0014 which can be considered as satisfying. The first part is to be treated by the extrapolation algorithm based on Prolate functions with the bandwidth parameter c = 7.

The numerical outcome for an extrapolation of 15% (up to t = 1.3) is shown in Fig. 6. Prolate coefficients are equally spread among even and odd PSWF. The reconstruction error inside the interval [-1, 1] is very low (of the order of 10 -16 ); the extrapolation error is bigger, but remains at acceptable levels (of the order of 10 -3 ). These results can be considered as being satisfying.

Second composite problem

As a second test-case, we consider the quite intricate smooth function:

∀t ∈ [-1, 1], f (t) = exp(-2t 2 ) exp(3t) sin(πt) cos(3πt) f lo (t) + 0.5[sin(5πt) -cos(7πt)] f hi (t)
. results are displayed on Fig. 7. One can see that reconstruction errors inside [-1, 1] are still very small, but extrapolation errors for 1 < t ≤ 1.3 are bigger compared to the preceding example since they are now of the order of 10 -2 . Few Prolate coefficients are needed to represent the function ( 7) despite the fact it is not rigorously band-limited because of the Gaussian function involved. Results are still satisfying though.

4 Extrapolation of higher frequencies with Compressed Sampling

The class of compressible signals

Roughly speaking, compressible signals correspond to a subclass for which there exists a convenient Hilbert basis allowing to represent each of them with a very limited number of coefficients; in this case, one says that these signals are "well compressed" in this basis. This notion is clearly dependent on the choice of the aforementioned basis since, in the case one considers Fourier or wavelet bases, smoothness in the time variable implies compressibility in both these bases. With the Fourier transform, it is possible to give a precise meaning to the subclass of "signals being compressible in the Fourier basis":

Proposition 1 Let f ∈ C p (R) and suppose that f (i) belongs to L 1 (R) or to L 2 (R), then for all ξ ∈ R, | f (ξ)| ≤ C p (1 + |ξ|) -p for a suitable constant C p ∈ R + .
Consequently, the class of compressible signals for which the approximate recovery theorems hold contains all the signals having a certain number of derivatives when expressed within the temporal variable. A property of the same flavor would also hold for wavelet bases.

Hence we can introduce a discrete analogue of this notion:

Definition 2 Let x ∈ R d ;
we say that it is a compressible vector in a given basis if its coefficients decay according to a power law. More precisely, there exists p ≥ 1 such that for any k ≤ d, the modulus of the k th largest coefficient of x is bounded by C p k -p .

Being correctly oversampled, bandlimited functions clearly yield discrete compressible signals in the Fourier basis since their Fourier transform has compact support. The first result concerning recovery of randomly sampled but compressible signals has been given in [START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies[END_REF]:

Theorem 5 (Optimal recovery of weak-ℓ p , [START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies[END_REF]) Let x ∈ R d a compressible vector in the sense of Definition 2 with p > 1 or x ℓ 1 ≤ C 1 for p = 1, and let α > 0 be a sufficiently small number. Assume we are given N random measurements Φ x, then with probability 1, the minimizer x opt of (10) is unique. Furthermore, with probability at least 1 -O(d -ρ α ), the following error bound holds:

x -x opt ℓ 2 ≤ C p,α C p log d N p-1 2 . ( 17
)
Remark 3 Theorem 5 doesn't refer to the RIP despite an earlier but similar concept appears in its proof. Its conclusion is probabilistic, as were the ones of the results proved in [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF]: it says that given O(s log d) measurements, one is able to retrieve through ℓ 1 minimization an approximation which is as good as if the s largest entries of x were known in advance. Moreover, the recovery algorithm (10) doesn't need to know the quantities p and C p .

Regularized Orthogonal Matching Pursuit (ROMP)

Both theorems stated in §2.2.2 rely on the numerical solving of a ℓ 1 minimization program. Despite the fact powerful methods now exist for completing this task, we prefer to turn to lighter iterative solvers, like the so-called ROMP routine introduced in [START_REF] Needell | Uniform Uncertainty Principle and signal recovery via Regularized Orthogonal Matching Pursuit[END_REF][START_REF] Needell | Signal recovery from incomplete and inaccurate measurements via Regularized Orthogonal Matching Pursuit[END_REF]. In this section, we recall two recovery theorems for ROMP, adapted to the kind of problems we want to process; the measurement N × d matrix is still denoted by Φ.

Theorem 6 (Needell et al. [START_REF] Needell | Uniform Uncertainty Principle and signal recovery via Regularized Orthogonal Matching Pursuit[END_REF]) Consider y ∈ R N , N ≤ d and the matrix Φ which satisfies δ 2s ≤ 0.03/ √ log s. For any x ∈ R d with y = Φ x, ROMP outputs a set I such that supp( x) ⊂ I, |I| ≤ 2s in at most s iterations.

This case corresponds to exact sparse recovery; we note that Theorem 2.1 in [START_REF] Needell | Uniform Uncertainty Principle and signal recovery via Regularized Orthogonal Matching Pursuit[END_REF] gives rigorous bounds in terms of RIP for widely used random measurement matrices. Concerning the stability when measurements or the original (supposedly sparse) signal is perturbed, we state the result in the context of compressible signals for which no algorithmic changes are required (contrary to [START_REF] Christensen | An introduction to frames and Riesz bases[END_REF] and ( 11)):

Theorem 7 (Needell et al. [START_REF] Needell | Signal recovery from incomplete and inaccurate measurements via Regularized Orthogonal Matching Pursuit[END_REF]) Consider y ∈ R N , N ≤ d and the matrix Φ which satisfies δ 8s ≤ 0.01/ √ log s. For any compressible vector x ∈ R d in the sense of Definition 2 with y = Φ x + ε, ROMP produces an approximation to x satisfying:

x romp -x ℓ 2 ≤ log s Cp s p-1 2 + C ε ℓ 2 . ( 18
)
We used the same letter p in both Proposition 1 and Definition 2 on purpose since the discrete Fourier transform converts smoothness into fast coefficients decay. Hence, assuming we deal only with smooth signals, and under a strong RIP restriction, Theorem 7 gives us an error control which compares with [START_REF] Khare | Sampling theory approach to prolate spheroidal wave functions[END_REF] onto the output of the ROMP algorithm in a case where the original vector x isn't sparse. The authors of [START_REF] Needell | Signal recovery from incomplete and inaccurate measurements via Regularized Orthogonal Matching Pursuit[END_REF] conjecture that the logarithmic factor in front of ε is unnecessary, and that a bound like the one of Theorem 2 actually holds. We point out that a Matlab code for ROMP is available at the following location: http://www.math.ucdavis.edu/ dneedell/romp.m

The choice of convenient random matrices

The first articles [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies[END_REF] insist on the fact that the spectacular results of CS stem from the combined use of random sensing on the one hand and sparse recovery through ℓ 1 minimization on the other hand. In order to achieve good numerical quality, the randomness of Φ is crucial; however, in the context of extrapolation, it is not completely obvious how to meet this requirement. Since Φ is obtained by randomly selecting N = O(s log d) rows inside the unitary d×d IDFT matrix and further renormalizing its columns so as to be of unit ℓ 2 norm, one can easily understand that not all random sequences can be used. Actually, we only know a fraction of d obs < d values; moreover, the d -d obs unknown values are all packed at the end of the observations interval in a strongly non-random manner. Hence, in order to set up correctly the CS algorithms, we must work with sensing matrices Φ obtained with those particular random sequences of N rows among d which don't perceive the unknown d -d obs values located on the right of the observations interval. Clearly, this puts a limitation on the size of the ratio α = d/d obs ≥ 1 as it will become more and more difficult to find such sequences for increasing α's while keeping constant the number N = O(s log d). We stress that according to the results of [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies[END_REF], the simple idea of selecting a random sequence of N rows among the d obs first ones (corresponding to the observed values) isn't convenient; indeed, the selected chosen measurement ensemble of N values has really to be equally likely as any other one (see [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF], page 492 and Remark 1).

Remark 4 It is interesting to observe that compared to the classical extrapolation techniques based on variants of Gerchberg-Papoulis (GP) algorithm [START_REF] Gerchberg | Super-resolution through error energy reduction[END_REF][START_REF] Papoulis | A New Algorithm in Spectral Analysis and Band-Limited Extrapolation[END_REF][START_REF] Jain | Extrapolation algorithm for discrete signal with application in spectral estimation[END_REF][START_REF] Strohmer | On discrete band-limited signal extrapolation[END_REF], the use of CS techniques allows to handle cases which would be intractable otherwise. In [START_REF] Donoho | Uncertainty principles and signal recovery[END_REF], a bound is given for uniqueness and perfect recovery using GP; within our notation, it reads:

N > d 1 - 1 2s .
This is highly more restrictive than the restriction N ≥ O(s log d) which suffices to apply the theorems in [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies[END_REF] holding for ℓ 1 minimization and random Fourier measurements.

Usual random sequences possibly contain the abscissas to be extrapolated, which is a problem because they correspond to unknown values of the signal (these values we precisely want to discover through the extrapolation process). Hence it may be tempting to restrict random sequences to the set of abscissas corresponding to the observed part only, which is the most straightforward way to proceed. However, such a choice leads to numerical results showing a sharp transition in terms of absolute errors when the edge of the observations interval is crossed. It resulted more efficient to generate global random sequences and, as a post-processing, to discard the subset of which perceiving something inside the unknown part. This seemingly differs in terms of numerical accuracy because the transition zone between observed and extrapolated parts of the signal is less noticeable; see for instance the figure 10 for which the signal is extrapolated in the area t ∈]1, 1.3].

Finally, a simple extrapolation algorithm based on Compressed Sensing techniques reads:

• Evaluate the sparsity s of the Fourier transform of the signal x ∈ R d to be reconstructed, for which there are d obs known values and d -d obs values to be extrapolated (generally located in one or two connected sets). • Compute the length N = O(s log d) of the random sequence to be used to under-sample

x in order to produce y; decreasing the ratio α = d/d obs in case no convenient sequence can be found may be necessary. • Build up the sensing matrix Φ = 1 √ N (exp(i(k -1)(j -1)/d)) k random of length N, j=1,...,d and y = (x k ) k random of length N .

• Compute the optimal x opt by ℓ 1 minimization or x romp through the ROMP algorithm.

• In general, one has to multiply this vector by √ N in order to produce the extrapolated signal by carrying out its FFT afterward.

Numerical results

The simple test-case

We repeat the simple test-case of producing a 15% (i.e. α = 1.15) extrapolation of the function [START_REF] Gosse | Analysis and short-time extrapolation of stock market indexes through projection onto discrete wavelet subspaces[END_REF] observed in the interval [-1, 1] using the ROMP algorithm. The sparsity level has been fixed as s = 8 and a Fourier measurement matrix was used with N the integer part of s log(K) where K = 256 is the number of discretization points in [-1, 1]. Results are shown in Fig. 8 which can be directly compared to Fig. 5 interval since the absolute errors in this region are of the order of 10 -3 ; the absolute errors in the extrapolation region are around 10 -2 . We stress that the Fourier basis we used isn't optimal for such a signal which isn't periodic; even if the difference |f pro (1.3) -f pro (-1)| isn't big, the numerical Fourier transform won't be rigorously sparse. Hence we believe that part of the extrapolation error comes from this fact because the ROMP algorithm seeks an extrapolation endowed with a s-sparse representation in the Fourier basis. This phenomenon should decrease when considering signals of higher frequency since the accuracy of Compressed Sensing doesn't depend on it (on the contrary of Prolate functions which heavily depends on the bandwidth). This experiment confirms our choice of extrapolating lower frequencies by means of Prolate functions; compare Fig. 8 with Fig. 5.

Remark 5 There is a special case for which the Fourier basis is especially well adapted, namely the one involving signals being periodic with a period equal to the Lebesgue measure of the whole interval (in our context |(-1, 1.3)| = 2.3). For instance, we briefly consider f per (t) = cos(8π(t + 1)/2.3) + sin(14π(t + 1)/2.3); from observations of f inside the interval [-1, 1], we aim at extrapolating in ]1, 1.3] using the ROMP algorithm with N = 1.5s log d random measures where s = 6, d = 295. Many convenient random sequences exist and the numerical results are displayed in Fig. 9. The salient feature is that for this particular class of signals, the transition between the absolute errors inside the observations interval and the ones in the extrapolation domain is hardly noticeable. This particular case shows how crucial is the choice of the basis for the quality of recovery and extrapolation.

First composite problem

We now deal with the oscillatory part of (15) which reads:

∀t ∈ [-1, 1], f hi (t) = 0.15[sin(13πt) -0.5 cos(23πt) + 0.3 cos(7πt)].
We carry out a 15% extrapolation of this strongly band-limited function by means of the In particular, finding convenient random sequences as explained in §4.3 is possible only for N < 1.2s log(d) ≃ 82. The results displayed in Fig. 10 show that these CS techniques handle such an extrapolation problem quite well; the absolute recovery errors are spread rather uniformly inside the interval [-1, 1.3] with a slight increase close to its edges. When comparing with Fig. 8, one sees that the Fourier basis isn't really adapted to the present case as well; however, the results are still rather good because the salient feature of the signal is its frequencies. And since there's no limitation in Compressed Sensing related to the bandwidth (only the sparsity of the Fourier transform matters), higher frequencies allows to get a slightly oversampled signal as an input.

Second composite problem

The oscillatory part of ( 16) has lower frequencies involved and reads:

∀t ∈ [-1, 1], f hi (t) = 1 2 [sin(5πt) -cos(7πt)].
We treat its extrapolation from [-1, 1] to ]1, 1.3] by means of ROMP with s = 7 and a slightly lower number of random observations, namely N = s log(295) ≃ 40. Results are shown in Fig. 11. This case contains frequencies lower than the preceding one, but its structure is a bit simpler as only 2 trigonometric functions are involved. Fig. 11 shows that the reconstruction is rather good uniformly in the computational domain, but there is a slight increase of absolute errors close to the edges of the interval [-1, 1] of observations. Starting from t ≃ 0.8, there are some small high-frequencies oscillations which don't perturb very much the overall quality of the extrapolated signal.

Putting the pieces altogether

By doing a byproduct of Theorem 3.11 in [START_REF] Shkolniskya | Approximation of bandlimited functions[END_REF] (see also Theorem 3 in [START_REF] Boyd | Approximation of an analytic function on a finite real interval by a bandlimited function and conjectures on properties of prolate spheroidal functions[END_REF]) and Theorem 7, it would be possible to derive an error estimate for the recovery of the original signal f by means of both algorithms put together inside the observations interval only. Indeed, it seems that no rigorous error estimates exist concerning the extrapolation of a band-limited signal through a truncated Prolate series of the type (8), even if in our context, ( 17) and ( 18) provide a reliable indication of the discrepancy coming from the Compressed Sensing part. This may constitute a direction for future research.

First composite problem

We now present the results of combined extrapolation algorithms on problem [START_REF] Gubner | A Simple Method for Computing Projections onto Subspaces of Prolate Spheroidal Wave Functions[END_REF] that the overall quality is quite good from both viewpoints of recovery (inside the interval [-1, 1]) and extrapolation (inside ]1, 1.3]). Absolute errors are uniformly distributed inside the whole interval [-1, 1.3] and their amplitudes vary roughly between 10 -3 and 10 -1 ; especially, there is no sharp transition when passing the edge t = 1. From Fig. 6, we know that absolute errors inside [-1, 1] come in very large part from the ROMP algorithm since the reconstruction errors for Prolate functions are negligible. In [-1, 1.3], both algorithms are endowed with roughly the same amount of numerical errors, and they accumulate when adding lower and higher parts of the signal. However, the numerical outcome can be considered as globally satisfying.

Second composite problem

The extrapolation of function ( 16) is shown in Fig. 13. This test-case is slightly more difficult because the correlation between f lo and f hi is stronger, which means that the difference between the lower and higher frequencies parts in the signal f is less sharp. As a consequence, it is less straightforward to make a clean separation between its "trend" and its "fluctuation" (compare with Figures 6,10 and 12). There are spurious oscillations of high frequency around the location t = 0.5 and this comes from the recovery of ROMP since the reconstruction errors of Prolate functions are once again negligible (see Fig. 7). Absolute errors oscillate between 10 -2 and 10 -1 which amounts to a relative error around 10% except in this instability area where it comes close to 1 2 . It is quite satisfying that, similarly as the preceding example, no visible transition appears when crossing the abscissa t = 1 separating measures from extrapolation.

Conclusion and outlook

We have presented in this paper a practical methodology to carry out the extrapolation of smooth and compressible signals relying on modern numerical techniques, namely the computation of Slepian's Prolate functions which are well suited for processing the reasonably low frequency signals and the Compressive Sampling techniques recently developed by Candès et al. which are independent of the size of the bandwidth hence are endowed with enough robustness to treat higher frequencies, especially when they can be expressed by means of a sparse Fourier representation. Each method has been tested independently and then, the absolute errors obtained by gluing together both algorithms have been studied. It came out that in order for this numerical method to deliver good results, several restrictions should be met: the correlation (1) between between low and high frequencies should be weak, the high frequencies should have a sparse representation in the Fourier basis (see §2.2.1 and Definition 2), and the size of the extrapolation domain should be reasonable (see §4.3). When these three points are met, numerical results are far better than those which could be obtained by standard methods, like e.g. Gerchberg-Papoulis methods at least because they wouldn't be able to treat such high frequencies (see Remark 4).

Fig. 1 .

 1 Fig. 1. Slepian functions (left) with c = 7 and repartition of eigenvalues λ n (right).

Fig. 2 .

 2 Fig. 2. Slepian functions (left) with c = 21 and repartition of eigenvalues λ n (right).

Fig. 3 .

 3 Fig. 3. Sparse spectrum recovery (left), absolute errors (middle) and ℓ 2 minimization (right).

Fig. 4 .

 4 Fig. 4. Nodes (left) and weights (right) for Gauss-Legendre quadrature and K = 256.

Fig. 5 .

 5 Fig. 5. Extrapolation (left), absolute errors (middle) and Prolate coefficients (right).

Fig. 6 .

 6 Fig. 6. Extrapolation (left), absolute errors (middle) and Prolate coefficients (right).

)Fig. 7 .

 7 Fig. 7. Extrapolation (left), absolute errors (middle) and Prolate coefficients (right).

Fig. 8 .

 8 Fig. 8. Extrapolation with ROMP (left) and absolute errors (right).

Fig. 9 .Fig. 10 .

 910 Fig. 9. Periodic extrapolation with ROMP (left) and absolute errors (right).

Fig. 11 .

 11 Fig. 11. Extrapolation with ROMP (left) and absolute errors (right).

Fig. 12 .

 12 Fig. 12. First case: global extrapolation (left) and absolute errors (right).

Fig. 13 .

 13 Fig. 13. Second case: global extrapolation (left) and absolute errors (right).

  1.3. §2.2 is devoted to surveying some important results in the theory of Compressed Sensing; especially the so-called Restricted Isometry Property, a crucial feature of certain random matrices, is recalled in §2.2.1 and recovery theorems are stated in §2.2.2. §2.2.3 shows a numerical illustration of the perfect recovery phenomenon for a sparse signal. In §3, we explain in detail how to set up a practical algorithm to conduct extrapolation of any band-limited (preferably low frequency) signal relying on the projection of the Paley-Wiener subspace of L 2 (R) generated by Prolate functions. Coefficients are computed using Gauss-Legendre quadrature. Several numerical test-cases are shown in §3.3. Extrapolation of higher frequencies relying on Compressed Sensing techniques are dealt with in §4. In §4.1, we introduce the useful class of "compressible signals" and state recovery results for them. In §4.2, we present a simple but efficient iterative recovery algorithm called ROMP endowed with good stability properties. The issue of finding convenient random sequences of observations is discussed in §4.3 where also a comparison between CS and standard Gerchberg-Papoulis recovery methods is made (seeRemark 4). Numerical test-cases are presented in §4. In §5, we display the numerical outcome of an overall extrapolation process where the signal to be treated is first split between a "trend" and a "fluctuation" being processed by Prolate functions and Compressed Sensing respectively; two test-cases are considered. Results are to be checked on Figures12 and 13. Conclusive remarks are given in §6.

one could even consider f in a space of tempered distributions S ′ on R, like for instance the Zakaï class[START_REF] Zakai | Band-Limited Functions and the Sampling Theorem Information[END_REF][START_REF] Lee | Approximate Interpolation and the Sampling Theorem[END_REF][START_REF] Cambanis | Zakai's class of bandlimited functions and processes: its characterization and properties[END_REF] but this is out of the scope of the present paper.