function CWT = cwt(s,n_voice,wavelet)

% CWT=cwt(s,n_voice,wavelet);
% Compute continuous wavelet transforms 
% of a signal whose length is n=2^p. 
% The wavelets used:
%    "Morlet", "Gabor", "Gauss", 
%     "Sombrero#", "Spline#".
% CWT      matrix n by n_scales 
%             n_scales = n_voice .* n_octaves
% The values of the scale parameter are such that
% log_2(a) is included between 2 and p.  


b=min(6,length(wavelet));
b1=min(8,length(wavelet));
x=s(:)';
n=length(x);
xhat=fft(x);
xi   = [ (0: (n/2)) (((-n/2)+1):-1) ] .* (2*pi);
	Gauss    = 0;
	Sombrero = 0;
	Morlet   = 0;
	Gabor    = 0;
    Spline   = 0;
     
if	strcmp(wavelet,'Gauss'),
		Gauss    = 1;
elseif	strcmp(wavelet(1:b1),'Sombrero'),
		Sombrero = 1;
		enne = str2num(wavelet(9:length(wavelet)));
elseif	strcmp(wavelet,'Morlet'),
		Morlet   = 1;	
		omega0  = pi*sqrt(2/log(2)); 
elseif  strcmp(wavelet(1:b),'Spline'),
		Spline	 = 1;
		enne = str2num(wavelet(7:length(wavelet)));
		a = 0.697066; sigma = 0.561145;
		omega0 = 0.409177/(2*pi);
		beta   = 2*sigma*(enne+1);
elseif	strcmp(wavelet,'Gabor'),
		Gabor	 = 1;  w0=5.33644;
else error('Use cwt: Gabor|Gauss|Morlet|Sombrero#|Spline#|');
	end	

n_octaves=floor(log2(n))-1;
n_scales= n_voice .* n_octaves;
	
k_scale=1;
scale=4;
CWT = zeros(n,n_scales);


for j=1:n_octaves,
	    for l=0:n_voice-1,
		   q_scale=scale .* (2^(l/n_voice));
		   omega=xi ./q_scale;


if     Gauss,
  psi_hat=exp(-omega.^2 ./2);
elseif Sombrero,
  psi_hat=((omega/2).^enne) .* exp(-(omega/2).^2 ./2);
elseif Morlet,
  psi_hat=exp(-(omega - omega0).^2 ./2) - exp(-(omega.^2 + omega0.^2)/2);
elseif Spline
  psi_hat=real(((1-exp(-i*(omega+.01)/16))./(i*(omega+.01)/16)).^enne);
elseif Gabor
  psi_hat=exp(-(omega-w0).^2 ./2);
end
w_hat=xhat.*psi_hat; 
w=ifft(w_hat);	
CWT(1:n,k_scale)=(w)';
k_scale=k_scale+1;
        end

        scale=scale .*2;
end 



% "Complex and Chaotic Nonlinear Dynamics. 
% Advances in Economics and Finance, 
% Mathematics and Statistics" 
% T.Vialar, Springer. 2009.
% Copyright(c).