Time Warping and Interpolation Operators for Piecewise Smooth Maps
A warping operator consists of an invertible axis deformation applied either in the signal domain or in the corresponding Fourier domain. Additionally, a warping transformation is usually required to preserve the signal energy, thus preserving orthogonality and being invertible by its adjoint. Initially, the design of such operators has been motivated by the idea of suitably generalizing the properties of orthogonal time-frequency decompositions such as wavelets and filter banks, hence the energy preservation property was essential. Recently, warping operators have been employed for frequency dispersion compensation in the Fourier domain or the identification of waveforms similarity in the time domain. For such applications, the energy preservation requirement can be given up, thus making warping a special case of interpolation. In this context, the purpose of this work is to provide analytical models and efficient computational algorithms for time warping with respect to piecewise smooth warping maps by transposing and extending a theoretical framework which has been previously introduced for frequency warping. Moreover, the same approach is generalized to the case of warping without energy preservation, thus obtaining a fast interpolation operator with analytically defined and fast inverse operator.
Time Warping, Interpolation, Perfect Reconstruction, Frames.
Signal processing spans a large variety of theoretical and application frameworks. Time-frequency tools address the necessity of representing signals in suitable domains in order to highlight specific features or properties. The design and introduction of these tools can be driven by novel abstract concepts and paradigms, emerging applications or the sake of computability and implementability. Warping techniques have been introduced as a new theoretical concept pursuing the idea of generating new Time-Frequency (TF) unitary representations, starting from already known ones such as wavelet analysis or filter banks [1, 2]. In this framework, rather than designing TF analysis tools specifically suited to a class of signals or aiming to extract specific TF features, an orthogonal deformation is applied to either the time or the frequency axis of the input signal followed by a standard TF analysis tool. This operation is equivalent to composing the warping operator with the TF operator. Since the composition of two subsequent orthogonal operators is still an orthogonal operator, a new orthogonal operator featuring the prescribed properties is obtained. This idea has been partially exploited but has not revealed any disruptive innovation [3, 4, 5]. Although this approach is conceptually appealing, many limitations arises from both application and computational point of view. In more detail, when dealing with Frequency Warping (FW), operations over the frequency axis implies the knowledge of the entire signal, hence working on finite length or windowed signals is mandatory. When dealing with finite spaces, warping operators can be designed only as tight frames rather than unitary operators , meaning that the obtained representation features some redundancy. Finally, with regard to computation, some issues arise when moving from the definition of warping in continuous domains to its implementation in discrete domains.
Despite the above considerations, warping has been effectively applied in various signal processing areas. For instance, current successful and popular implementations of warping attempt to compensate for inherent deformation occurring on either the time or the frequency axis because of physical phenomena such as dispersive propagation [7, 8, 9, 10, 11, 12, 13]. In those applications, theoretical and computational aspects become secondary as the knowledge of the physical parameters or functions, upon which the design of the warped operator is based, is often not accurate enough to require properties such as orthogonality and perfect reconstruction [14, 15]. Moreover, as the main goal is to invert an axis deformation, the energy preservation requirement does not always apply. Other applications include the Constant-Q Transform, where the invertibility is a major issue  and image compression .
Previous works by the authors about FW [18, 19, 20] mainly aimed at solving implementation issues. In fact, a thorough investigation of computational aspects and signal representation issues involved in the design of FW as a frame has been pursued. More specifically, with the aim of obtaining an approximation of FW being more accurate than the one obtained by simply sampling the frequency axis, an interesting theoretical framework about the representation of signals having non-smooth Fourier transforms has been identified together with a fast algorithm for FW computation. In addition, the same framework has been exploited for obtaining an analytical representation of the dual frame allowing perfect and fast reconstruction.
The work we here propose aims to move to Time Warping (TW) the mathematical framework which has been identified for FW, i.e. direct and inverse transform with respect to piecewise smooth warping maps. In this context, moving from the Fourier domain to the time domain arises theoretical and computational differences about the continuous representation of the warped axis and about design constraints which will be dealt with in the paper. As a natural completion of this framework, we also extend the same theory to the identify a time interpolation operator and its inverse.
The paper is organized as follows. Section 2 reviews basic concepts about warping operators and describes the problems being tackled in this work. In Section 3 the transition from FW to TW is detailed whereas in Section 4 warping is revisited as interpolation. Finally, some experimental results and conluding remarks are shown in Section 5 and Section 6 respectively.
2 Warping Review and Goals
First a review of the basics of warping are presented [21, 22], introducing a new perspective in order to focus on the difference which arises when moving from FW to TW. Moreover we also introduce the mathematical model for employing the warping framework to perform interpolation and its inverse.
2.1 Warping and interpolation in continuous spaces
The core transformation of warping is described by the following deformation operator applied to either the time or the frequency axis
where is a bijection, hence its derivative is always finite and positive, thus the inverse map also exists. Operator acts as warping if the integration is performed with respect to the second variable of the kernel
whereas it acts as unwarping if applied to the first variable, i.e. when the adjoint operator is considered (operator (1) is real so the transpose is equal to the adjoint)
hence the direct operator built from the inverse map is equal to the adjoint (or inverse) operator built from the direct map. By setting and equal to the operator corresponding to , we have . This property becomes more intuitive if one gives up to the orthogonalization factor and considers the neat interpolation problem. To show this, let us introduce the generalized operator
with representing a generic power , such that the simple interpolation can be referred as . The natural way for inverting would be to apply . Less intuitively, can be equivalently inverted by applying , i.e.
being the identity operator. In fact, the composition can be computed by
Obviously for these continuous operators it must hold . Instead, when coping with discrete-time signals, considering the approximated operator derived from rather than from brings to different results and accuracies. Section 4 addresses the problem of approaching the inverse of the interpolation operator.
2.2 Warping in periodic spaces
The continuous approach described in 2.1 has two major limitations. The warping operator (1) and its orthogonality rely on an infinite continuous axis. As a first step, one has to switch from an infinite axis to a limited interval in order to deal with the fact that input signals are known on a limited time interval for TW or are band-limited signals for FW. Relaxing the continuity is more delicate and will be dealt with in 2.2.
As far as FW is concerned, working with a limited axis is quite natural. In fact, a time-continuous band-limited signal can be safely sampled its spectrum will be continuous and periodic. In addition, it is required that the number of non-zero samples is finite, so that the signal is actually manageable. So, the FW operator designed for discrete-time signals always corresponds to a periodic-wise map. The periodic-wise map, with respect to a normalized period of length verifies the following property
For FW, the warping function must be also odd in order to transform real signals into real signals . The key concept introduced about FW is that , although the map is designed on a single period, its properties have to be evaluated on its periodic-wise extension. So, a map being smooth on a single period becomes only piecewise smooth when considered on the entire frequency axis.
With reference to TW, let us now consider a discrete-time signal in time domain. A continuous representation could be obtained by ideal interpolation, i.e. by applying an ideal rectangular filter to its spectrum. A limited time interval could be isolated by time domain windowing. This approach would be formally correct but the resulting time-continuous representation would depend on samples not belonging to the considered time-window. In order to be able to represent TW as an algebraic operation with respect to an input signal of finite dimension, we consider as continuous-time representation the circulant interpolation. By doing so, as for FW, TW can be described as a deformation of the entire time axis performed by means of a suitable periodic-wise map. The time axis deformation is represented in Fig. 1 with respect to a single time domain sinusoidal component. A time map being smooth with respect to a single period is guaranteed to be only (continuous with non-continuous derivatives) with respect to the entire time axis. Hence, the orthogonalizing factor in equation (1) and so the warped signal cannot be guaranteed to be unless the warping map features special regularities conditions at the period boundaries. Hence, TW maps are generally assumed to be piecewise smooth. As an example, in Fig. 2, the process of warping the same signal as in Fig. 1 by means of a globally smooth warping map is represented. The resulting spectrum has a fast decay, thus the signal is practically bandlimited and the formal description in the periodic space does not involve any advantage. As previously done for FW in [21, 22], the aim of this paper is to work with piecewise smooth maps like splines as they allow for a much flexible design.
2.3 Warping in finite-dimensional spaces
The operations described so far applies to finite-support domains but still involve continuous integration as in equations (2) and (3), which obviously is not computable with a finite procedure or representable with a finite series unless the continuous function is isomorphic with a finite space. The intuitive way to decrease the space dimension from infinite to finite is to perform a sufficiently dense sampling on . This approach will be referred to as Sample Without Filtering (SWF). If the warping map is only piecewise smooth over the whole axis, the function spectrum has a polynomial decay because of the warping function singularities, as represented in Fig. 1. Hence, from a pure theoretical point of view is not band-limited and cannot be sampled without aliasing. The less intuitive but more accurate way to perform the space dimension reduction is to first apply a filter over followed by a suitable dense sampling. This strategy will be referred to as Sample After Filtering (SAF). The SWF can be simply implemented by means of Non Uniform Fast Fourier Transform (NUFFT) algorithms [23, 24, 25, 26], whereas the SAF requires some further processing .
For FW, the function is complex and its spectrum is actually a time domain interval. Conversely, for TW is in time domain, hence its spectrum is actually a proper frequency spectrum. In 3.3 we will detail afterwards how to select the filtering band for TW and FW.
According to what we have described so far, TW can be modelled in the following way. First we introduce supporting Fourier and windowing/filtering operators. As a convention, all operators applied to discrete-time domains will be represented with vector notation, whereas warping operators applied to continuous-time domains have been represented by gothic letters. We first introduce the Fourier series
and the Discrete Fourier Transform operator
where by now is a suitably defined set of contiguous integer (details will be given in 3.3). We also introduce the warped Fourier series
and the time windowing/filtering operator
Finally, the approximate TW operator can be computed by performing the warped circulant interpolation , evaluating its Fourier series by F, selecting spectrum component, , and going back to the time domain by , that is
where the subscript stays for time.
As far as FW is concerned, the problem is slightly simpler as there is no need to perform a circulant interpolation and to go back to the warped domain. We have
where represents a warped Fourier transform for a discrete-time signal and the inverse Fourier transform.
2.4 Paper goals and results
For FW, a way to approach SAF by compensating aliasing on SWF has been modelled and analysed in . A decomposition allowing for a fast computation has been also provided. In , the dual operator such that has been identified by a analytical model. The main goal of this work is to transpose the results obtained for FW into TW together with some generalizations and expansions. In more detail, in Section 3 we will focus on (i) how to apply to the model and algorithm obtained for , (ii) how to transpose the algorithm for into and (iii) how the choice of input and output domains differently impact on TW and FW.With respect to (iii), we will also extend FW to input and output domains which have not been covered in previous works. Furthermore, in Section 4, according to the notation used in equation (4), we will also detail how to generalize operator , which could be also referred to as , into the generic operator , with special reference to the pure interpolating operator and its inverse .
3 From Frequency Warping to Time Warping
Rather than first providing a review of FW and then redefining the model for TW, we provide results obtained for FW together with their modifications for TW and generalizations.
3.1 Model transition for the direct operator
The core idea of the method which have been proposed in  is to obtain the operator corresponding to SAF approach by correcting the operator obtained by SWF. To describe both TW and FW with the same approach, we introduce the following shift-variant operator
which features a center of symmetry for and as it can be easily verified that . Operator (13) is a unitary operator, i.e. it satisfies the requirements for getting the exact inversion by means of its adjoint as it represents a pure transposition of (1) into discrete-time domains according to what has been described in 2.2. According to (12), can be obtained by windowing both its input and its output. By noticing , we can write
and is usually taken in the interval in order to include the center of symmetry. In previous work on FW only the symmetrical case with even has been considered. The extension to a generic is dealt with in 3.3. In order to relate W to , we highlight that , then equation (11) can be rewritten as
Hence, we point out the following principle. Every result obtained for FW involving W can be employed for TW by simply applying a conjugation to W and including it between discrete Fourier transform and its inverse of size and respectively.
The SWF approximation for FW is obtained by sampling the integral in (13) in points and will be referred to as
and can be described by
where is the Warped Discrete Fourier Transform
where and have been neglected as both input and output of are already band-limited. The above result can be considered obvious, as it represents the obvious way to perform a warped circulant interpolation. Operator and , as well as and , are computed in a similar way but have quite different characteristic as and behave as circulant warped interpolators whereas and appear as dispersive delaying operators.
Following the equivalence which has been highlighted in (3), the sampled TW operator corresponding to the inverse map is also introduced
As it has been anticipated in Section 2, although for the operators in continuous spaces corresponding to and it holds , for the sampled operators we have . Nevertheless, can be still used as an approximation of the inverse operator, i.e. . An accuracy comparison will be provided later in this Section.
The computational strategy for consists in first finding a model for the decaying tails of W for , then using this model to compute aliasing to correct operator . The error operator which is introduced when switching from operator to operator is represented by
and aliasing is obtained by periodic sum over E
such that . This decomposition is schematically represented in Fig. 3, where also the typical sparsity pattern and decay of matrix W is highlighted. We refer to the singularities of the considered piecewise smooth map as . Under certain conditions which will be detailed in 3.3, E and A can be factorized as follows
where and are diagonal matrices having as main diagonal and respectively, is a lower triangular kernel matrix depending on the various order derivatives of the warping map on (see equations (41)-(43)) while V, Y and U are real matrixes obtained by sampling fixed functions not depending on the warping map. V, Y and U will be reviewed in 3.3 according to the general input and output indexing (14) introduced in this work (see equations (37), (38) and (40)).
Hence, by following the same principle as in (15) we have
By taking advantage of (23), the generic time aliasing component referred to can be rewritten as
so that U and V, being constant matrixes which may be precomputed, can be replaced with their Fourier transforms performed along columns and rows respectively, whereas and , being modulations, are turned into circular shifts.
3.2 Model transition for the dual operator
Operator is designed in such a way that provides a very good approximation of the dual operator. Nevertheless, the exact dual operator of has been found in  by exploiting the Neumann series
and recalling that . By further calculations one gets
where H is a block column matrix whose items are and Z is deterministically obtained by H, Y and . In a similar fashion, we have
which can be obtained again by applying the principle highlighted in (15)
The generic -th block item of square block matrix results
where the same matrix used for aliasing in (25) has been highlighted, meaning that can be computed by using the same basis matrixes employed for aliasing compensation.
In order to exemplify the accuracy which can be achieved according to the various considered inversion strategies, reconstruction errors with respect to a random input signal have been plotted in Fig. 4. A complete analysis of the error behaviour should take into account parameters such as the redundancy (see [27, 22]). Nevertheless, Fig. 4 shows that is about as accurate as , but suggests that the corresponding error matrix may not have a low rank. A thorough error analysis will be provided in Section 5.
Finally, as far as the computation is concerned, the whole FW to TW transition can be summarized by the following operations:
replace basis vectors U and V with their Fourier transforms
replace kernels S and Z with their conjugates
replace modulation matrixes P and Q with corresponding circulant shifts.
Other issues regarding convergence and constraints on windowing/filtering will be discussed in next subsections.
3.3 Input and output domains indexing
The aim of this subsection is to focus on how differently the sets and have to be chosen in FW and TW with respect to the decomposition and algorithm involved in the SAF approach. When introducing (11) and (12) and more remarkably (22) and (23), constraints for invertibility and representation convergence have not been given. So, the feasibility conditions will be reviewed and updated to the generalized setting (14).
A necessary condition for invertibility is expressed by
This condition can be empirically obtained by imposing that the sampling (16) of the continuous axis in (13) is capable of catching the band enlargement brought by warping. However, by analysing the structure of matrix W (see Fig. 3), it turns out that condition (29) is not sufficient, as a similar condition has to be satisfied separately for positive and negative indexes
As far as FW is concerned, it is worth highlighting how the indexing of the input domain affects the way the temporal component of the input signal are treated. From a qualitative point of view, samples corresponding to positive indexes are processed in a causal way, whereas samples corresponding to negative samples are processed in a non-causal way. So, if , FW behaves qualitatively as a causal transformation. Conversely, from a quantitative point of view, every warping map belonging to the class of piecewise smooth functions acts on positive indexes in a non-causal way because of the slow decay of along (this is true also for smooth maps except the map , where is a suitable parameter , being the only purely causal map). So, the non-causal behaviour is inherent with the considered warping maps and implies that the output domain has to be chosen accordingly. However, it is interesting to analyse if the proposed computational approach could be used with input and output domains as close as possible to and .
As far as TW is concerned, the input and output domains have to be chosen to guarantee that equation (15) produces a real signal. The symmetrical choice and intuitively suits this requirements as and operate in the frequency domain, meaning that the indexes in and are frequencies and for each positive frequency a negative one has to be present. This and other constraints for TW will be detailed later in this Section.
The set represents a sampling over a period of the time domain for FW and of the frequency domain for TW. In digital signal processing a major difference holds between sampling an interval in even or odd number of points. Therefore, it is sensible to introduce the period boundaries not dependant on the particular sampling choice. To do this, the interval covered by is thought at as a continuous interval of length , whose left and right boundaries are
As instance, when is even, the covers the interval , whereas when is odd covers the interval . Then, the relative shift representing how much the set differs from the set is introduced
When is equal to , is equal to for both even and odd . So, represents the symmetric indexing case, while represents the causal case. We also highlight the following
So, conditions (30) can can be rewritten in the following way
These conditions are necessary for both SWF and SAF. However, for the SAF approach to be implementable by taking advantage of decompositions (22) and (23), conditions (36) are necessary but not sufficient. To show this and obtain a sufficient condition we review and update the definition of V, Y, U and S. Following, the approach pursued in previous works, V, Y and U are defined as sampling of continuous non-divergent functions. In more detail, the -th column of matrix V is obtained by sampling the polynomial , , defined on , hence it must be normalized with respect to its maximum which occurs on one of the boundaries and is identified by (34), hence we set
The -th row of matrix Y is obtained by sampling the function , , defined on , whose maximum occurs in (35), so it follows
Finally, The -th row of matrix U is obtained by sampling the derivatives of the following function
In  the kernel referred to a generic singularity has been also decomposed in two terms, one depending on the map only and the other one depending on the parameter setting, i.e. the relative redundancy . The first term remain unchanged with respect to the generalized indexing (14), whereas the second one has to be redefined in order to include the relative shifts and . So, is given by
where represents the element by element product, is a lower triangular matrix evaluated either in or in which will be detailed later
ans is a scaling matrix whose elements are given by
whose decaying behavior is mainly affected by
When , is guaranteed by (29). For different choices, it might be required to increase . It is worth pointing out that condition (45) refers to the derivative in singularity points, i.e. and not to . If occurs in one of the singularities, then (45) implies both (36).
The parameter , being responsible for the vanishing of , also determines how to truncate the dimension of in order to obtain a numerically accurate representation. So, for a kernel vanishing in the same way as for can be obtained only with a larger relative redundancy . As instance, with , the number of output samples needs to be times larger. In a real application, one can assume that and are set by specifications, while and can be adjusted in order to meet a target , which implies constant . From eq. (35) and (31)-(32) one gets . So, any choice of and leaves unchanged. With respect to the pure causal case and , one would have , thus but still and equal to the case and . These considerations are exemplified in Fig. 5, where different choices of the output domain are considered.
This short analysis clarifies the fact that for FW there is a certain degree of freedom in the choice of the input indexing, whereas the most convenient output indexing is always the same as if the input indexing was enlarged and made symmetrical. This makes the employment of the generalized indexing quite unpractical unless computational complexity is not an issue and the focus is just on getting a correct numerical representation of the transformation.
We now go back to TW and the way it is conceived starting from FW. By comparing (11) and (12), the choice of the indexing regards the output of , as it is fed into operator , which is not shift-invariant. Hence, and must be chosen in such a way that if the input of is a spectrum of a real signal, then its output has to be a spectrum of a real signal too. The property is easily verified for the core operator W. In fact, by introducing operator R performing the indexing reversal
one has to impose that, with respect to an input signal s such that , the output is also invariant with respect to the application of reversing and conjugation, that is
which is verified if
The above condition is the same as stating that , which has already pointed out. Conversely, for the truncated operator , the condition is rewritten as
where and are now square matrixes of size and . For this condition to be verified the following choices are required
Symmetry on the time domain comes trivially from the Fourier transform. Less obviously, only odd and are allowed. Critically sampled signals, i.e. even, could be still dealt with by resampling them on frequency points by splitting the -th frequency coefficient over and . In a similar fashion, an even could be obtained by forcing an additional oversampling after warping.
Finally, we can summarize the results about domain constraints by the following statement. The employment of SWF for FW comes with the cost of employing a non strictly causal transformation. This drawback can still be overtaken in some time-frequency analysis application such as Constant-Q Transform (CQT), while it represents a major limitation for the applications where FW represents a way to compensate for physical phenomena as dispersive propagation. Conversely, the employment of SWF for TW comes with the minor constraints of taking only odd values for and .
4 Warping as Interpolation
As discussed in 1 and 2.1, the entire design of warping operator is driven by orthogonality, but, in certain cases, the presence of the orthogonalizing factor is not practical either for TW and FW, whereas the availability of the inverse operator is still necessary. The suppression of the orthogonalizing factor, i.e. considering the interpolation operator , does not allow any more to employ the transpose operator for recovering the original signal and requires instead an approximation of the continuous operator , as shown in 2.1. The discretization of brings to , whereas the discretization of brings to . To clarify the difference between these approaches, operators , and have been depicted in Fig. 6 together with their sparsity pattern in the frequency domain. In the time domain, both and , being pure interpolators, appear as a constant warped diagonal matrix. Operator differs from for an amplitude factor, hence in the frequency domain they share the same sparsity pattern. The sparsity pattern of originates from aliasing along rows, thus limiting the reconstruction performances as it will be shown in Section 5. For this reason, we focus on extending the results obtained for to .
With reference to SWF for TW, the difference between and is merely an amplitude factor. Conversely, when considering SAF, the difference between and involves reviewing the decompositions (22) and (23). From a practical point of view this can be accomplished by redefining the matrix (42) employed to express the kernel (41). Rather than reviewing the whole strategy for obtaining the decompositions, we just recall that the factorization model comes from the possibility of expliciting the dependency of on by
where stays for the differential value and , that is
The model (46) is responsible for matrix Y (38). By expressing the derivatives in a symbolic way with respect to the derivative order, matrix U (37) can be explicated and the value of coefficients in (42) found. As already done in (4), we here replace the power with a generic power and show that the resulting expression is compliant to the one obtained for . In 2.1 we already pointed out that the value of being interesting from a practical point of view are and . For these two values, the adopted model would not work for reasons which will explained in 4.2. The workaround employed here consists in solving the problem symbolically also with respect to . This is also the reason why we have considered the soft orthogonalizing factor in (4). Moreover, we here detail the symbolic algorithm for finding the coefficients .
Before dealing with the decomposition (46), we aim to show that, given the possibility of expressing and with a SAF approach, then the analytical dual operators can also be defined.
4.1 The generalized dual operator
As a starting point, the basic operator W in (13) is redefined by changing the exponent of the orthogonalization factor from to a generic exponent
All related operators such as A and E will be referred to in the same way. As shown in (5), inversion can be still obtained by applying operator where
By restricting the input by and employing the splitting between and E as done in (20) we get
which can be used as in (26). It is worth noting that
because truncations are operated separately, nevertheless is still an approximation of the identity as accurate as , hence its inverse can still be calculated by means of the Neumann series
Most importantly, in 4.2 and 4.3 it will be shown that and share the same decomposition as reported in (22) but with different kernels. As a consequence, equation (27) can be rewritten by just suitably replacing H with either or and kernel Z with , which is properly recomputed by using the pair :
Finally, the extension to TW is simply obtained by applying the same procedure as in 3.2, i.e.
The computation of can be easily inferred by .
4.2 Polynomials identification
The problem of expressing the derivative of the composition of two functions, being here the exponential function and the warping function with some additional coefficients and factors as represented in (46) and (47), has a notable solution known as Faà di Bruno’s formula, which also has a simpler form in case the first function is the exponential. Nevertheless, to the best of our knowledge the decomposition proposed here is remarkably different since the inherent recursive structure of the exponential function derivatives is exploited. Moreover, the proposed approach is more effective from the algorithmic point of view with respect to our purposes. In fact, we specifically aim to express the derivative of (47) as a power series of , which allows for the decomposition (22) and for an effective characterization of convergence . We also point out that this symbolic computation cannot be performed by any of the symbolic softwares available.