Forums

How to implement Continuous Wavelet Transform (CWT) using FFT

Started by tt- February 19, 2012
Hi all, 

I'm having some troubles implementing the Continuous Wavelet Transform
(CWT) 
using the Fast Fourier Transform (FFT).

I understand that the idea to do this is to first define

  psi~(a, t) = (1 / sqrt(a)) * conj(psi(-t / a))

where conj denotes a complex conjugate and psi is the original wavelet
function. Now the CWT can be expressed as the convolution of the original
signal and psi~. Then I can calculate the convolution using the convolution
theorem and FFT.

So the algorithm should like 

 for i = 1:numscales
     % take next scale
     a = scales(i);
     
     % calculate psi the but how?

     % convolution theorem
     X(i, :) = ifft(fft(signal) .* fft(psi));
 end;

The problem I'm having is how to calculate the values for the wavelet
function psi.

Suppose I want to use the Morlet wavelet, psi = exp(-t^2/2) * cos(5*t),
which is defined in the interval -4 <= t <= 4.

The first step is to replace the t with (-t /a) as suggested by the 
definition of psi~, so the equation becomes psi = exp(-(-t/a)^2/2) *
cos(5*(-t/a)).

This is as far as I get. If I do something like 
 t = linspace(-4, 4, 4096)
to calculate values of psi I'm getting the wrong results. 

How do I calculate the psi to get the code working?



On Feb 19, 12:47&#2013266080;pm, "tt-" <tuomas.p@n_o_s_p_a_m.gmail.com> wrote:
> Hi all, > > I'm having some troubles implementing the Continuous Wavelet Transform > (CWT) > using the Fast Fourier Transform (FFT). > > I understand that the idea to do this is to first define > > &#2013266080; psi~(a, t) = (1 / sqrt(a)) * conj(psi(-t / a)) > > where conj denotes a complex conjugate and psi is the original wavelet > function. Now the CWT can be expressed as the convolution of the original > signal and psi~. Then I can calculate the convolution using the convolution > theorem and FFT. > > So the algorithm should like > > &#2013266080;for i = 1:numscales > &#2013266080; &#2013266080; &#2013266080;% take next scale > &#2013266080; &#2013266080; &#2013266080;a = scales(i); > > &#2013266080; &#2013266080; &#2013266080;% calculate psi the but how? > > &#2013266080; &#2013266080; &#2013266080;% convolution theorem > &#2013266080; &#2013266080; &#2013266080;X(i, :) = ifft(fft(signal) .* fft(psi)); > &#2013266080;end; > > The problem I'm having is how to calculate the values for the wavelet > function psi. > > Suppose I want to use the Morlet wavelet, psi = exp(-t^2/2) * cos(5*t), > which is defined in the interval -4 <= t <= 4. > > The first step is to replace the t with (-t /a) as suggested by the > definition of psi~, so the equation becomes psi = exp(-(-t/a)^2/2) * > cos(5*(-t/a)). > > This is as far as I get. If I do something like > &#2013266080;t = linspace(-4, 4, 4096) > to calculate values of psi I'm getting the wrong results. > > How do I calculate the psi to get the code working?
Are you sure the convolution integral and the continuous wavelet integral are of the same form?
Yes, I'm quite sure. I used the material available at 
http://www.mathworks.se/help/toolbox/wavelet/gs/f3-1000759.html#bsu1pdf
for reference. So this is probably the same way it is implemented in
Matlab. 

>Are you sure the convolution integral and the continuous wavelet >integral are of the same form?
On Feb 20, 1:30&#2013266080;pm, "tt-" <tuomas.p@n_o_s_p_a_m.gmail.com> wrote:
> Yes, I'm quite sure. I used the material available at > http://www.mathworks.se/help/toolbox/wavelet/gs/f3-1000759.html#bsu1pdf > for reference. So this is probably the same way it is implemented in
Have you tried using the psi that Matlab uses? http://www.mathworks.com/help/toolbox/wavelet/ref/cwtftinfo.html Dale B, Dalrymple