Laurent Condat

Reconstruction from non-uniform samples: A direct, variational approach in shift-invariant spaces

Keywords: Non-uniform sampling, variational reconstruction, interpolation, shift-invariant spaces, splines

We propose a new approach for signal reconstruction from non-uniform samples, without any constraint on their locations. We look for a function that minimizes a classical regularized least-squares criterion, but with the additional constraint that the solution lies in a chosen linear shift-invariant space-typically, a spline space. In comparison with a pure variational treatment involving radial basis functions, our approach is resolution dependent; an important feature for many applications. Moreover, the solution can be computed exactly by a fast non-iterative algorithm, that exploits at best the particular structure of the problem.

Introduction

In modern digital data processing and communication systems, signals and numerical data are usually available as a sequence of discrete-eventually nonuniform and/or noisy-samples. For the purpose of deriving numerical algorithms, it is sometimes desirable to represent the data by a continuously-defined parametric function. This is particularly relevant for resampling tasks: a function is fitted to the data and resampled at new locations. Non-integer sampling rate conversion, arbitrary time delay, or edge detection in images are examples where such a treatment may be required.

Numerous methods have been proposed in the literature for reconstructing a signal from non-uniform samples. Since the pioneering work of Shannon [START_REF] Shannon | Communication in the presence of noise[END_REF], it is known that any π/T -bandlimited function s(t) ∈ L 2 (R) can be perfectly reconstructed from its uniform samples s(T n), n ∈ Z. This is still true when considering non-uniform samples, assuming strong limitations on the samples locations [START_REF] Gröchenig | On Landau's necessary density conditions for sampling and interpolation of band-limited functions[END_REF]. Perfect reconstruction from non-uniform samples is also possible under the weaker assumption that s(t) belongs to a linear shift-invariant (LSI) space [START_REF] Chen | A sampling theorem for shift-invariant subspaces[END_REF], with similar constraints on the sampling set [START_REF] Aldroubi | Beurling-Landau-type theorems for nonuniform sampling in shift invariant spline spaces[END_REF]. Practical iterative algorithms have been proposed for reconstruction in the bandlimited case [START_REF] Feichtinger | Efficient numerical methods in non-uniform sampling theory[END_REF][START_REF] Strohmer | Numerical analysis of the non-uniform sampling problem[END_REF] and in LSI spaces [START_REF] Liu | Irregular sampling for spline wavelet subspaces[END_REF][START_REF] Aldroubi | Nonuniform sampling and reconstruction in shift-invariant spaces[END_REF]. Strong conditions have to be met for these algorithms to converge [START_REF] Chen | Maximal gap of a sampling set for the exact iterative reconstruction algorithm in shift-invariant spaces[END_REF].

In this article, we propose a novel approach for one-dimensional signal reconstruction from non-uniform samples, without any constraint on their locations. We aim at reconstructing a continuous-time function f T (t) that best fits the data, while being resolution dependent, i.e. completely determined by its values f T (T k), k ∈ Z, where the parameter T can be chosen arbitrarily. In other words, we seek a reconstructed function that is constrained to lie in a linear shift-invariant space. This ensures that the solution is parameterized by coefficients attached to the uniform reconstruction grid having the desired resolution, and not to the data locations. In this reconstruction space, we formulate reconstruction as a regularized least-squares problem; that is, the solution minimizes a cost depending on two terms: the sum of the squared errors at the sampling locations on one side, and a quadratic variational functional that enforces the solution to be smooth on the other side. This formulation circumvents the strong limitations-typically, a restriction on the maximal gap between samples-of formulations where only the fit to the data is considered. The main advantage of our formulation is that the coefficients, that parameterize the reconstructed function, can be computed using a simple, fast, and non-iterative algorithm. In essence, it performs a two-pass time-varying recursive filtering of the data.

The paper is organized as follows. We first give the necessary mathematical background and formulate the problem in Section 2. In Section 3, we derive its solution and discuss the novelty of our approach in comparison with previous related works. We then propose a practical algorithm for computing the solution in Section 4. In Section 5, we discuss the properties of the reconstructed function and illustrate our method with experimental results. In the last section, some applications are discussed.

Mathematical preliminaries 2.1 Definitions and notations

Throughout the paper, parentheses are used for continuous-time signals, e.g., f (t), and brackets for the samples of a discrete time signal like u = (u[n]) n∈Z . Continuous and discrete convolutions are denoted by * . We define the z-transform of a discrete signal u by U (z) = n∈Z u[n]z -n , and the convolution inverse u -1 of u as the sequence whose z-transform is 1/U (z). We also introduce the autocorrelation function a ϕ of a function ϕ: a ϕ (t) = ( φ * ϕ) (t), using the flip operator φ(t) = ϕ(-t).

Linear shift-invariant spaces

In this work, we make use of functions that belong to linear shift-invariant (LSI) spaces. Let the real number T > 0 denote our choice of sampling step. The LSI space V T (ϕ) is the functional space spanned by the T -shifts of a generating function ϕ( t T ):

V T (ϕ) = k∈Z c T [k] ϕ( t T -k) : c T [k] ∈ R ∀k ∈ Z . (1)
The potential of LSI spaces has been recognized for quite some time. They have been used in finite elements and approximation theory [START_REF] Boor | The structure of finitely generated shift-invariant spaces in L 2 (R d )[END_REF][START_REF] Jia | Shift-invariant spaces and linear operator equations[END_REF] and for the construction of multiresolution approximations and wavelets [START_REF] Daubechies | Ten lectures on wavelets[END_REF]. In sampling theory, they are used extensively [START_REF] Unser | Sampling-50 Years after Shannon[END_REF][START_REF] Aldroubi | Nonuniform sampling and reconstruction in shift-invariant spaces[END_REF]. We assume in this work that ϕ is bounded with compact support, and that the functions {ϕ

( t T -k)} form a Riesz basis of V T (ϕ)∩L 2 (R), which ensures that each function f T ∈ V T (ϕ) has a unique expansion of the form f T (t) = k∈Z c T [k] ϕ( t T -k).
This last condition is equivalent to the requirement that there exist two constants 0 < C 1 and C 2 < +∞ (the lower and upper Riesz bounds) such that C 1 < k∈Z a ϕ (k)e -jωk < C 2 a.e.

A classical example of LSI space is the set of π/T -bandlimited functions, obtained with ϕ(t) = sinc(t) [START_REF] Shannon | Communication in the presence of noise[END_REF]. Particularly useful LSI spaces are the polynomial spline spaces [START_REF] De Boor | A Practical Guide to Splines[END_REF][START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF]. They are obtained by choosing ϕ = β d , the centered B-spline of degree d, that is symmetric and has compact support of width d + 1. The B-splines are constructed by successive convolution:

β d = β 0 * β d-1 , from the indicator function β 0 defined by β 0 (t) = 1 if t ∈ [-1 2 , 1
2 ), 0 otherwise .

3 Variational reconstruction

Problem statement

We assume that we are given a finite number

N of measurements (v[n]) n∈[0,N -1] at locations (x[n]) n∈[0,N -1] within an interval I = [0, S]; that is, v[n] = s(x[n])
where s(t) is some unknown process defined on I. Let us choose a resolution step T > 0 (whose choice will be discussed in Sections 5, 6). Our goal is to reconstruct a continuous-time function f T (t) defined on I, that modelizes the data and belongs to the reconstruction space V T (ϕ). That is, we look for a function having the form

f T (t) = k∈Z c T [k]ϕ( t T -k), (2)
where the discrete coefficients c T [k] are the unknowns to be determined. We assume that the generator ϕ has compact support, included in the interval (-W, W ), with W ∈ N, and that S = KT , for some K ∈ N. Then, for every Section 3.3). This allows to use simple notations indexed by k ∈ Z, but we have to keep in mind that only a finite number of coefficients are not zero, e.g. in [START_REF] Gröchenig | On Landau's necessary density conditions for sampling and interpolation of band-limited functions[END_REF].

t 0 ∈ I the value f T (t 0 ) is determined by only a few coefficients c T [k], for k in the interval ( t0 T -W, t0 T + W ). Therefore, the function f T is completely determined on the interval I by the c T [k], k ∈ [-W + 1, K + W -1]. So, for convenience, we only consider functions of V T (ϕ) such that c T [k] = 0 for every k / ∈ [-W + 1, K + W -1] in (1) (except in
We define our reconstruction problem as a variational problem, whose f T is the solution:

f T = argmin g∈VT (ϕ) N -1 n=0 g(x[n]) -v[n] 2 + λ I |g (r) (t)| 2 dt . (3)
where g (r) is the r th derivative of g, for some integer r ≥ 1. This variational criterion is composed of two antagonist terms, one controlling the closeness to the data, the other one enforcing the solution to be smooth. Interestingly, these terms are similar to the so-called external forces (the function is attracted by the data) and internal forces (the bending energy of the function tends to be minimized) in the theory of snakes and active contours used in computer graphics [START_REF] Jacob | Efficient energies and algorithms for parametric snakes[END_REF]. The quantity λ > 0 is a Lagrangian parameter working as a tradeoff factor between these two terms. The integer r controls the smoothness of the reconstruction. The values r = 1 and r = 2 are the most frequently used, and correspond to searching a function that has maximum flatness, and minimum curvature, respectively. The resolution parameter T > 0 controls the coarseness of the representation. Finally, the generator function ϕ(t) controls the space where the reconstructed function lives. Signal processing practitioners often rely on localized kernels such as B-splines [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF] or MOMS [START_REF] Blu | Moms: Maximal-order interpolation of minimal support[END_REF]. Other compactly-supported kernels were recently introduced in [START_REF] Madani | Optimized compactsupport interpolation kernels[END_REF]. The influence of this set of parameters is discussed and illustrated in Section 5.

Before determining the solution of (3) in the remainder of Sect. 3, we have to introduce some mathematical safeguards for the problem to be well posed. First, we suppose that ϕ is such that R |ϕ (r) (t)| 2 dt < ∞. For instance, if choosing a B-spline ϕ = β d , then d ≥ r is necessary. Moreover, we require that the samples are at r or more distinct locations. We will see that these conditions allow the solution of our problem to exist and to be unique.

It is important to note that the Tikhonov-like formulation (3) is not limited to the choice of the squared norm of the derivative for the regularization term. The approach could be extended without difficulty to any other quadratic shiftinvariant penalty with impulse response of compact support. We exemplify the approach with the energy of the r-th derivative because it is a classical penalty for functions living in Sobolev spaces, but another regularizer may be more appropriate, for a given application. The reader may wonder if a non-quadratic criterion like the total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] could be considered; that is, replacing the integral in (3) by I |g (r) (t)|dt (without the square). In that case, an iterative algorithm for convex non-smooth optimization would be necessary to solve (3), whereas the efficiency of the proposed non-iterative algorithm is due to the linear structure of the problem. Therefore, the extension of the approach to non-quadratic regularization is beyond the scope of the paper and is left for future work.

Solution to (3)

Finding the reconstructed function f T (t) amounts to determining the sequence c T in (2), so that the cost function

Ψ(c T ) = N -1 n=0 f T (x[n]) -v[n] 2 + λ I |f (r) T (t)| 2 dt (4)
is minimized. First, we rewrite the data fidelity term as a function of c T :

N -1 n=0 f T (x[n]) -v[n] 2 = N -1 n=0 v[n] - k∈Z c T [k]ϕ( x[n] T -k) 2 . (5)
Second, we rewrite the variational term as a function of c T , with the integral over R and not I in a first time. Let us introduce the discrete sequence q ϕ,r defined by q ϕ,r [k] = (-1) r T a (2r) ϕ (k). Since ϕ (r) (-t) = (-1) r φ(r) (t), and differentiations commute with convolutions, we have:

R f (r) T (t) 2 dt = 1 T 2 R k∈Z c T [k]ϕ (r) ( t T -k) 2 dt = 1 T k,l∈Z c T [k]c T [l] R ϕ (r) (x -(k -l))ϕ (r) (x) dx = (-1) r T k,l∈Z c T [k]c T [l] ( φ(r) * ϕ (r) ) (k -l) = (-1) r T k,l∈Z c T [k]c T [l] a (2r) ϕ (k -l) = k∈Z c T [k] (c T * q ϕ,r )[k]. (6)
In the case of spline reconstruction (ϕ = β d ), we can give the general form of the filter q ϕ,r . First, B-splines verify the simple relation a β d (t) = β 2d+1 (t). Moreover, the derivative of a spline is also a spline of lower degree, since

β d (1) (t) = β d-1 (t + 1 2 ) -β d-1 (t - 1
(z) = 1, B 1 (z) = 1, B 2 (z) = 1 8 z + 3 4 + 1 8 z -1 , B 3 (z) = 1 6 z + 2 3 + 1 6 z -1 .
Then the filter q ϕ,r has the following form:

Q β d ,r (z) = 1 T (-z + 2 -z -1 ) r B 2d+1-2r (z). (7)
For example, if ϕ = β 1 and r = 1, then

Q ϕ,r (z) = (-z + 2 -z -1 )/T . If ϕ = β 3 and r = 2, then Q ϕ,r (z) = (z 3 -9z + 16 -9z -1 + z -3 )/6T .
In order to express the cost Ψ(c T ) in terms of matrices and vectors, we introduce the following quantities (• T is the transpose operator):

c = c T [-W + 1] • • • c T [K + W -1] T , v = v[0] v[1] • • • v[N -1] T , M = [(M [n, k]) n,k ] with M [n, k] = ϕ( x[n] T -k), (8)
for

n ∈ [0, N -1], k ∈ [-W + 1, K + W -1], Q = [(Q[k, l]) k,l ] for k, l ∈ [-W + 1, K + W -1]. (9)
Note that performing a matrix-vector product like Mc is equivalent to applying a time-varying filter to the signal c T . Using these matrices, the cost Ψ(c T ) can be rewritten as:

Ψ(c) = Mc -v 2 + λc T Qc, (10)
where the values of the matrix Q are given by ( 6):

Q[k, l] = q ϕ,r [k -l] = (-1) r T a (2r) ϕ (k -l), (11)
except for the first and last rows of Q that contain particular values because Ψ(c T ) is defined with the integral over I, and not R as in [START_REF] Strohmer | Numerical analysis of the non-uniform sampling problem[END_REF]. These special values are in squares of size (2W -1) 2 in the lower-left and upper-right corners of the matrix. To compute them for the left boundary (this would be the same for the other boundary), we have to develop the left-hand side of the following equality, and identify the coefficients with its right-hand side:

2W -1 0 W -1 k=-W +1 c T [k]ϕ (r) ( t T -k) 2 dt = W -1 k,l=-W +1 Q[k, l]c T [k]c T [l]. (12)
For instance, if ϕ = β 1 , r = 1 and ϕ = β 3 , r = 2, Q takes the respective forms:

Q = 1 T × Q = 1 6T ×           1 -1 0 0 0 • • • -1 2 -1 0 0 . . . 0 -1 2 -1 0 . . . 0 0 -1 2 -1 . . . 0 0 0 -1 2 . . . . . . . . . . . . . . . . . . . . .           ,           2 -3 0 1 0 • • • -3 8 -6 0 1 . . .
When considering the rows of this system, we obtain a set of equations that can be found directly from ( 5) and ( 6) by setting the partial derivatives ∂Ψ/∂c

T [k] to zero, for each k ∈ [-W + 1, K + W -1]: l∈Z N -1 n=0 ϕ( x[n] T -k)ϕ( x[n] T -l) + λ Q[k, l] c T [l] = N -1 n=0 ϕ( x[n] T -k)v[n]. (14)
We have to check that the solution to the linear system ( 13) is well defined. In fact, M T M and Q are positive semi-definite matrices, and A is symmetric and positive definite, as proved in Section 8. As a consequence, the linear system has a unique solution. Moreover, thanks to the hypothesis of compact support of ϕ, A is band-diagonal with only 4W -1 diagonals containing non-zero entries; more precisely, A[k, l] = 0 if |k -l| ≥ 2W . This will be the key for efficiently solving this linear system, as detailed in Section 4.

Reconstruction with mirror boundary conditions

There is a second possibility to deal with reconstruction on a finite interval: considering mirror boundary conditions that allow to extend a finite signal to an infinite one by symmetry and periodicity [START_REF] Brislawn | Classification of nonexpansive symmetric extension transforms for multirate filter banks[END_REF]. The main advantage is that there is no overhead of coefficients as with the previous method; that is, we only have to compute c T [k] for k ∈ [0, K]: the decomposition is said to be nonexpansive. In this subsection, we assume that ϕ is symmetric (ϕ(t) = ϕ(-t)).

The proposed method consists in solving the following problem:

f T = argmin g∈VT (ϕ)∩S N -1 n=0 g(x[n]) -v[n] 2 + λ I |g (r) (t)| 2 dt . ( 15
)
where S is the set of functions that are symmetric around t = 0 and t = T K. With these boundary conditions, the coefficients c T of f T verify the following symmetry and periodicity conditions:

c T [-k] = c T [k] and c T [k + 2pK] = c T [k] for every k, p ∈ Z. In other words, f T is completely determined by the coeffi- cients c T [k], k ∈ [0, K].
An interpretation of ( 15) is to consider that we want to fit a function on an infinite data set that is symmetric and periodic, as illustrated in Fig. 1.

The coefficients (c T [k]

) k∈[0,K] are solution of a linear system A c = y to be determined (the tildes are for distinguishing the matrices from the ones previously defined). Let us study the structure of this system on the left boundary

= c T [W -1] • • • c T [1] c T [0] c T [1] • • • T because of the boundary conditions c T [-k] = c T [k]. Each row of Mc can be expanded as • • • + M [k, -1]c T [-1] + M [k, 0]c T [0] + M [k, 1]c T [1] + • • • , that can be rewritten M [k, 0]c T [0] + (M [k, -1] + M [k, 1])c T [1] + • • • . Therefore, M c = Mc,
where M is derived from M by removing its first W -1 columns and adding their values to the ones of the mirror columns with respect to the W th column (with the similar manipulation at the right boundary). This is equivalent to "folding" the matrix M, as in the following example with W = 2, where the first column is added to the third one and removed:

M =   a b c d • • • e f g h • • • . . . . . . . . . . . . . . .   -→ M =   b c + a d . . . f g + e h . . . . . . . . . . . . . . .   . ( 16
)
This folding operation has to be performed on the first and last rows and columns of the matrix Q and A in order to get Q and A, respectively. Similarly, the first and last values of y have to be folded. For instance, if ϕ = β 3 and r = 2 (thus W = 2), this yields:

Q = 1 6T × Q = 1 6T ×           2 -3 0 1 0 • • • -3 8 -6 0 1 . . . 0
          -→          8 -9 0 1 . . . - 9
         . ( 17
)

Related works

Most of the works in the literature address the non-uniform reconstruction problem with restrictions on the samples locations, usually a density constraint or a maximum gap x[n + 1] -x[n], so as to achieve perfect reconstruction of an unknown signal s(t) from its non-uniform noiseless samples [START_REF] Gröchenig | On Landau's necessary density conditions for sampling and interpolation of band-limited functions[END_REF][START_REF] Aldroubi | Beurling-Landau-type theorems for nonuniform sampling in shift invariant spline spaces[END_REF][START_REF] Gröchenig | Fast local reconstruction methods for nonuniform sampling in shift invariant spaces[END_REF]. If noise is present, some of these methods yield a reconstructed function that approaches the data at best in the least-squares sense. The variational formulation adopted in the present paper allows to handle the problem in whole generality, whatever the distribution of samples. It is classical, in approximation theory and statistics, to reconstruct a function f (t) that minimizes the criterion of ( 3), but over the whole Sobolev space of order r instead of over V T (ϕ) [START_REF] Eubank | Nonparametric regression and spline smoothing[END_REF]. The solution to this problem can be expressed as [START_REF] De Boor | A Practical Guide to Splines[END_REF]. It is made of a polynomial p(t) of degree less than r and a linear combination of radial basis functions (RBF) |t| 2r-1 positioned at the sampling locations x[n]. This implies that the solution f (t) is a non-uniform polynomial spline of degree 2r -1 with knots at the x[n]. These spline spaces have other basis functions, the non-uniform B-splines, which are compactly supported and, therefore, much better conditioned than the radial basis functions [START_REF] De Boor | A Practical Guide to Splines[END_REF]. Non-uniform B-splines revert to classical B-splines when the samples are uniform. In comparison with the RBF framework, our approach enforces the solution to lie in V T (ϕ), and thus to be resolution-dependent. Our solution can be expressed as a linear combination of translates of a single generator ϕ, while in the RBF framework, a non-uniform B-spline has to be determined for each sample. Moreover, the choice of ϕ is free in our approach, while in the RBF framework, it is induced by the differential regularization operator, and thus constrained to be a spline of odd degree.

f (t) = N -1 n=0 c[n] |t -x[n]| 2r-1 + p(t)
The idea of minimizing a regularized least-squares criterion in a LSI space has also been proposed recently in 2-D for image reconstruction [START_REF] Arigovindan | Variational image reconstruction from arbitrarily spaced samples: A fast multiresolution spline solution[END_REF][START_REF] Vázquez | Reconstruction of nonuniformly sampled images in spline spaces[END_REF], but the associated algorithms, although relatively fast in comparison with previous works, are still iterative and not adapted to real-time implementations. In fact, our algorithm is limited to the 1-D case. This is the price to pay for having the direct and very efficient implementation based on time-varying recursive filtering, detailed in Section 4.

Another related approach has been proposed in [START_REF] Barrutia | Non-uniform to uniform grid conversion using least-squares splines[END_REF]; it also provides a resolution-dependent reconstruction f T (t), in the noise-free case. To this end, the non-uniform spline f (t) interpolating the data is computed, and projected orthogonally onto V T (ϕ) in the L 2 sense, so as to obtain f T (t). This method is more complicated than ours, since the non-uniform interpolation has to be performed first, and the projection step is relatively involved. Moreover, this formulation gives more credit, in the definition of f T (t), to the estimate f (t) than to the data (v[n]). It is more reliable to directly formulate the solution f T (t) from the available data (v[n]), as proposed in this work.

A fast reconstruction algorithm 4.1 Strategy

We have seen that our problem boils down to solving a linear system. It is possible to perform this operation using standard linear algebra libraries (such as LAPACK), that implement routines to factorize matrices in various ways and to solve linear systems. Such libraries are highly optimized, but they are generic. The practitioner interested by our approach may want to have a dedicated algorithm, that exploits at best the particular structure of the problem. So, we detail in this section the best strategy for solving the linear system Ac = y. This system is symmetric, positive definite, and band-diagonal. Thus, the best way to solve it relies on the Cholesky factorization of A and consists in the following three steps:

1. The Cholesky decomposition of A is performed; that is, we compute the unique lower triangular matrix L such that A = LL T [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF][START_REF] Golub | Matrix computations[END_REF]. L is also band-diagonal, with 2W diagonals containing non-zero entries:

L[k, l] = 0 only if -2W < k -l ≤ 0.
The Cholesky decomposition can be performed row by row by a fast version of Crout's algorithm [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] that takes advantage of the band-diagonal structure of A.

The lower triangular system L

• c = y is solved by forward substitution: for

k from k min to k max , • c T [k] = 1 L[k, k] y[k] - -1 i=-2W +1 L[k, k + i] • c T [k + i] , (18)
3. The upper triangular system L T c =

• c is finally solved by backward substitution: for k from k max down to k min ,

c T [k] = 1 L[k, k] • c T [k] - 2W -1 i=1 L[k + i, k]c T [k + i] , (19)

Practical algorithm

The practical algorithm that computes the coefficients c

T [k], k ∈ [-W + 1, K + W -1]
, consists in two passes: the first pass implements the Cholesky decomposition and solves the first system [START_REF] Madani | Optimized compactsupport interpolation kernels[END_REF], while the second pass solves the second linear system [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. We define the auxiliary variables a

[i] = A[k, k + i] and u[k, i] = L[k + i, k].
Then our algorithm can be written in pseudo-code as:

• First pass:

for k from -W + 1 to K + W -1 { i min := max(-2W + 1, -W + 1 -k); i max := min(2W -1, K + W -1 -k); for i from 0 to i max , a[i] := N -1 n=0 ϕ( x[n] T -k)ϕ( x[n] T -k -i) + λ Q[k, k + i]; u[k, 0] := a[0] - -1 i=imin u[k + i, -i] 2 1/2 ; • c T [k] := 1 u[k, 0] N -1 n=0 ϕ( x[n] T -k)v[n] - -1 i=imin u[k + i, -i] • c T [k + i] ; for i from 1 to i max , u[k, i] := 1 u[k, 0] a[i] - -1 j=max(i-2W +1,-W +1-k) u[k + j, -j] u[k + j, i -j] ; } • Second pass: for k from K + W -1 down to -W + 1 { i max := min(2W -1, K + W -1 -k); c T [k] := 1 u[k, 0] • c T [k] - imax i=1 u[k, i]c T [k + i] ; }
In fact, the sums indexed by n have to be evaluated only for the few samples located in the interval (T (k -W ), T (k + W )) that have a non-zero contribution. If the samples are ordered such that x[n + 1] ≥ x[n] for every n, the data set (x[n], v[n]) can be accessed progressively. Hence the first pass can be computed on the fly, as the data are made available. For example, if x[n] is the time of the measurement v[n], the first pass can be performed in real-time with a delay of at most W time units. Moreover, during the second backward pass, the coefficients c T [k] can be computed in place, replacing the intermediate values

• c T [k]
. Note that this algorithm can be interpreted as a two-pass time-varying recursive filtering on the sequence y, with a filter updated on the fly by the Cholesky decomposition.

If the data locations are not sorted, or if N ≫ K, the computation time will be consumed mostly for the evaluation of the sums indexed by n. In this case, • First pass:

for k from -W + 1 to K + W -1 { • c T [k] = 0; for i from 0 to min(2W -1, K + W -1 -k), u[k, i] := λ Q[k, k + i]; } for n from 0 to N -1, for k from ⌊ x[n] T + 1 -W ⌋ to ⌈ x[n] T -1 + W ⌉ { for i from 0 to min(2W -1, K + W -1 -k), u[k, i] := u[k, i]+ϕ( x[n] T -k)ϕ( x[n] T -k-i); • c T [k] := • c T [k] + ϕ( x[n] T -k)v[n]; } • Second pass: for k from -W + 1 to K + W -1 { i min := max(-2W + 1, -W + 1 -k); i max := min(2W -1, K + W -1 -k); u[k, 0] := u[k, 0] - -1 i=imin u[k + i, -i] 2 1/2 ; • c T [k] := 1 u[k, 0] • c T [k]- -1 i=imin u[k + i, -i] • c T [k+i] ; for i from 1 to i max , u[k, i] := 1 u[k, 0] u[k, i] - -1 j=max(i-2W +1,-W +1-k) u[k + j, -j] u[k + j, i -j] ; } • Third pass: for k from K + W -1 down to -W + 1 { i max := min(2W -1, K + W -1 -k); c T [k] := 1 u[k, 0] • c T [k] - imax i=1 u[k, i]c T [k + i] ;
} These algorithms are for the strategy without boundary conditions, proposed in Section 4.2. The implementation for the method with mirror boundary conditions is very similar. In that case, the "folding" operation on the matrices is implemented by assigning each contribution ϕ(

x[n] T -k)ϕ( x[n] T -k -i) to its folded place directly, for example to u[-k, -i] instead of u[k, i] if k < 0.

Computation time and storage requirements

The computation time of the proposed algorithm can be modeled as O(W 2 N ), O(W 2 K), O(W N ) for the calculation of the elements in A, L and y respectively, and O(KW ) for the forward and backward substitutions. So, the total time reduces to O(W 2 (N + K)); it is linear in N and K, which is the best one could have hoped for. If the reconstruction is to be performed on an interval I of fixed size S = KT , the total time may be rewritten O(W 2 (N + S/T )), so as to let appear the linear dependence in 1/T . Experimental computation times are reported in Fig. 2 for an implementation in C language of the algorithm proposed in Section 4.2 (second variant), running on a 1.6 GHz laptop PC. The computation time is asymptotically linear in K and N , as predicted.

Apart of the memory required to store the coefficients c T [k], auxiliary memory of size 2W (K + 2W -1) units (or 2W (K + 1) if using mirror boundary conditions) is needed to store the values u[k, i], which are generated during the forward pass of our algorithm and used in the backward pass.

Instead of the Cholesky decomposition, it is possible to use a LU-factorization that does not exploit the symmetry of the matrix. Although this decomposition requires two times more computation, the square-root operator is not needed, and one diagonal is saved in the storage of L; that is, (2W -1)(K + 2W -1) instead of 2W (K + 2W -1) memory units are used. Another variant is the LDL T factorization, that requires one more pass on the data, but also avoids the square-root operator.

Note that our algorithm uses a Cholesky decomposition that is "extremely stable numerically" [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF]. However, even if our linear system is positive definite, its condition number clearly depends on the sampling locations. In fact, our method amounts to performing the deconvolution of a time-varying filter. In large gaps without samples, this filter reduces to q ϕ,r , which has roots on the complex unit circle. An inverse filter with poles on the unit circle is said to be marginally stable, because the impulse response corresponding to these poles does not decay, but does not grow either. So, if a round-off error occurs on a coefficient c T [k] during the computation, it can propagate to its neighbors inside a region without samples, but its amplitude, limited to the machine accuracy, will not grow. Therefore, this is not a problematic issue.

Choice of the parameters

In this section, we discuss the influence of the parameters ϕ, r, T , λ.

Since the reconstruction is performed in a LSI space, this space has to be chosen before hand. The "best" space depends on the characteristics of the signal to be modeled. If prior knowledge on the process that gave the samples is available, it can be used to choose a particular generator ϕ [START_REF] Ramani | Sampling in practice: Is the best reconstruction space bandlimited?[END_REF]. For instance, if the reconstructed function is required to be continuous and continuously differentiable, these properties will be enforced on ϕ. In concrete problems, spline spaces have shown to be particularly adequate for representing signals [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF]. Their optimal approximation properties have been demonstrated theoretically [START_REF] Blu | Quantitative Fourier analysis of approximation techniques: Part I-interpolators and projectors-and part II-wavelets[END_REF][START_REF] Blu | Moms: Maximal-order interpolation of minimal support[END_REF] and confirmed by practical experiments [START_REF] Lehmann | Survey: Interpolation methods in medical image processing[END_REF][START_REF] Thévenaz | Interpolation revisited[END_REF][START_REF] Meijering | Quantitative evaluation of convolution-based methods for medical image interpolation[END_REF]. Part of their interest lies in the fact that a B-spline has the maximal possible approximation order, given the size of its support. The ability to reproduce high order polynomials is of primary concern for approximating arbitrary signals. Besides, for a given parameter r, choosing ϕ = β 2r-1 ensures that our solution f T coincides with the RBF solution if the samples are uniform, at locations x[n] = T n. This means that spline functions are the natural choice when minimizing a criterion based on a derivative.

The parameter r controls the kind of smoothness that is enforced on the solution. Three configurations are illustrated with a simple synthetic example in Fig. 3. In (a), the choice ϕ = β 1 , r = 1 results in a piecewise linear recon- struction, with knots at the T k, k ∈ Z, which means that f T (t) is linear on each interval [T k, T (k + 1)]. Note that with ϕ = β 1 , there is no other possible choice than r = 1. Moreover, the two strategies for handling the reconstruction on a finite interval (no boundary conditions as in Section 4.2 or mirror boundary conditions as in Section 3.3) are equivalent in this case. In (b), (c), (d), ϕ = β 3 yields a smoother reconstruction, that is twice continuously differentiable. If r = 1, the variation |∇f T | 2 is minimized, and the solution tends to behave like a straight line in large gaps, as in (b). If r = 2, the second derivative modelizes a curvature energy, and in large gaps, the solution f T behaves like a polynomial of degree three, as in (c); in this case, f T can go beyond the dynamics of the initial samples, which may be a disadvantage. Note that, in large gaps, we do not have information about the signal to reconstruct. So, fitting a polynomial using the data at the boundaries of the gap is, in essence, the best we can do. The choice of the boundary conditions is illustrated in Fig. 3 (c) and (d): mirror boundary conditions yield a reconstruction whose first derivative is constrained to be zero at the boundaries of the reconstruction interval. In this case, f T is parameterized by 11 coefficients c T [k], k ∈ [0, 10]. In (d), no boundary conditions are enforced, but the solution is now parameterized by 13 coefficients

c T [k], k ∈ [-1, 11].
The parameter T controls the coarseness of the representation. When reconstructing a signal over an interval [0, S], we obtain a parametric solution with K + 1 = S/T + 1 degrees of freedom. If the parsimony of the representation modeling the data is an important criterion, for instance in coding applications or if the computation time is limited, then T will be chosen relatively large. Conversely, when T → 0, the solution f T becomes closer and closer to the nonuniform solution in the RBF framework, since the space V T (ϕ) becomes dense in the whole Sobolev space of order r. The influence of T is illustrated in Fig. 4 In practical applications, the parameter T will be matched to the cut-off frequency of the reconstruction lattice, as discussed in the next section. In fact, for each function f T ∈ V T (ϕ), there is a one-to-one correspondence between its coefficients c T [k] and its point values w

T [k] = f T (T k) at locations T k. That is why f T ∈ V T (ϕ) is said to have resolution 1/T .
The regularization factor λ is a key parameter: an excessive value will oversmooth the solution, while a small value will provide a solution that is close to the data, but may have large disturbing variations. Let us consider the behavior of the reconstructed function f T in the limit case when λ → 0. The solution c = (M T M + λQ) -1 M T s has a well-defined limit. This limit corresponds to setting exactly λ = 0 in the equations only if the matrix M T M is invertible. This is the case, for example, if 0 < x[n + 1] -x[n] < T for every n. In the general case, the limit solution has the following interpretation: when λ → 0, the regularization term in (3) becomes neglictible in comparison with the fit-todata term. So the limit function f T minimizes the criterion n g

(x[n])-v[n] 2 ,
and if the solution g(t) is not unique, the remaining degrees of freedom are used to minimize I |g (r) (t)| 2 . Therefore, the behavior of the reconstructed function is the following: if λ is large, f T is smooth, whatever the locations x[n] of the data, as in Fig. 5 (a). If λ is very small, there are two cases: if locally the samples set is sparse (i.e. x[n + 1] -x[n] > T for a few successive values of n), then f T almost passes through the samples, and is smooth inbetween (Fig. 5 (b), (c), (d), for t > 3). Thus large gaps are "in-painted" in a smooth way. Conversely, in a region where x[n + 1] -x[n] < T for a few successive values of n, f T approaches the data at best in the least-squares sense (Fig. 5 (b), (c), (d), for t < 2). Thus, in the noise-free case, it is tempting to choose λ very small. However, this may result in large unexpected oscillations, going far beyond the initial range of the signal, as shown in Fig. 5 (c), (d): the exact interpolation of the samples is a too strong constraint.

In practice, the measurements are often noisy or lacking precision and it is not suitable to reconstruct a function that interpolates the data exactly. It is therefore appropriate to achieve a tradeoff between the closeness of fit and the smoothness of the solution through the parameter λ. The best value for the problem at hand has to be determined on a case by case basis. There is no rule giving an optimal value, and empirical adjustment is the most reliable technique. It is also possible to learn λ from the data by cross-validation [START_REF] Craven | Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation[END_REF][START_REF] Hutchinson | Smoothing noisy data with spline functions[END_REF]. If we consider a stochastic framework where the data come from a random stationary process with additive random noise, both with known spectral characteristics, it is suggested in [START_REF] Eldar | Non-ideal sampling and interpolation from noisy observations in shift-invariant spaces[END_REF] to choose λ inversely proportional to the signal-to-noise ratio. This may serve as a heuristic in the deterministic case. This is confirmed in the example shown in Fig. 6, where λ = σ yields the minimum root-mean-square error, when approximating an unknown signal from its samples contaminated by additive noise with standard deviation σ.

Applications

There are plenty of problems where it is useful to fit a parametric curve on discrete data. Numerous methods have been proposed, generally issued from statistical estimation theory [START_REF] Wahba | Spline models for observational data[END_REF]. Our approach allows to reconstruct a function f T that is resolution dependent, a feature that offers many advantages in practical applications. For example, if a non-uniform signal is to be rendered on a display device with point-spread-function Γ(t), it is straightforward to apply the proposed approach: we choose ϕ = Γ and match T to the resolution of the device, so that V T (ϕ) is the set of all continuously-defined signals that can be rendered by the device. Therefore, by minimizing the criterion in (3), we ensure that our solution is optimal, given the available data.

Another potential field of applications is image analysis using multiscale "pyramidal" representations. The proposed work can be used for this task, if not only a single function but a whole collection of functions {f T }, with different values of T , is computed. For instance, the performances of procedures such as edge detection or image registration can be improved by processing from coarser to finer levels [START_REF] Thévenaz | A pyramid approach to subpixel registration based on intensity[END_REF]. To this purpose, our approach is more general than classical dyadic representations, that only apply to uniform signals and discrete dyadic resolutions (T = 2 n , n ∈ N).

Our approach is also particularly adapted to resampling problems involving rate conversions, such as image resizing. Let us present some generic problems that could benefit from our approach.

Non-uniform to uniform resampling

Given the non-uniform samples v[n] at locations x[n], suppose we want to obtain resampled values w T [k] located on the uniform reconstruction lattice (T k) k∈Z . We can simply compute the function f T (t) and resample it on the uniform lattice: this yields w T [k] = f T (T k) for every k. In fact, once the c T [k] have been computed, the signal w T is directly obtained by digital filtering [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF]:

w T = c T * b -1 , where b[k] = ϕ(k) for every k ∈ Z.
With this method, the representation capabilities of the target lattice, where the resampled signal lives, are exploited optimally: f T (t) retains at best the information contained in the non-uniform samples and representable on this lattice. Conversely, the information that is not representable is canceled out, and no irrelevant structure is introduced. As a consequence, the aliasing issue is automatically handled because, qualitatively, the frequencies higher than the Nyquist rate of the target lattice are not representable in V T (ϕ). That is why T is chosen so that f T has a resolution matched to the cut-off frequency of the reconstruction lattice. Note that we assume that ϕ is lowpass, and that λ is correctly chosen so as not to distort the frequency content of the data. This approach is also useful for resampling a uniform (x[n] = n) signal v. If v is to be resampled at coarser resolution or magnified with a non-integer factor, interpolation followed by resampling can introduce severe distortions [START_REF] Unser | Enlargement or reduction of digital images with minimum loss of information[END_REF]. Instead, we look for a function having the resolution of the target resampled signal, and not of the source signal. This formulation has been proposed previously in [START_REF] Aldroubi | Discrete spline filters for multiresolution and wavelets of l 2[END_REF], in the more restrictive case where T is an integer.

Interestingly, our approach amounts to formulate resampling as an inverse problem: we seek the uniform signal w T that, when interpolated by the function f T (t) and resampled back on the source lattice (x[n]), is the closest to the initial signal (v[n]). This is the opposite of the forward approach that fits a function on the source signal, typically by interpolation, and then resample it, independently of the target lattice.

Warping

Another application is signal or image warping. Suppose we have the uniform samples v[n] = s(n) of an unknown function s(t) at our disposal, and we want to compute the uniform samples (s(W -1 (n))) n∈Z of the warped function s(W -1 (t)) for some reversible transform W(t). The classical method used for warping, called backward mapping, consists in interpolating the uniform samples v[n] with a function f (t), that is then warped and resampled, providing the warped signal (f (W -1 (n))) n∈Z . Artifacts may appear with this method because the spectrum of the warped function f (W -1 t) extends potentially beyond the Nyquist frequency of the target lattice, resulting in aliasing when this function is sampled. In fact, the samples v[n] = s(W -1 (W(n))) of the warped function at the non-uniform locations x[n] = W(n) are known, so we are back to a non-uniform to uniform resampling problem, solved as previously. For example, arbitrary time-delay is a particular case of warping, with W(t) = t + τ , for a delay τ ∈ R.

The proposed algorithm can be used for image warping as well, as long as the 2-D warping W(x, y) is separable, by applying the resampling algorithm on each row and then on each column of the image. It is more rigorous to use a truly 2-D approach with a 2-D variational criterion, as proposed in [START_REF] Arigovindan | Variational image reconstruction from arbitrarily spaced samples: A fast multiresolution spline solution[END_REF][START_REF] Vázquez | Reconstruction of nonuniformly sampled images in spline spaces[END_REF], but one has to give up our direct and fast algorithm in that case. As an example, the image Barbara is warped with our approach in Fig. 7 with W(x/511, y/511)/511 = (2.2x-3.6x 2 + 2.4x 3 , 0.5y 2 + 0.5y), for (x, y) ∈ [0, 511] 2 . As can be seen in Fig. 7 (b), (e), aliasing appears on the trousers with backward mapping, because the stripes in the initial image are not representable any more in the warped image: their frequency has become higher than the Nyquist frequency of the target lattice. Our method is free from this drawback, and renders a flat zone without aliasing.

Conclusion

In this article, reconstruction from non-uniform samples is stated as a variational problem in a shift-invariant space. This formulation yields a well-defined solution, whatever the sample locations, which may be completely arbitrary. Moreover, the reconstructed function has a given resolution that can be tuned, for example to match the representation capabilities of a target lattice for resampling purpose, or the desired rate for coding. Taking into account the particular structure of the problem, we have proposed an efficient algorithm that computes the exact solution of the optimization problem, by solving a band-diagonal linear system, without having to represent and store any matrix explicitly.

8 Proof that A is positive definite First, M T M is positive semi-definite, i.e. u T M T Mu ≥ 0 for every vector u, since u T M T Mu = Mu 2 ℓ2 ≥ 0. Q is also positive semi-definite: let us choose a vector u and define g(t) = K+W -1 k=-W +1 u[k]ϕ( t T -k). Then u T Qu = I |g (r) (t)| 2 dt by construction of Q, and this integral is positive. Now, let us show that A is positive definite, i.e. u T Au = 0 ⇒ u = 0. We suppose that u T Au = 0. Then u T M T Mu = 0 and u T Qu = 0. This yields g(x[n]) = 0 for every n ∈ [0, N -1] and I |g (r) (t)| 2 dt = 0. Therefore, g (r) (t) = 0

2 ).

 2 We define b d as the discrete centered B-spline of degree d: b d [k] = β d (k) for every k ∈ Z and d ∈ N. We recall the expression of the first few discrete B-splines in the z-domain: B 0

  us define A def = M T M + λQ and y def = M T v. Minimizing the cost Ψ(c) amounts to solving the linear system Ac = y.

Figure 1 :

 1 Figure 1: A finite data set of N = 3 samples extended to an infinite one using mirror boundary conditions

Figure 2 :

 2 Figure 2: Computation time for the calculation of the coefficients c T [k], k ∈ [0, K], in log-log scale, for a reconstruction problem from N samples randomly located in the interval [0, 100], as in Fig. 6 below (T = 100/K, λ = 0.01). (a): K is fixed to 100 and N is varying. (b): N is fixed to 10000 and K is varying. The dotted line is for linear spline reconstruction (ϕ = β 1 , r = 1), while the solid line is for cubic spline reconstruction (ϕ = β 2 , r = 2).

Figure 3 :

 3 Figure 3: Uniform splines with knots at the integers (T = 1, λ = 0.01) fitted on 7 point samples in the interval [0, 10], with different polynomial degrees and values of the regularization parameter r. (a): ϕ = β 1 , r = 1. (b): ϕ = β 3 , r = 1. (c), (d): ϕ = β 3 , r = 2. Mirror boundary conditions are used for (a), (b), (c), not for (d).

Figure 4 :

 4 Figure 4: Uniform linear splines with different resolutions (ϕ = β 1 , r = 1, λ = 0.01) fitted on 7 point samples in the interval [0, 10]. (a): T = 0.1. (b): T = 1. (c): T = 2. (d): T = 5. The splines have their knots at the T k, k ∈ [0, 10/T ].

Figure 5 :

 5 Figure 5: Uniform cubic splines with knots at the integers (ϕ = β 3 , r = 2, T = 1, mirror boundary conditions) fitted on 7 point samples in the interval [0, 10], for different values of the smoothing parameter λ. (a): λ = 1.0. (b): λ = 0.001. (c): λ = 0.0001. (d): limit case when λ → 0.

Figure 6 :

 6 Figure 6: Reconstruction of the function s(t) = sin(( t 30 ) 3 ) (dotted line) from 100 noisy samples (additive Gaussian noise with standard deviation σ) at random locations in the interval [0,100]. (a): uniform cubic spline with knots at the integers (ϕ = β 3 , r = 2, T = 1) fitted on the noisy data (σ = 0.1, λ = 0.1). (b): RMSE e = ( 1 100 100 0 |s(t) -f T (t)| 2 dt) 1/2 in function of the noise level σ and the smoothing parameter λ, in log-log scale.

Figure 7 :

 7 Figure 7: Image Barbara (a) warped using: backward mapping (b), and proposed approach with λ = 10 -4 (c) (see text). (d), (e), (f): zoom-ins on (a), (b), (c) respectively.

  ,

	1.4												1.4
	1.2												1.2
	1												1
	0.8												0.8
	0.6												0.6
	0.4	0	1	2	3	4	5	6	7	8	9	10	0.4	0	1	2	3	4	5	6	7	8	9	10

within I, and g(t) is a polynomial of degree less than r in this interval. This polynomial has N roots at the x[n] (with at least r distinct roots by hypothesis on the samples locations). Then g = 0 and u[k] = 0 for every k.