% Continuous wavelet transform 
%
%
% WPsif(a,b) = \int f(x) 1/\sqrt{|a|} Psi((x-b)/a) dx 
%
% Radu Balan, 16 April 2012

% Define the wavelet
% Psi = -d2/dx2 exp(-pi x^2) = d/dx [2 pi x exp(-pi x^2)] = 2pi (1-2pi x^2)exp(-pi x^2)
%
deltax = 0.001;
x = (-10:deltax:10);
Psi00 = 2*pi*(1-2*pi*x.*x).*exp(-pi*x.*x);

% Signal of interest
alpha = 2;
f = exp(2*pi*i*alpha*x);
logav = -3:0.01:1.1;
bv = -2:0.01:2;
WPsif = zeros(length(bv),length(logav));

for ia = 1:length(logav)
    loga = logav(ia);
    a(ia) = 10^(loga);
    for ib = 1:length(bv);
        b = bv(ib);
        y = (x-b)/a(ia);
        Psiab = 2*pi*(1-2*pi*y.*y).*exp(-pi*y.*y);
        WPsif(ib,ia) = (1/sqrt(a(ia)))*sum(f.*Psiab)*deltax;
    end
end

figure(1)
plot(x,Psi00,x,real(f),x,imag(f))
legend('\Psi','real(f)','imag(f)')

figure(2)
mesh(a,bv,abs(WPsif))
xlabel('a'),ylabel('b')

figure(3)
ix = ceil(length(bv)/2.231);
semilogx(a,real(WPsif(ix,1:length(a))),a,imag(WPsif(ix,1:length(a))),a,abs(WPsif(ix,1:length(a))))
xlabel('a'),ylabel('W\Psi f')
legend('real','imag','abs')
title(sprintf('Continuous Wavelet Transform at b = %f',bv(ix)))

% Find the maximum
[WPsimax,imax] = max(abs(WPsif(ix,:)));
amax = a(imax);
ac = sqrt(5/(4*pi))/alpha;
fprintf(1,'Maximum occurs at: amax = %f , WPsif max = %f \n',amax,abs(WPsif(imax)));
fprintf(1,'Theoretical value: ac = %f \n',ac);
