Local cubic splines on non-uniform grids and real-time computation of wavelet transform

In this paper, local cubic quasi-interpolating splines on non-uniform grids are described. The splines are designed by fast computational algorithms that utilize the relation between splines and cubic interpolation polynomials. These splines provide an efficient tool for real-time signal processing. As an input, the splines use either clean or noised arbitrarily-spaced samples. Formulas for the spline’s extrapolation beyond the sampling interval are established. Sharp estimations of the approximation errors are presented. The capability to adapt the grid to the structure of an object and to have minimal requirements to the operating memory are of great advantages for offline processing of signals and multidimensional data arrays. The designed splines serve as a source for generating real-time wavelet transforms to apply to signals in scenarios where the signal’s samples subsequently arrive one after the other at random times. The wavelet transforms are executed by six-tap weighted moving averages of the signal’s samples without delay. On arrival of new samples, only a couple of adjacent transform coefficients are updated in a way that no boundary effects arise.


Introduction
Since their introduction in [14], splines have become one of the most powerful tools in mathematics and computer aided geometric design. In recent decades, splines have served as a source for the constructions of wavelets, wavelet packets and wavelet frames. Splines and spline-based wavelets, wavelet packets and frames are extensively used in signal and image processing applications (see, for example, [3,4]).
Interpolating splines possess exclusive approximation properties. In particular, an interpolating spline of order p, which consists of pieces of polynomials of degree p − 1, restores polynomials of the same degree. Due to this property, these splines generate biorthogonal wavelets with p and, in some cases, p + 1 vanishing moments ( [3,4]). However a drawback in the design and manipulation of interpolating splines is that, for their computation, a system of equations that involves all the available grid samples has to be solved. This fact prevents the usage of these splines for real-time processing. 1 Therefore, the idea to have splines that can be designed and manipulated directly without resorting to systems of equations, while their approximation accuracy is close to that of the interpolating splines, is attractive. A method for the design of such splines on a uniform grid was presented in the pioneering spline paper [14]. When the grid is non-uniform, the design and estimation of the approximation properties of such splines is more complicated, especially for higher-order splines. A number of investigations in this field were carried out in the 70's (see, for example, [13]). These splines are called local because the computation of a spline's value at a fixed point requires to use only a few adjacent grid samples. Nevertheless, there exist local splines which provide the same approximation order as the interpolating splines. That is, there exist local splines of any order p which restore polynomials of degree p − 1. Such local splines are referred to as quasi-interpolating splines.
In this paper, we describe a procedure to design and analyze local cubic quasiinterpolating splines on arbitrary grids. Typically, local splines are designed via their B-spline representation. An alternative approach presented in this paper, which is based on the relation between quasi-interpolating splines and cubic interpolating polynomials, results in a "simple" algorithm for spline computation. When samples are given on a limited interval, which is a typical situation, the spline is extended to the boundaries of the interval without loss of its approximation accuracy. In addition, a method for accurate extrapolation of the spline beyond the sampling interval is presented. Only six adjacent samples are needed to compute the spline's value at a certain point. Therefore, the design of the spline can be implemented in a real-time mode when samples of a signal arrive dynamically and sequentially at random times. Due to the extension formulas, the design can be carried out without delay up to the latest sample arrival. In addition, the extrapolation algorithm can be applied to a prediction-correction processing of signals evolving in time.
Similarly to cubic interpolating splines, the described splines restore cubic polynomials at inner parts of the sampling interval and so also near the boundaries and in the extrapolation process. Moreover, the spline representation via interpolating polynomials provides sharp estimations for the approximation errors.
The described splines are used for the construction of wavelet transforms of arbitrarily sampled signals. Most of the existing wavelet transforms operate on uniformly sampled signals. A few works that describe wavelets on non-uniform grids appeared in recent years. Equidistant wavelets, which are utilized for denoising of non-uniformly sampled signals, appear in [6]. Scaling functions on different levels of enclosed irregular grids, which appear in [7], are designed as limit functions of subdivision and wavelets are designed as their linear combinations. In [11], a local decomposition of the space of splines, which are constructed on a fine irregular grid, onto a coarse-grid spline space and its orthogonal complement was presented. The bases of these subspaces are formed by quasi-interpolating splines and compactly supported pre-wavelets. A general method for constructing wavelets on irregular lattices in R d and a general atomic frame decomposition of the space L 2 (R d ) are described in [2].
In the current paper, the quasi-interpolating splines are utilized as a source of finite impulse response "filters" to generate wavelet transforms of discrete-time arbitrarily sampled signals. When sampling is uniform, these are filters in a proper sense of the word. Neither scaling functions nor continuous wavelets are involved.
A natural way to design and implement wavelet transforms of signals is the Lifting Scheme introduced in [18], which consists of subsequent applications of prediction and updating operators to signal's samples. The scheme was extended in [19] to non-uniformly sampled signals. The lifting wavelet transform on a non-uniform grid, where wavelets have two vanishing moments, is applied in [20] to signal denoising. Our design scheme is based on quasi-interpolating local splines. The idea which, in the case of uniform sampling, is explored in [4], is to predict odd samples of the signal to be transformed by values of a spline constructed on the even samples of the signal. Then, the detail coefficients are derived by subtraction of the predicted odd samples from the original ones. The next step consists of updating the even samples by values of the spline designed on the detail coefficients. Since quasi-interpolating splines restore cubic polynomials, the detail wavelet transform coefficients are zero if the signal is a sampled cubic polynomial (at least locally). In that sense, we claim that the arising discrete-time wavelet transforms have four discrete vanishing moments. Typically, wavelet transforms of signals (images) defined on limited intervals (areas) require an extension of the object beyond its boundaries ( [5]), in order to reduce boundary effects. However, by using the definition of the splines near the boundaries of the sampling interval, the wavelet transforms are implemented without this type of extension. The design of splines make them useful for real-time (with no delay) execution of wavelet transforms in the situation when samples of a signal subsequently arrive at random time moments.
A real-time denoising method, which is based on the lifting wavelet transform of uniformly sampled signals, is presented in [12]. The online wavelet transform is implemented using a moving window (see also [21]), which produces some delay with respect to the sample acquisition and requires an artificial extension of the data beyond the moving window. Our algorithm produces wavelet coefficients with no delay. No extension of the data is needed. Arrival of a new sample leads to the production of a new coefficient and to updating of a couple of already produced coefficients.
The paper is organized as follows. Section 2 introduces the necessary notation and recalls definitions of the divided differences and interpolating polynomials. Section 3 introduces local cubic quasi-interpolating splines, describes algorithms for spline computation on unlimited and finite intervals and an extrapolation method. Approximation properties of the splines are investigated and numerical examples are provided. The Lifting Scheme of the wavelet transforms is presented in Section 4.1. In Section 4.2, spline-based prediction and updating operators are explicitly described. A scheme for the real-time wavelet transforms computation is outlined in Section 4.3 and numerical examples are provided.
In order for the paper to be self-contained, some results about local quasiinterpolating splines, which appeared in [17,23], are included.

Preliminaries
The following notations are used throughout the paper: The grid on the real axis is denoted by t := {t[k]}. The steps of the grid are denoted by h[k] The n-order divided difference (DD) ( [1,9,16]) of a sequence f = {f [k]} with respect to the grid t is where If a function f (t) belongs to C n [t[k], t[k + n]], then the DD of the sequence

Local cubic splines
The cubic B-spline is It is supported on the interval (t[k], t[k + 4]). The following observation will be used for further design.
where q = {q[ν]} is a sequence of real numbers, which determines the spline's properties.
Typically, for a spline that approximates a function f (t), the coefficients q[ν] are derived from the samples {f [k] := f (t [k])} of f (t) on the grid t. In that case, the spline is denoted by s[f ]. For interpolating splines, the coefficients q[ν] are obtained by solving a tridiagonal system of equations, which involves all the available samples and some boundary conditions. An alternative is provided by the so-called local splines.  It is well-known that cubic interpolating splines restore cubic polynomials, thus having the approximation order 4. Definition 3.4 A local cubic spline, whose approximation order is 4, is referred to as the quasi-interpolating local spline.
To get the simplest cubic spline on the grid t, which approximates the function Denote by where h[k] is the grid step.

Proposition 3.5 ([22])
If the coefficients in Eq. 6 are chosen as then the spline

restores cubic polynomials. For the computation of the spline's values
Proof Proved by direct computation of the spline's values of the functions f (t) ≡ 1, f (t) := t r , r = 1, 2, 3.

Remark 3.6
In the rest of the paper, we deal exclusively with the quasi-interpolating splines s 1 [f ]. Therefore, we drop the index · 1 . Thus, in the sequel s[f ] := s 1 [f ].

Computation of quasi-interpolating cubic splines
The quasi-interpolating cubic splines s[f ] can be represented in an alternative form, which is based on a relation between the splines s[f ] and cubic interpolation polynomials. Such form is computationally efficient because it utilizes the wellestablished algorithms for the computation of interpolation polynomials, and facilitates extension of the spline to the boundaries of the intervals. In addition, that representation makes it possible to derive sharp estimates of the approximat errors for the splines s[f ]. Originally, this representation was introduced in [23]. For the paper to be self-contained, we place the proof of the following statement into Appendix.

Theorem 3.7 For t = t[k] + h[k] τ, τ ∈ [0, 1], the cubic spline s[f ](t), which restores cubic polynomials, is expressed by
where the coefficients Proof In Appendix. When the grid t is uniform and t[k] = hk, k ∈ Z, the finite differences are used instead of the DDs to express the spline s[f ](t). In this case, when t = h (k+τ ), τ ∈ [0, 1], the spline is represented explicitly by For the computation of the spline's values Thus, the spline can be computed in real time with a delay of two samples.

Approximation properties of quasi-interpolating cubic splines
Approximation properties of the spline s[f ](t) are close to the properties of the cubic spline s i [f ](t) which interpolates the function f (t) on the grid t such that If the function f (t) belongs to If the grid t is uniform, at least locally (that means that f (4)
Proof In Appendix.

Remark 3.9
The estimates in Eq. 14 is sharp in the sense that it becomes an identity for the function f (t) = t 4 when t = h(k + 1/2).

For comparison, a similar sharp estimate for interpolating splines
whereh is the maximal grid step at the interval I . That means that the estimates on a certain interval depend only on the function's behavior in a close vicinity of this interval.

Remark 3.11
The sharp estimate in Eq. 14 for the uniform grid was established in [23] by the technique different from that used in the proof of Theorem 3.8. The inequality in Eq. 12 makes it possible to establish sharp estimates of the approximation errors for functions from the classes C m [I [k]], m < 4. In the case when the available samples f [k] are noisy, the estimation in Eq. 12 enables us to evaluate the noise contribution into the approximation error.

Quasi-interpolating cubic splines at finite intervals
Assume that only a finite number of samples .
The constant ( Proof In Appendix.   3 . Then, the difference becomes

Remarks on the real-time spline computation
Due to the fact that no more than 6 adjacent grid samples are needed for the computation of the spline s[f ] at a point t, the spline can be computed in real time. It means that the spline follows the samples that arrive one after another at random times.  ]. This fact substantiates our claim that the spline follows the arriving samples.
The above scheme is illustrated in Example 2 in Section 3.6.

Examples
Example 1: Restoration of the sine function from randomly located samples: Figure 1 illustrates results from the restoration and extrapolation experiments of the function sin(t) + sin 2t from 30 randomly spaced samples (top plot) and from uniformly spaced samples (bottom plot). In both cases, the first and the last samples

Spline-based wavelet transform
The Lifting Scheme, introduced in [18,19], is a method that constructs bi-orthogonal wavelet transforms and provides their efficient implementation. The main feature of the lifting scheme is that all the constructions are derived directly in the spatial domain and therefore can be custom designed to more general, irregular settings such as non-uniformly spaced data samples and bounded intervals. In addition, the lifting scheme fits well to the real-time execution of the wavelet transforms. In this section, we outline the lifting scheme and describe how to use the local quasi-interpolating cubic splines, designed in Section 3, for the construction of wavelet transforms of non-equally sampled discrete-time signals and in the situation when signal's samples arrive subsequently at random time moments.

Discrete lifting wavelet transform
The lifting wavelet transform of a discrete-time signal can be implemented either in a primal or in a dual mode. We outline only the primal mode.

Primal decomposition
The lifting wavelet decomposition of a signal from the space l 1 consists of four steps:

Primal reconstruction
Reconstruction of the signal f from the arrays y 0 and y 1 is implemented in a reverse order: 1. Undo Normalization: a = y 0 / √ 2, d = √ 2y 1 .

Undo Lifting:
The even array is restored by e = a − U d.

Undo Predict:
The odd array component is restored by o = d + P e. 4. Undo Split: Restoration of the signal from its even and odd arrays by f = Merge{e, o}.
The lifting transform is perfectly invertible with any choice of the operators P and U . The direct and inverse transforms can be symbolically represented in a matrix form: where I is the identical operator and the analysis "polyphase matrix"M of the first decomposition level isM The inverse transform is represented by where the synthesis "polyphase matrix" M of the first decomposition level is It is readily seen that, once the operators P and U commute with each other, we have MM = I. Therefore, subsequent applications of the operatorsM and M to the vector · · · y 1 } is implemented in a reverse order.

Spline-based prediction and updating operators
Wavelet transforms generated by the lifting scheme are determined by the choice of the prediction P and the updating U operators. Splines provide flexible tools for the design of such operators. The idea is to construct a spline on the even sub-grid {t[2k]}, which either interpolates or quasi-interpolates the samples of the signal e. The values of this spline at the odd grid points {t[2k + 1]} are used for the prediction of the odd samples o. The next step is to construct a spline on the odd grid points, which (quasi-)interpolates the samples of the prediction-error signal d. Then, the values of this splines at the even grid points are computed. These values are used for updating the even samples e.
In the case of equally spaced samples, calculations are reduced to low-pass filtering of the corresponding arrays [4]. Application of wavelet transforms to signals sampled on finite grids and, in particular, to images, requires to extend the signals beyond their boundaries, otherwise distortion appears near the boundaries [5]. Various extension schemes have been developed to deal with the boundary effects of finite-length signals: zero padding, periodic extension and symmetric extension are basic extension methods. However, quasi-interpolating splines on finite intervals designed in Section 3.3, make it possible to implement wavelet transforms of signals sampled on bounded intervals without extending the signals beyond their boundaries.
The prediction and updating operations are reduced to the design of the splines s[f ] on different grids and computation of their values at intermediate points.

Prediction and updating operators on unlimited grids
To highlight the structure of the operators, we use a standard representation of the quasi-interpolating spline via the B-splines given in Eq. 8: where the coefficients β i [k] are defined by Eq. 7 and the cubic spline , then the sum in Eq. 30 comprises only 6 terms:

Prediction and updating operators on limited grids
Assume that the function f is sampled on a bounded interval, thus Assume that N is an odd integer (Fig. 4): A similar design is performed when N is an even integer. To implement the next step of the wavelet transform, the above operations are applied to the smooth array y 0 . As a result, the arrays y 0 [2] and y 1 [2] are produced. This process repeats itself.
The spline-based wavelet transform possesses the local "discrete vanishing moments" property, which is formulated next.  Figure 5 displays a randomly sampled signal and coefficients from its four-level wavelet transform. Note that the inverse wavelet transform from these coefficients provides a restoration of the signal where its maximal deviation from the original, which is displayed in Fig. 6, is 1.5 10 −12 .

Real-time execution of the wavelet transform
The transform scheme described in Section 4.2.2 makes it possible to execute the wavelet transform of a signal in real time. It means that the samples f [ν] of a signal  At the same time, the above procedures are applied to the smooth coefficient array y 0 [ν] , ν = 0, ..., (N − 1)/2, to obtain the arrays y 0 [2] [ν] and y 1 [2] [ν] , and so on. The transforms were applied to two signals x 1 and x 2 such that :  The 10-level lifting wavelet transform was applied to the signals x 1 and x 2 , then signals were restored from the transform coefficients by the inverse wavelet transform. The transforms were executed by non-optimized Matlab codes on the PC with the Intel(R) Core(TM) i7-3770 CPU@3.40 GHz processor. Herewith, the CPU time was 3.357637 seconds for the direct transform and 3.374883 seconds for the inverse transform.

CPU example
The differences between the original x 1 and x 2 signals and the restored X 1 and X 2 ones are displayed in Fig. 8. In both cases the differences do not exceed 1.5 · 10 −15 . We emphasize that no boundary effects arrived in process of the signals' restoration

Conclusion
The local cubic splines, which were designed in Section 3 and supplied with a simple fast computational algorithms, can serve as an efficient tool for the real-time signal processing. As an input, the splines use either clean or noisy arbitrarily-spaced samples. Sharp estimates of the approximation error were established. An important application of the designed splines is the real-time wavelet analysis of non-uniformly sampled signals.

Summary of the paper's contributions
• Simple algorithm for the real-time computation of the local cubic quasiinterpolating splines including a smooth extension of the splines to the boundaries of the sampling interval. • Formulas for the extrapolation of the splines beyond the sampling interval, while retaining the approximation accuracy. The formulas can be used for the prediction-correction signal processing. • Sharp estimates of the approximation errors for the local cubic quasiinterpolating splines. • Design of wavelet transforms of the arbitrarily sampled signals based on the local cubic quasi-interpolating splines. The wavelet transforms have four local "discrete vanishing moments" in the sense of Proposition 4.2.
• Real-time scheme of the execution of the wavelet transforms of signals whose samples arrive dynamically and sequentially at random times.
In future work, the approach will be extended to the higher-order splines. The designed splines and spline-based wavelet transforms will be applied to online denoising and compression of practical signals such as, for example, ECG signals. Another potential field of application of the developed techniques is the real-time process control including detection of specific transient events.