The Wavelet transform using the Fourier transform

We start with the set of scalar products $c_0(k)=<f(x),\phi(x-k)>$. If $\phi(x)$ has a cut-off frequency $\nu_c\le {1\over 2}$[#starck1<#14704,#starck2<#14705,#starck3<#14706,#starck4<#14707], the data are correctly sampled. The data at the resolution j=1 are:

$$c_1(k)=<f(x),\frac{1}{2}\phi(\frac{x}{2}-k)>$$ (14.47)

and we can compute the set c₁(k) from c₀(k) with a discrete filter $\hat h(\nu)$:

$$\hat h(\nu)= \left\{ \begin{array}{ll} {\hat{\phi}(2\nu)\over \ha... ...u_c \\ 0 & \mbox{if } \nu_c \leq \mid \nu \mid < {1\over 2} \end{array}\right.$$ (14.48)

and

$$\forall \nu, \forall n \mbox{ } \hat h(\nu + n) = \hat h(\nu)$$ (14.49)

where n is an integer. So:

$$\hat{c}_{j+1}(\nu)=\hat{c}_{j}(\nu)\hat{h}(2^{j}\nu)$$ (14.50)

The cut-off frequency is reduced by a factor 2 at each step, allowing a reduction of the number of samples by this factor.

The wavelet coefficients at the scale j+1 are:

$$w_{j+1}(k)=<f(x),2^{-(j+1)}\psi(2^{-(j+1)}x-k)>$$ (14.51)

and they can be computed directly from c_j(k) by:

$$\hat{w}_{j+1}(\nu)=\hat{c}_{j}(\nu)\hat g(2^{j}\nu)$$ (14.52)

where g is the following discrete filter:

$$\hat g(\nu)= \left\{ \begin{array}{ll} {\hat{\psi}(2\nu)\over \ha... ...u_c \\ 1 & \mbox{if } \nu_c \leq \mid \nu \mid < {1\over 2} \end{array}\right.$$ (14.53)

and

$$\forall \nu, \forall n \mbox{ } \hat g(\nu + n) = \hat g(\nu)$$ (14.54)

The frequency band is also reduced by a factor 2 at each step. Applying the sampling theorem, we can build a pyramid of $N+{N\over 2}+\ldots+1=2N$ elements. For an image analysis the number of elements is ${4\over 3}N^2$. The overdetermination is not very high.

The B-spline functions are compact in this directe space. They correspond to the autoconvolution of a square function. In the Fourier space we have:

$$\hat B_l(\nu)={\sin\pi\nu\over\pi\nu}^{l+1}$$ (14.55)

B₃(x) is a set of 4 polynomials of degree 3. We choose the scaling function $\phi(\nu)$ which has a B₃(x) profile in the Fourier space:

$$\hat{\phi}(\nu)={3\over 2}B_3(4\nu)$$ (14.56)

In the direct space we get:

$$\phi(x)={3\over 8}[{\sin{\pi x\over 4}\over {\pi x\over 4}}]^4$$ (14.57)

This function is quite similar to a Gaussian one and converges rapidly to 0. For 2-D the scaling function is defined by $\hat \phi(u,v) = {3\over 2}B_3(4r)$, with $r = \sqrt(u^2+v^2)$. It is an isotropic function.

The wavelet transform algorithm with n_p scales is the following one:

1.

We start with a B3-Spline scaling function and we derive $\psi$, h and g numerically.

2.

We compute the corresponding image FFT. We name T₀ the resulting complex array;

3.

We set j to 0. We iterate:

4.

We multiply T_j by $\hat g(2^ju,2^jv)$. We get the complex array W_j+1. The inverse FFT gives the wavelet coefficients at the scale 2^j;

5.

We multiply T_j by $\hat h(2^ju,2^jv)$. We get the array T_j+1. Its inverse FFT gives the image at the scale 2^j+1. The frequency band is reduced by a factor 2.

6.

We increment j

7.

If $j \leq n_p$, we go back to 4.

8.

The set $\{w_1, w_2, \dots, w_{n_p}, c_{n_p}\}$ describes the wavelet transform.

If the wavelet is the difference between two resolutions, we have:

$$\hat \psi(2\nu) = \hat \phi(\nu) - \hat \phi(2\nu)$$ (14.58)

and:

$$\hat g(\nu) = 1 - \hat h(\nu)$$ (14.59)

then the wavelet coefficients $\hat w_j(\nu)$ can be computed by $\hat c_{j-1}(\nu) - \hat c_j(\nu)$.