3
$\begingroup$

I want to write continuous wavelet transform codes manually by matlab. And I want to use complex morlet function. Here are some background:
Continuous wavelet transform definition: $C(S,T;f(t),\psi (t))=\frac{1}{\sqrt{S}}\int_{S}^{b}f(t)\psi^{\ast }(\frac{t-T}{S})$
S is scale vector.for example 1:60;$T$ is time slid. And $\psi(t)= \frac{1}{\sqrt{\pi f_{b}}}e^{2\pi f_{c} t}e^{\frac{-t^{2}}{f_{b}}}$

psi = ((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*t).*exp(-t.^2/fb); % for example fb=15;fc=1;

my discrete signal has $N$ points. discrete version of that integral for the first element in the scale vector is: for s(i)
$C_{T}(s)=\sum_{n=0}^{N-1}f(n)\psi^{*}(\frac{T-n}{S})$ this must calculate for all scales,1:60;
at the end it must return me a complex matrix with $N\times S$ size. I am confused about wrting this codes. Anyone can help?
P.S I don't want to uses matlab function conv2 to calculate that convolution.
Thanks in advance. here is my first attempt, but it's not working at all.

%% user CWT
clear all
N=300;                %sample point numbers
t=linspace(0,30,N);
%% signal
x=5*sin(2*pi*0.5*t);  % signal with freq of 0.5 HZ
%% cwt
    fc=1;fb=15;
% psi=((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*...
%     t).*exp(-t.^2/fb);
%% convolution Psi([N-n]/S)*x(n) so we calculate convolution(psi(n/s),x(n))
    for s=1:60   %scale vector s=[1:1:60]
for i = 1:N      % number of discrete times
       for k = 1:i 
               if ((i-k+1)<N+1) && (k <N+1)
                PSI = ((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*...
       (t/s)).*exp(-(t/s).^2/fb);   
                        c(i,s) = c(i)+ x(k)*PSI(i-k+1);
                end
       end
end
    end
$\endgroup$
2
$\begingroup$

Here is the Python code that I use for CWT. You'll need a short 16bit mono 44.1Khz test.wav soundfile as input.

In order to compute the CWT, we need to compute the convolution between the input x[n] and the morlet wavelet. An efficient way to do this is to do ifft(x_ft * morlet_ft) ( I understand this as $f * g(t) =\frac{1}{2 \pi} \int_{-\infty}^{+\infty} \hat{f}(w) \hat{g}(w) e^{i t w} d w$ : it is more efficient to do a multiplication of the fourier transforms and THEN to do an inverse fourier transform, than to compute a convolution directly)

from scipy.io.wavfile import read
from pylab import *
import matplotlib.pyplot as plt
import numpy as np
import scipy

# read the file                              
(sr, samples) = read('test.wav')    
x = (np.float32(samples)/((2 ** 15)-1))
N = len(x)
x_fft = np.fft.fft(x)

# scales
J = 200    
scales = np.asarray([2**(i * 0.1) for i in range(J)])

# pre-compute the Fourier transform of the Morlet wavelet 
morletft = np.zeros((J, N))
for i in range(J):
  morletft[i][:N/2] = sqrt(2 * pi * scales[i]) * exp(-(scales[i] * 2 * pi * scipy.array(range(N/2)) / N - 2)**2 / 2.0)

# compute the CWT 
X = empty((J, N), dtype=complex128)
for i in range(J):
  X[i] = np.fft.ifft(x_fft * morletft[i])

# plot
plt.imshow(abs(X[:,scipy.arange(0,N,100)]), interpolation='none', aspect='auto')
plt.show()

Here is the scalogram plot that this code gives :

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for the code, I know ifft convolution is efficient and the computational burden is lower, but do you have any idea how can I get the convolution with normal definition? tnx much @basj $\endgroup$ – Electricman Nov 25 '13 at 21:44
  • $\begingroup$ If x[n] is defined for n=0,1,2,..N-1 you can use cwt[a,b]=sum(x*psi[(here scaled and shifted)]) , doesn't it work ? ;) $\endgroup$ – Basj Nov 25 '13 at 22:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.