function [wave,period] = wavelet(x,dt);
%WAVELET  Wavelet transform
%       [wave,period] = wavelet(x,dt)
%       Computes the wavelet transform of the vector x (length n),
%       with sampling rate dt.
%       wave is the complex wavelet transform, na-by-n complex matrix,
%       where na=4*log(n*dt)/log(2dt) is the number of scales. 
%       period is a vector (length na) of the Fourier periods corresponding
%       to the chosen scales.
%
%       By default, the Morlet wavelet (k0=6) is used.
%       The wavelet basis is normalized to have total energy=1 at all scales.

f = fft(x);
n = length(x);

k = [1:fix(n/2)];
k = k.*((2.*pi)/(n*dt));
k = [0., k, -k(fix((n-1)/2):-1:1)];

a0 = 2.*dt;
oct = fix(log(n*dt/a0)/log(2)+0.999);
voice = 4;
na = oct*voice + 1;
ascale = a0*2.^((0:na-1)/voice);
period = ascale;

wave = [];

% loop through all scales and compute transform
for a1 = 1:na

%  Morlet wavelet
        a = ascale(a1);
        k0 = 6.;
        expnt = -(a.*k - k0).^2/2.*(k > 0.);
        factor = sqrt(a*k(2)*2.*sqrt(pi));
        morlet = factor/sqrt(2*pi)*exp(expnt);
        morlet = morlet.*(k > 0.);
        period(a1) = period(a1).*(4*pi)/(k0 + sqrt(2 + k0^2));

%  Wavelet transform
        wave  = [wave ; ifft(f.*morlet)];

end


