      function [A1, A2, A3]= gaussq(s,n_voix,q)

%              [A1,A2,A3]=gaussq(x,n_voix,q)
%
%Cette fonction restitue la transformée en ondelettes continue
%           (a,b) ---> (1/sqrt (a)) * cwt_psi [S] (a,b)    
%d'un signal analogique réel S discrétisé en un signal digital de longueur
%n=2^p (donné en colonne) lorsque l'ondelette est l'atome psi(q ;.), dérivé
%à l'ordre q de l'atome Gaussien. En sortie, on obtient ainsi une matrice
%A3 à n colonnes et nb_voix * (p-1) lignes ; les valeurs du paramètre a
%d'échelle sont telles que log_2(a) varie entre 2 et p, chacune des p-1
%octaves dyadiques étant divisée en nb_voix canaux. Les lignes de A1 sont
%déduites des lignes de A3 en interpolant les valeurs prises sur chaque
%ligne aux maxima locaux de |A3| ; l'image positive A2 affiche, ligne par
%ligne, les valeurs aux maxima locaux de l'image |A3| sur la ligne correspondante.  

nn    = length(s);
shat = fft(s');
r=mod(q,2);
xi   = [ (0: (nn/2)) (((-nn/2)+1):-1) ] .* (2*pi); 
n_octaves=floor(log2(nn))-1; 	
n_echelles=n_voix.* n_octaves;
k_echelle=1;
echelle=4;
B1= zeros(nn,n_echelles);
B2= zeros(nn,n_echelles);
B3=zeros(nn,n_echelles);
for j = 1:n_octaves,
	    for l = 0:n_voix-1,
		   q_echelle = echelle.* (2^(l/n_voix));
		   omega  = xi ./ q_echelle ;
                   psi_hat=(omega.^q).* exp(-(omega.^2) ./2);
		   w_hat = psi_hat .* shat;
		   w    = ifft(w_hat);
                   if (r==0) 
                      u1=locmax(real(w));
                      u3=real(w);
                   else
                      u1=locmax(imag(w));
                      u3=imag(w);
                   end
                      u2=localmax(abs(w));	
                   B1(1:nn,k_echelle)=u1';
                   B2(1:nn,k_echelle)=u2';
                   B3(1:nn,k_echelle)=u3';
		   k_echelle = k_echelle+1;
	    end
	    echelle  = echelle *2;
    	end 
 B2(nn,:)=zeros(1,n_echelles);
 B11=interpol1(B1);
 A1=B11';
 A2=B2';
 A3=B3'; 