inquiries for writing continuous wavelet transform codes manually

Question

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

Answer 1

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 :