

BSpline based Continous Wavelet Transform 
Introduction  Description  Example 
1 Introduction
The Continuous Wavelet Transform (CWT) of a signal in the basis is defined as
It can be interpreted as the correlation of the input signal with a mirrored version of , which is rescaled by . For a 1D input, the result is a 2D timefrequency description of the signal. The scale is proportionally inverse to the frequency studied, and represents the time localization at which we analyze the signal. The bigger the scale the wider the analyzing function , and hence the smaller the corresponding frequency. The big advantage over FT is that the frequency description is time localized, and hence evolves in time. Fast implementations of the Wavelet Transform take into account multiresolution properties of wavelet functions to compute the CWT at dyadic scales and timeshifts . Using Mallat's algorithm or derivatives, they get an overall complexity. Other techniques compute the Wavelet Transform at dyadic scales and at linear timepoints, or, by using mother wavelets closely dilated, at Madic scales. Other methods compute the CWT at integer scales ( ).
We present a new method, using expansion of both the input signal and the mother
wavelet in BSpline basis, that computes the CWT at any real scale, allowing
to zoom in the frequency domain. This method also gets rid of recursivity scheme
present in the previous methods, and fits well to parallel implementation.
2 Review of Bspline properties
Basic Definitions
Continuous Operators :
Using the exchange rules for the onesided power function and for the shift with the rescaling operator given in , and noting that the rescaled finite differences operator can be expressed as
Finally, we define the mixted convolution of a continuous signal with a discrete one with intersample distance equals to as:
Bspline Interpolation
Splines are one of the most famous family of interpolators : given a sampled
signal, know for example at integer points, one easily interpolates it at any
real point, conserving essential properties such as smoothness. Unser describes
the power and the efficiency of Bspline interpolation, which allows to get
rid of all the matrician inversion calculated in conventional splineinterpolation.
Given an input signal, its Bspline description is easily obtained via simple
digital filtering, using IIR filters with fast recursive implementations.
Among all existing wavelets, that verify the admissibility conditions, Bspline wavelets have the advantage of being explicitly known, and of not depending on recursive definitions. The famous Haarwavelet is just the sum of two Bsplines of degree 0. Other wavelets, such as the first derivative of a Gaussian () or the second derivative (Mexican wavelet) can be closely approximated by Bsplines of higher degrees . The corresponding approximation errors of a Gaussian are negligible at small scales, and tend to when the Bspline degree grows. This expression (14) of wavelets into Bspline basis (through either orthogonal or oblique projections) allows to compute the inherent convolution products of the CWT in Bspline basis efficiently, taking advantage of properties such as compact support (convolutions are computed on finite and small number of points). According to all these considerations, we may now suppose that the mother wavelet is expressed in a Bspline basis of order :
3 Real CWT Computation in Bspline basis
General case
Using (10), and thanks to commutation properties of the convolution operator, the CWT becomes :
Expanding the modified finite differences operator in the previous results yields:
(with ), which is a particular case (finite sum) of a mixted convolution defined in (13). Bspline description of the analyzed function If we consider that our input signal belongs to a Bspline space of degree ,
where the interpolation coefficients were determined by , with representing the discrete samples of the Bspline of degree (Bspline kernel), and representing sampled at integer points. We can rewrite as :
thus : with and . Going back to the mixtedconvolution expression of the CWT, the previous relation yields :
where . Finally, we use the compact support property of Bsplines : vanishes outside . Consequently, the only nonnull coefficients of are given by:
By calling , we obtain the final expression of the CWT :
Comments According to (23),(25) and (26), the weights are equal to :
The 'big' number of dimensions (3) of this filter can be explained easily. The
index ''corresponds to the decomposition of the initial wavelet into
weighted and rescaled Bsplines ; corresponds to the decomposition
of each of these rescaled Bsplines. Finally, corresponds to the
convolution of Bsplines from the signal decomposition with the ones from the
rescaled wavelet description.
Usually, corresponds to the time locations of the sampled signals and
takes thus integer values. The shift by in cancels
out with the corresponding shift of That means that the values
of
remain the same for , and
fixed, whatever . This reduces considerably the number of weights since
only
need to be computed.
It is also interesting to study the localization of the weights. We first freeze
the index . One can observe that, for fixed, the weights are
modulated values of Bspline, uniformly distributed with a step of 1 (we call
it a 'convolution cluster'). Then, when we increment of decrement ,
values of the Bspline may or may not be equal, depending on the value of ,
but the central value and the lower bound are
simply translated by . Hence for fixed, the CWT consists in
filtering with 'clusters' of length each
cluster being separated from its neighbors by a distance . This can
be seen as a kind of modified 'a trous' filter, as shown in figure 2.
Moreover, when is incremented or decremented, the whole set of clusters
corresponding to a index is also shifted by since
and play the same role in (27).
Fast Implementation
We have developed a fast algorithm for the computation of the CWT at any real scale and any integer timelocalization by using equation (26) and the study of the filter As initialization step, we compute the values of the weights . This coefficients may be stored for precalculated cases (predefined wavelets and their Bspline expansion coefficients and expansion degree , interpolation degree and predefined values of the scales ). In a first step, the Bspline expansion coefficients of the sampled signal are calculated, and the finite differences operator is applied times. The intermediate result does not depend on the scale The second step, repeated for all the values of , consists in filtering with the stored values of . The dependence with regards to only concerns the values of chosen and the intercluster distance for the filtering, but the complexity per point is fixed by the values of and . It has to be noticed that the independence between scales allows straightforward parallel implementations. The algorithm is described in figure 3.
4 Complex CWT (GaborLike CWT, CCWT)
The Continuous Wavelet Transform we used so far supposed that the mother wavelet was a real function. Such a transform can be considered as a correlator, where the correlation pattern is rescaled by This method is very efficient for edge detection, singularities characterization (Hölder exponent), ... , but does not really give us a timefrequency description of the signal (in the sense of a Fourier Transform indicating which sinusoid frequencies are present in the signal). An efficient method for the characterization of the period of evolving signals is the Complex Continuous Wavelet Transform (GaborLike CWT). This CCWT is based on a complex wavelet, namely a complex exponential of frequency windowed by a Gaussian :
This wavelet transform is relatively similar to the Short Time Fourier Transform,
but since the gaussian window is a function of , it will grow in width
for large scale and hence keep a constant ratio between analyzing period and
concerned frequency, keeping a constant precision over scales. Unser
() described a fast implementation of this transform using Bsplines in order
to approximate the gaussian function. The key idea is that the central Bspline
converges to a Gaussian uniformly when tends to infinity.
For , the variance product is already within of the limit
specified by the uncertainty principle. In order to compute the Complex CWT,
Unser proposes a modulation approach : the input signal is first modulated by
a complex exponential
, then filtered through the window
at scale , and finally demodulated by an other exponential
, as shown in figure 4. At a scale the CCWT
is thus equivalent to the correlation of the input signal with a windowed complex
exponential whose frequency is that means, for an initial
frequency of (common case), the scale simply corresponds
to the period. Usually, the modulus of the CCWT is displayed in a scalogram,
a 2Dplot where the xaxis is the time , and the yaxis corresponds
to (Figure 4).
The link with our method is rather obvious. Excepted the pre and postmodulations by complex exponentials, we may consider the filtering window as a wavelet described in a Bspline basis, where the coefficients of (14) are one for and zero otherwise. The Gaussian is approximated by a Bspline of degree So we used the same algorithm, but replaced the real cwt method used at integer scales by ours in order to be able to compute it at any scale.

EPFL Swiss Federal Institute of Technology Lausanne BIG Biomedical Imaging Group Address: DMT/IOA, BMEcublens, CH1015 Lausanne, Switzerland 
Imaging Web Demonstration 
Contact: Webmaster 
21 June 2000ont> 