Biomedical Imaging Group: B-Spline based Continous Wavelet Transform (description)

B-Spline based Continous Wavelet Transform

Introduction

Description

Example

1 Introduction

The Continuous Wavelet Transform (CWT) of a signal $f$ in the $\psi$ basis is defined as

$\begin{displaymath} W_{\psi ,f}(a,b)=\frac{1}{\sqrt{a}}\int ^{+\infty }_{-\infty... ...{\sqrt{a}}f(.),\psi \left( \frac{b-.}{a}\right) \right\rangle \end{displaymath}$

(1)

It can be interpreted as the correlation of the input signal with a mirrored version of $\psi _{a}$ , which is $\psi$ rescaled by $a$ . For a 1-D input, the result is a 2-D time-frequency description of the signal. The scale $a$ is proportionally inverse to the frequency studied, and $b$ represents the time localization at which we analyze the signal. The bigger the scale $a,$ the wider the analyzing function $\psi _{a}$ , 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 $a\, =\, 2^{i}$ and time-shifts $b=2^{i}k,\, k\in Z$ . Using Mallat's algorithm or derivatives, they get an overall $O(N)$ complexity. Other techniques compute the Wavelet Transform at dyadic scales and at linear time-points, or, by using $M$ mother wavelets closely dilated, at M-adic scales. Other methods compute the CWT at integer scales ( $a=k$ ).

We present a new method, using expansion of both the input signal and the mother wavelet in B-Spline 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 B-spline properties

Basic Definitions

The B-spline (``B'' standing for ``Basic'') $\beta ^{0}$ of degree $0$ is defined as

$\begin{displaymath} \beta ^{0}(x)=\left\{ \begin{array}{c} 1\, \, if\, \, -0.5<x\leq 0.5\\ 0\, \, otherwise \end{array}\right. \end{displaymath}$

(2)

The B-spline $\beta ^{n}$ of degree $n$ is defined as the $n-th$ convolution of $\beta ^{0}$ with itself :

$\begin{displaymath} \beta ^{n}(x)=\left( \underbrace{\beta ^{0}*\beta ^{0}*...*\beta ^{0}}_{n\, times}\right) (x) \end{displaymath}$

(3)

Now, we enumerate a series of definitions that are necessary in our formalism :

Continuous Operators :

Derivative :

$\begin{displaymath} Df(x)=\frac{df(x)}{dx} \end{displaymath}$ (4)
Integral :

$\begin{displaymath} D^{-1}f(x)=\int _{-\infty }^{x}f(t)dt \end{displaymath}$ (5)
One-sided power function :

$\begin{displaymath} x_{+}^{n}=\left\{ ^{x^{n}\, if\, x\geq 0}_{0\; otherwise}\right. \end{displaymath}$ (6)

which is equivalent to the (n+1)-fold integral operator since :

$\begin{displaymath} D^{-(n+1)}f(x)=\frac{x_{+}^{n}}{n!}*f(x) \end{displaymath}$ (7)

Discrete Operators

Finite differences :

$\begin{displaymath} \Delta =\delta (x)-\delta (x-1) \end{displaymath}$ (8)
Inverse finite differences operator :

$\displaystyle \Delta ^{-1}$ $\textstyle =$ $\displaystyle \sum _{n\geq 0}\delta (x-n)$ (9)

Using these definitions, a B-spline $\beta ^{n}$ of degree $n$ can be written as:

$\begin{displaymath} \beta ^{n}(x)=\left( \Delta ^{n+1}*\frac{x_{+}^{n}}{n!}*\delta \left( \cdot +\frac{n+1}{2}\right) \right) (x) \end{displaymath}$

(10)

Intuitively, this means that the convolution with a B-spline of degree $n$ is equivalent to shift the signal by $\frac{n+1}{2}$ , to integrate it $n$ times, and to apply $\left( n+1\right)$ times a finite differences operator to the result of the integration.
Using the exchange rules for the one-sided power function and for the shift with the rescaling operator given in , and noting that the rescaled finite differences operator $\Delta ^{n}_{a}$ can be expressed as

$\begin{displaymath} \Delta _{a}^{n}\left( x\right) =\left\vert a\right\vert \sum... ...c} n\\ k \end{array}\right) (-1)^{k}}_{q_{n}(k)}\delta (x-ak) \end{displaymath}$

(11)

while a B-spline of order $n$ and rescaled by a scale $a$ writes :

$\displaystyle \beta _{a}^{n}(x)$

$\textstyle =$

$\displaystyle \left( \Delta _{a}^{n+1}*\frac{1}{a^{n+1}}\frac{x_{+}^{n}}{n!}*\delta \left( \cdot +a\frac{\left( n+1\right) }{2}\right) \right) (x)$

(12)

Finally, we define the mixted convolution of a continuous signal $v(x)$ with a discrete one $q(k)$ with intersample distance equals to $a$ as:

$\displaystyle g(x)$

$\textstyle =$

$\displaystyle v(x)*\sum _{k\in Z}q(k)\delta (x-ak)$

(13)

**Figure 1:** Graphical Interpretation of the mixted convolution : weighted sum of shifted versions of the function v(x)
$\begin{figure} \setlength {\unitlength}{1973sp}%\par\begingroup\makeatletter\i... ...nt{6}{7.2}{\rmdefault}{\mddefault}{\updefault}=}}} \end{picture}\par\end{figure}$

B-spline 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 B-spline interpolation, which allows to get rid of all the matrician inversion calculated in conventional spline-interpolation. Given an input signal, its B-spline description is easily obtained via simple digital filtering, using IIR filters with fast recursive implementations.

B-spline Wavelets

Among all existing wavelets, that verify the admissibility conditions, B-spline wavelets have the advantage of being explicitly known, and of not depending on recursive definitions. The famous Haar-wavelet is just the sum of two B-splines of degree 0. Other wavelets, such as the first derivative of a Gaussian () or the second derivative (Mexican wavelet) can be closely approximated by B-splines of higher degrees . The corresponding approximation errors of a Gaussian are negligible at small scales, and tend to $0$ when the B-spline degree grows.

This expression (14) of wavelets into B-spline basis (through either orthogonal or oblique projections) allows to compute the inherent convolution products of the CWT in B-spline 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 $\psi$ is expressed in a B-spline basis of order $n_{1}$ :

$\begin{displaymath} \psi (x)=\sum ^{K}_{k=-K}d_{k}\beta ^{n_{1}}(x-k) \end{displaymath}$

(14)

3 Real CWT Computation in B-spline basis

General case

Using (10), and thanks to commutation properties of the convolution operator, the CWT $W_{\psi ,f}(a,b)=\frac{1}{\sqrt{a}}\left( f(\cdot )*\psi _{a}(\cdot )\right) (b)$ becomes :

$\begin{displaymath} W_{\psi ,f}(a,b)=\left( \frac{1}{a^{n_{1}+\frac{3}{2}}}\sum ... ... \cdot +a\left( \frac{n_{1}+1}{2}-k\right) \right) \right) (b) \end{displaymath}$

(15)

This shows that the wavelet transform can be seen, for each weight $d_{k}$ , as the application of a modified finite differences operator to the $\left( n_{1}+1\right)$ -th integral of a shifted version of the analyzed function $f$ .

Expanding the modified finite differences operator $\Delta ^{n_{1}+1}_{a}$ in the previous results yields:

$\displaystyle W_{\psi ,f}(a,b)$	$\textstyle =$	$\displaystyle \underbrace{A(a)\sum ^{K}_{k=-K}d_{k}D^{-(n_{1}+1)}f\left( \cdot ... ...) \right) }_{v(x)}*\sum ^{n_{1}+1}_{l=0}q_{n_{1}+1}(l)\delta \left( b-al\right)$	(16)
	$\textstyle =$	$\displaystyle v(x)*\sum ^{n_{1}+1}_{l=0}q_{n_{1}+1}(l)\delta (b-al)$	(17)

(with $A(a)=\frac{\left\vert a\right\vert }{a^{n_{1}+\frac{3}{2}}}$ ), which is a particular case (finite sum) of a mixted convolution defined in (13).

B-spline description of the analyzed function $f$

If we consider that our input signal $f(x)$ belongs to a B-spline space of degree $n_{2}$ ,

$\displaystyle f(x)$	$\textstyle =$	$\displaystyle \sum _{i\in Z}c_{f}(i)\cdot \beta ^{n_{2}}(x-i)$	(18)
	$\textstyle =$	$\displaystyle c_{f}*\beta ^{n_{2}}(x)$	(19)

where the interpolation coefficients $c_{f}$ were determined by $c_{f}=\left( b^{n_{2}}_{1}\right) ^{-1}*f_{\delta }$ , with $\left( b^{n_{2}}_{1}\right)$ representing the discrete samples of the B-spline of degree $n_{2}$ (B-spline kernel), and $f_{\delta }$ representing $f$ sampled at integer points.
We can rewrite $v(x)$ as :

$\begin{displaymath} v(x)=A(a)c_{f}*\left[ \sum ^{K}_{k=-K}d_{k}D^{-(n_{1}+1)}\be... ...}}\left( x+a\left( \frac{n_{1}+1}{2}-k\right) \right) \right] \end{displaymath}$

(20)

But :

$\displaystyle D^{-\left( n_{1}+1\right) }\beta ^{n_{2}}(x)$

$\textstyle =$

$\displaystyle \Delta ^{-\left( n_{1}+1\right) }*\beta ^{\left( n_{1}+n_{2}+1\right) }\left( x-\frac{n_{1}+1}{2}\right)$

(21)

thus :

$\displaystyle v(x)$	$\textstyle =$	$\displaystyle A(a)c_{f}\left[ \sum ^{K}_{k=-K}d_{k}\Delta ^{-(n_{1}+1)}\beta ... ..._{2}+1\right) }\left( x+(a-1)\left( \frac{n_{1}+1}{2}\right) -ak\right) \right]$	(22)
	$\textstyle =$	$\displaystyle A(a)\underbrace{\Delta ^{-(n_{1}+1)}c_{f}}_{C}\sum ^{K}_{k=-K}d_{k}\beta ^{n}\left( x+\tau _{k}\right)$	(23)

with $\tau _{k}=\left( \left( a-1\right) \left( \frac{n_{1}+1}{2}\right) -ak\right)$ and $n=\left( n_{1}+n_{2}+1\right)$ .

Going back to the mixted-convolution expression of the CWT, the previous relation yields :

$\displaystyle W_{\psi ,f}(a,b)$	$\textstyle =$	$\displaystyle A(a)C\left[ \sum ^{K}_{k=-K}d_{k}\beta ^{n}\left( x+\tau _{k}\right) \right] \sum ^{n_{1}+1}_{l=0}q_{n_{1}+1}(l)\delta (b-al)$	(24)
	$\textstyle =$	$\displaystyle \sum ^{K}_{k=-K}\sum ^{n_{1}+1}_{l=0}\sum _{p\in Z}A(a)d_{k}q_{n_{1}+1}(l)C(p)\beta ^{n}\left( \tau _{k,l}'-p\right)$	(25)

where $\tau _{k,l}'=b-al+\tau _{k}$ .

Finally, we use the compact support property of B-splines : $\beta ^{n}$ vanishes outside $\left[ \frac{n+1}{2},\frac{n+1}{2}\right)$ . Consequently, the only non-null coefficients of $\beta ^{n}\left( \tau _{k,l}'-p\right)$ are given by:

$\begin{eqnarray*} -\frac{n+1}{2}\leq & \tau _{k,l}'-p & <\frac{n+1}{2} \end{eqnarray*}$

By calling $p_{o}=\left\lceil \tau _{k,l}'-\frac{n+1}{2}\right\rceil$ , we obtain the final expression of the CWT :

$\begin{displaymath} W_{\psi ,f}(a,b)=\sum ^{K}_{k=-K}\sum ^{n_{1}+1}_{l=0}\sum ^... ...1}(l)\beta ^{n}\left( \tau _{k,l}'-p\right) }_{w_{a,b}(k,l,p)} \end{displaymath}$

(26)

Comments

According to (23),(25) and (26), the weights $w_{a,b}(k,l,p)$ are equal to :

$\begin{displaymath} w_{a,b}(k,l,m)=A(a)d_{k}q_{n_{1}+1}(l)\beta ^{n}\left( b-p-a... ...) +\left( a-1\right) \left( \frac{n_{1}+1}{2}\right) \right) , \end{displaymath}$

(27)

where $p$ is an integer value between $0$ and $n.$

The 'big' number of dimensions (3) of this filter can be explained easily. The index ' $k$ 'corresponds to the decomposition of the initial wavelet into $K$ weighted and rescaled B-splines ; $'l'$ corresponds to the decomposition of each of these rescaled B-splines. Finally, $'p'$ corresponds to the convolution of B-splines from the signal decomposition with the ones from the rescaled wavelet description.

Usually, $b$ corresponds to the time locations of the sampled signals and takes thus integer values. The shift by $b$ in $\tau _{k,l}'$ cancels out with the corresponding shift of $p_{0}.$ That means that the values of $w_{a,b}(k,l,p)$ remain the same for $a$ , $k$ and $l$ fixed, whatever $b$ . This reduces considerably the number of weights since only $w_{a,0}(k,l,p)$ need to be computed.

It is also interesting to study the localization of the weights. We first freeze the index $k$ . One can observe that, for $l$ fixed, the weights are modulated values of B-spline, uniformly distributed with a step of 1 (we call it a 'convolution cluster'). Then, when we increment of decrement $l$ , values of the B-spline may or may not be equal, depending on the value of $a$ , but the central value $\tau _{k,l}'$ and the lower bound $p_{0}$ are simply translated by $a$ . Hence for $k$ fixed, the CWT consists in filtering $C$ with $n_{1}+1$ 'clusters' of length $(n+1),$ each cluster being separated from its neighbors by a distance $a$ . This can be seen as a kind of modified 'a trous' filter, as shown in figure 2. Moreover, when $k$ is incremented or decremented, the whole set of clusters corresponding to a $k$ -index is also shifted by $a,$ since $k$ and $l$ play the same role in (27).

**Figure 2:** Localization of the weights w(a,k,l,m) for k fixed, a=30, n1=3, n2=0
$\resizebox* {1\columnwidth}{!}{\rotatebox{90}{\includegraphics{WPosition.eps}}}$

Fast Implementation

We have developed a fast algorithm for the computation of the CWT at any real scale $a$ and any integer time-localization $b$ by using equation (26) and the study of the filter $w.$

As initialization step, we compute the values of the weights $w_{a}(k,l,m)$ . This coefficients may be stored for pre-calculated cases (pre-defined wavelets $\psi$ and their B-spline expansion coefficients $d_{k}$ and expansion degree $n_{1}$ , interpolation degree $n_{2}$ and pre-defined values of the scales $a$ ).

In a first step, the B-spline expansion coefficients $c_{f}$ of the sampled signal $f_{\delta }$ are calculated, and the finite differences operator $\Delta ^{-1}$ is applied $\left( n_{1}+1\right)$ times. The intermediate result $C$ does not depend on the scale $a.$

The second step, repeated for all the values of $a$ , consists in filtering $C$ with the stored values of $w$ . The dependence with regards to $a$ only concerns the values of $w$ chosen and the inter-cluster distance for the filtering, but the complexity per point is fixed by the values of $K,$ $n_{1}$ and $n_{2}$ . It has to be noticed that the independence between scales allows straightforward parallel implementations. The algorithm is described in figure 3.

**Figure 3:** Fast implementation of the fast CWT. (a) : Initialization, (b) : scale dependent step
$\begin{picture}(300,180) \put(0,140){(a)} \put(100,125){$ f_{\delta }(k) $} ... ...r(1,0){40}} \put(230,25){$ W_{\psi ,f}(a,k) $} \put(0,40){(b)} \end{picture}$

4 Complex CWT (Gabor-Like CWT, CCWT)

The Continuous Wavelet Transform we used so far supposed that the mother wavelet $\psi$ was a real function. Such a transform can be considered as a correlator, where the correlation pattern is $\psi$ rescaled by $a.$ This method is very efficient for edge detection, singularities characterization (Hölder exponent), ... , but does not really give us a time-frequency 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 (Gabor-Like CWT). This CCWT is based on a complex wavelet, namely a complex exponential of frequency $f$ windowed by a Gaussian $w$ :

$\begin{displaymath} \psi (x)=w(x)e^{-j2\pi fx} \end{displaymath}$

(28)

This wavelet transform is relatively similar to the Short Time Fourier Transform, but since the gaussian window is a function of $x$ , 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 B-splines in order to approximate the gaussian function. The key idea is that the central B-spline $\beta ^{n}$ converges to a Gaussian uniformly when $n$ tends to infinity. For $n=3$ , the variance product is already within $0.5\%$ 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 $e^{-j2\pi fx/a}$ , then filtered through the window $w_{a}$ at scale $a$ , and finally demodulated by an other exponential $e^{j2\pi fx/a}$ , as shown in figure 4. At a scale $a,$ the CCWT is thus equivalent to the correlation of the input signal with a windowed complex exponential whose frequency is $\frac{f}{a},$ that means, for an initial frequency $f$ of $1$ (common case), the scale $a$ simply corresponds to the period. Usually, the modulus of the CCWT is displayed in a scalogram, a 2D-plot where the x-axis is the time $b$ , and the y-axis corresponds to $a$ (Figure 4).

**Figure 4:** Computation of the Complex CWT by modulation, gaussian-windowing and demodulation

The link with our method is rather obvious. Excepted the pre- and post-modulations by complex exponentials, we may consider the filtering window as a wavelet $\psi$ described in a B-spline basis, where the $d_{k}$ coefficients of (14) are one for $k=0$ and zero otherwise. The Gaussian is approximated by a B-spline of degree $n\geq 3.$ 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, BM-Ecublens, CH-1015 Lausanne, Switzerland

Imaging Web Demonstration

Contact: Webmaster

21 June 2000ont>