/* @(#)wave1d.c	17.1.1.1 (ES0-DMD) 01/25/02 17:21:23 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/******************************************************************************
**                   Copyright (C) 1993 ESO
*******************************************************************************
**
**    UNIT
**
**    Version: 17.1.1.1
**
**    Author: Jean-Luc Starck
**
**    Date:  02/01/25
**    
**    File:  wave1d.c
**
*****************************************************************************
**
**    DESCRIPTION 1D wavelet transform and reconstruction
**    -----------    
**   
*****************************************************************************
** 
** wave_1d_mex (Signal, W_1D, N, Nbr_Voie, Nbr_Plan, Scale_0)
** float ***W_1D, *Signal, *Scale_0;
** int N, Nbr_Voie, *Nbr_Plan;
**
** wavelet transform by the mexican hat
** Signal = IN: data
** W_1D = OUT: wavelet transform
** N = IN: number of data in the Signal
** Nbr_Voie = IN: number of channels per octave (in general, we take 12)
** Nbr_Plan = OUT: number of scales
** Scale_0 = OUT: first scale
**
** the wavelet transform W_1D[0..Nbr_Plan-1][0..N] is computed
** W_1D[i] defines the scale i 
**   
****************************************************************************
**
** wave_1d_mex_rec (W_1D, Signal, N, Nbr_Voie, Nbr_Plan)
** float **W_1D, *Signal;
** int N, Nbr_Voie, Nbr_Plan;
**
** wavelet reconstruction from a mexican hat wavelet transform
**
** W_1D = IN: wavelet transform
** Signal = OUT: recosntructed signal
** N = IN: number of data in the Signal
** Nbr_Voie = IN: number of channels per octave (in general, we take 12)
** Nbr_Plan = IN: number of scales
**   
**************************************************************************** 
**   
** wave_1d_french (Signal, W_1D, N, Nbr_Voie, Nbr_Plan, Scale_0)
** float ***W_1D, *Signal, *Scale_0;
** int N, Nbr_Voie, *Nbr_Plan;
**
** wavelet transform by the french hat
**
** Signal = IN: data
** W_1D = OUT: wavelet transform
** N = IN: number of data in the Signal
** Nbr_Voie = IN: number of channels per octave (in general, we take 12)
** Nbr_Plan = OUT: number of scales
** Scale_0 = OUT: first scale
**
** the wavelet transform W_1D[0..Nbr_Plan-1][0..N] is computed
** W_1D[i] defines the scale i 
**
**************************************************************************** 
**   
** wave_1d_french_rec (W_1D, Signal, N, Nbr_Voie, Nbr_Plan)
** float **W_1D, *Signal;
** int N, Nbr_Voie, Nbr_Plan;
**
** reconstruction from the wavelet transform with the french hat 
**
** W_1D = IN: wavelet transform
** Signal = OUT: recosntructed signal
** N = IN: number of data in the Signal
** Nbr_Voie = IN: number of channels per octave (in general, we take 12)
** Nbr_Plan = IN: number of scales
**
****************************************************************************
** 
** wave_1d_spline1 (Signal, W_1D, N, Nbr_Plan)
** float ***W_1D, *Signal;
** int N, *Nbr_Plan;
**
** wavelet transform, a trous algorithm,  with a b1-spline 
** Signal = IN: data
** W_1D = OUT: wavelet transform
** N = IN: number of data in the Signal
** Nbr_Plan = OUT: number of scales
**
****************************************************************************
** 
** wave_1d_linear (Signal, W_1D, Np, Nbr_Plan)
** float ***W_1D, *Signal;
** int Np, *Nbr_Plan;
**
** linear wavelet transform, a trous algorithm
**
** Signal = IN: data
** W_1D = OUT: wavelet transform
** N = IN: number of data in the Signal
** Nbr_Plan = OUT: number of scales
**
****************************************************************************
** 
** wave_1d_spline3 (Signal, W_1D, N, Nbr_Plan)
** float ***W_1D, *Signal;
** int N, *Nbr_Plan;
**
** wavelet transform with a b3-spline, a trous algorithm
**
** Signal = IN: data
** W_1D = OUT: wavelet transform
** N = IN: number of data in the Signal
** Nbr_Plan = OUT: number of scales
** 
***************************************************************************
** 
** wave_1d_algo_trou_rec (W_1D, Signal, N, Nbr_Plan)
** int N, *Nbr_Plan;
** float **W_1D, *Signal;
**
** reconstruction from a a-trou algorithm 
**
** Signal = IN: data
** W_1D = OUT: wavelet transform
** N = IN: number of data in the Signal
** Nbr_Plan = OUT: number of scales
** 
***************************************************************************
**
** wave_1d_morlet (Signal,W_1D_re,W_1D_im,N,Nbr_Voie,Nbr_Plan,Nu_0,Scale0)
** float ***W_1D_re, ***W_1D_im, *Signal, Nu_0, *Scale0;
** int N, Nbr_Voie, *Nbr_Plan;
**
**  Morlet transform (Nu_0 > 0.8)
**
** Signal = IN: data
** W_1D_re = OUT: real part of wavelet transform
** W_1D_im = OUT: imaginary part of wavelet transform
** N = IN: number of data in the Signal
** Nbr_Voie = IN: number of channels per octave (in general, we take 12)
** Nbr_Plan = OUT: number of scales
** Nu_0 = IN: Morlet's parameter (Nu_0 must be > 0.8)  
** Scale_0 = OUT: first scale
**
** the wavelet transform W_1D[0..Nbr_Plan-1][0..N] is computed
** W_1D[i] defines the scale i 
** 
****************************************************************************
** 
** wave1d_transform(Signal,Nc,Type_Trans,Nbr_Voie,W_1D,Nbr_Plan,Scale_0,Nu_0)
** float *Signal, ***W_1D, *Scale_0, Nu_0;
** int Nc, Nbr_Voie, Type_Trans, *Nbr_Plan;
**
** 1D wavelet transform
**
** Signal = IN: data
** W_1D = OUT: wavelet transform
** Nc = IN: number of data in the Signal
** Nbr_Voie = IN: number of channels per octave (in general, we take 12)
** Nbr_Plan = OUT: number of scales
** Nu_0 = IN: Morlet's parameter (Nu_0 must be > 0.8)  
** Scale_0 = OUT: first scale
** Type_Trans = type of wavelet transform
**                    TO1_FRENCH 
**                    TO1_MEX
**                    TO1_LINEAR
**                    TO1_B1SPLINE
**                    TO1_B3SPLINE
**                    TO1_MORLET
**                    TO1_ROBUST
**                    TO1_D1GAUS
**
** the wavelet transform W_1D[0..Nbr_Plan-1][0..N] is computed
** W_1D[i] defines the scale i 
** for the Morlet's transform, W_1D[i] define the modulus and
** W_1D[i+Nbr_Plan] defines the phase of the transform
**
****************************************************************************
** 
** wave1d_recons (W_1d,Nc,Nbr_Plan, T_Transf, Nbr_Voie, Signal, Nu_0)
** float **Signal, **W_1D, Nu_0;
** int Nc, Nbr_Voie, T_Transf, Nbr_Plan;
**
** W_1D = IN: wavelet transform
** Signal = OUT: reconstructed signal
** Nc = IN: number of data in the Signal
** Nbr_Voie = IN: number of channels per octave (in general, we take 12)
** Nbr_Plan = IN: number of scales
** Scale_0 = IN: first scale
** T_Transf = IN: type of  transform used
** Nu_0 =  IN:Morlet's parameter 
** 
****************************************************************************/

static char sccsid[] = "@(#)wave1d.c 17.1.1.1 02/01/25 ESO @(#)";

static char Send[100];

#include <stdio.h>
#include <math.h>
#include <string.h>

#include "Def_Math.h"
#include "Def_Mem.h"

#include "Def_Wavelet.h"

/***************************************************************************/

static test_ind (ind, N)
int ind, N;
{
    int Val;
    
    if (ind < 0) Val = 0;
    else
    {
        if (ind >= N) Val = N - 1;
        else Val = ind;
    }
    return (Val);
}

/****************************************************************************/

wave_1d_mex_rec (W_1D, Signal, N, Nbr_Voie, Nbr_Plan)
float **W_1D, *Signal;
int N, Nbr_Voie, Nbr_Plan;
{
    int i,j,k,l,jm,jp;
    float D_Scale, Scale, x, Step;
    float C_psi_mex = PI;
    float tmp, Wave, Scale_0;

    Scale_0 = Scale = 1. / sqrt (3.);
    D_Scale = pow (2., (1. / (float) Nbr_Voie));
    Step = log(D_Scale);

    for (j = 0; j < N; j++) Signal[j] = 0.;

    for (i = 0; i < Nbr_Plan; i++)
    {
        k = 4. * Scale;
        for (j = 0; j < N; j++)
        {
            jm = ( (j-k)  >= 0) ? (j-k) : 0;
            jp = ( (j+k)  < N) ? (j+k) : (N-1);
            tmp = 0.;
            for (l = jm; l < jp; l++)
            {
                x = (float)(j - l)/ Scale;
                x *= x;
                Wave = (1. - x) * exp(-x * .5);
                tmp += Wave * W_1D[i][l];
            }
            Signal[j] += tmp / (Scale * C_psi_mex) * Step;
        }
        Scale *= D_Scale;
    }
}

/***************************************************************************/

wave_1d_mex (Signal, W_1D, N, Nbr_Voie, Nbr_Plan, Scale_0)
float ***W_1D, *Signal, *Scale_0;
int N, Nbr_Voie, *Nbr_Plan;
/* wavelet transform by the mexican hat */
{
    int i,j,k,l,jm,jp;
    float D_Scale, Scale, x, signal;

   /* The minimum scale is such that the step size is just the place
       where the wavelet changes sign.
   */
    *Scale_0 = Scale = 1. / sqrt (3.);

    /*  *Nbr_Plan is the number of scales, i.e. the number of voices per
       octave multiplyed by the number of octaves. The number of octaves
       is the integral part of the binary logarithm of the ratio
       scale_max / scale_min. Since the extension of the
       wavelet is $8 scale$, we have $scale_{max} = n/8$.
    */
    *Nbr_Plan = (float) Nbr_Voie * log( (float) N / (8. * Scale))  / log(2.);
    sprintf (Send,"Nbr_Plan = %d\n", *Nbr_Plan);
    SCTPUT(Send);

    D_Scale = pow (2., (1. / (float) Nbr_Voie));

    /* Allocation memoire pour la transformee */
    *W_1D = f_matrix_alloc (*Nbr_Plan, N);

    for (i = 0; i < *Nbr_Plan; i++)
    {
        k = 4. * Scale;
        for (j = 0; j < N; j++)
        {
            jm = j-k, jp = j + k;
            (*W_1D)[i][j] = 0.;

            for (l = jm; l < jp; l++)
            {
                x = (float)(j - l)/ Scale;
                x *= x;

                if (l <  0)        signal = Signal[-l];      /* Mirroring */
                if (l >= N)        signal = Signal[2*N-2-l]; /* Mirroring */
                if (l>=0 && l<N)   signal = Signal[l];       /* Default value */

                (*W_1D)[i][j] += (1. - x) * exp(-x * .5) * signal;

            }
            (*W_1D)[i][j] /= Scale;
        }
        Scale *= D_Scale;
    }
}

/***************************************************************************/

wave_1d_french (Signal, W_1D, N, Nbr_Voie, Nbr_Plan, Scale_0)
float ***W_1D, *Signal, *Scale_0;
int N, Nbr_Voie, *Nbr_Plan;
/* wavelet transform by the french hat */
{
    int i,j,k,l,jm,jp,jm1,jp1,k3;
    float D_Scale, Scale;

    /* The minimum scale is such that we have just two points on each
    side of the central point for the computation of the integral.
    Remember that the French hat goes from -3 scale to 3 scale (2/3).
    */
    *Scale_0 = Scale = .66;

    /* Since the extension of the
    wavelet is $6 scale$, we have scale_max = n/6.
    */

    /* Nombre d'echelle de la transformee */
    *Nbr_Plan = (float) Nbr_Voie * log((float) N / (6. * Scale)) / log(2.);
    sprintf (Send,"Nbr_Plan = %d\n", *Nbr_Plan);
    SCTPUT(Send);

    D_Scale = pow (2., (1. / (float) Nbr_Voie));
    Scale = .66;

    /* Allocation memoire pour la transformee */
    *W_1D = f_matrix_alloc (*Nbr_Plan, N);

    for (i = 0; i < *Nbr_Plan; i++)
    {
        k = Scale;
        k3 = 3. * Scale;
        for (j = 0; j < N; j++)
        {
            jm = ( (j-k)  >= 0) ? (j-k) : 0;
            jp = ( (j+k)  < N) ? (j+k) : (N-1);
            jm1 = ( (j-k3)  >= 0) ? (j-k3) : 0;
            jp1 = ( (j+k3)  < N) ? (j+k3) : (N-1);
            (*W_1D)[i][j] = 0.;
            for (l = jm1; l < jm; l++) (*W_1D)[i][j] -= Signal[l];
            for (l = jm; l <= jp; l++) (*W_1D)[i][j] += 2*Signal[l];
            for (l = jp+1; l <= jp1; l++) (*W_1D)[i][j] -= Signal[l];
            (*W_1D)[i][j] /= Scale;
        }
        Scale *= D_Scale;
    }
}

/***************************************************************************/

wave_1d_french_rec (W_1D, Signal, N, Nbr_Voie, Nbr_Plan)
float **W_1D, *Signal;
int N, Nbr_Voie, Nbr_Plan;
/* reconstruction */
{
    int i,j,k,l,jm,jp,jm1,jp1,k3;
    float D_Scale, Scale;
    float C_psi_french = 27.;
    float tmp, Step,Scale_0;

    Scale_0 = Scale = .66;
    D_Scale = pow (2., (1. / (float) Nbr_Voie));
    Step = log(D_Scale);

    for (j = 0; j < N; j++) Signal[j] = 0.;

    for (i = 0; i < Nbr_Plan; i++)
    {
        k = Scale;
        k3 = 3. * Scale;

        for (j = 0; j < N; j++)
        {
            jm = ( (j-k)  >= 0) ? (j-k) : 0;
            jp = ( (j+k)  < N) ? (j+k) : (N-1);
            jm1 = ( (j-k3)  >= 0) ? (j-k3) : 0;
            jp1 = ( (j+k3)  < N) ? (j+k3) : (N-1);

            tmp = 0.;
            for (l = jm1; l < jm; l++) tmp -=  W_1D[i][l];
            for (l = jm; l <= jp; l++) tmp += 2. * W_1D[i][l];
            for (l = jp+1; l <= jp1; l++) tmp -= W_1D[i][l];

            Signal[j] += tmp / (Scale * C_psi_french) * Step;
        }
        Scale *= D_Scale;
    }
}


/***************************************************************************/

wave_1d_linear (Signal, W_1D, Np, Nbr_Plan)
float ***W_1D, *Signal;
int Np, *Nbr_Plan;
/* linear wavelet transform */
{
    int i,j,indi1,indi2,Num_Plan;
    int Step;
    float *Data;
    int N = Np;

    /* number of scales of the transform */
    *Nbr_Plan = log((float) N / 4. * 3.) / log(2.);
    sprintf (Send,"Nbr_Plan = %d\n", *Nbr_Plan);
    SCTPUT(Send);

    /* memory allocation */
    *W_1D = f_matrix_alloc (*Nbr_Plan, N);

    Data = f_vector_alloc (N);
    for (i = 0; i < N; i++) Data[i] = Signal[i];
    for (Num_Plan = 0; Num_Plan < *Nbr_Plan - 1; Num_Plan++)
    {
          /* copy the data in the cube */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] = Data [i];

          Step = pow(2., (float) Num_Plan) + 0.5;
          for (i = 0; i < N; i ++)
          {
              indi1 = test_ind (i - Step, N);
              indi2 = test_ind (i + Step, N);
              Data[i] = 0.25 * ((*W_1D)[Num_Plan][indi1] 
                      + (*W_1D)[Num_Plan][indi2]) + 0.5 * (*W_1D)[Num_Plan][i];
          }

          /* calcul the wavelet coefficients */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] -= Data [i];
    }
    
    /* copy the low resolution signal in the cube */
    for (i = 0; i < N; i++) (*W_1D)[*Nbr_Plan-1][i] = Data [i];

    free ((char *) Data);
}

/***************************************************************************/

wave_1d_spline1 (Signal, W_1D, N, Nbr_Plan)
float ***W_1D, *Signal;
int N, *Nbr_Plan;
/* wavelet transform with a b1-spline */
{
    int i,indi1,indi2,Num_Plan;
    int Step;
    float *Data;

    /* The size of the wavelet for the i th scale is 2^i. We
       want to have a correct computation of the wavelet coefficients
       for at least a quater of the signal at the largest scale.
       Therefore we must have a number \verb|nx| of scales such that
       2^{Nbr_Plan} = 3 n / 4
    */

    /* number of scales of the transform */
    *Nbr_Plan = log(3. * (float) N / 4.) / log(2.);
    sprintf (Send,"Nbr_Plan = %d\n", *Nbr_Plan);
    SCTPUT(Send);

    /* memory allocation */
    *W_1D = f_matrix_alloc (*Nbr_Plan, N);

    Data = f_vector_alloc (N);
    for (i = 0; i < N; i++) Data[i] = Signal[i];

    for (Num_Plan = 0; Num_Plan < *Nbr_Plan - 1; Num_Plan++)
    {

          /* copy the image in the cube */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] = Data [i];

          Step = pow(2., (float) Num_Plan) + 0.5;

          for (i = 0; i < N; i ++)
          {
              indi1 = test_ind (i - Step, N);
              indi2 = test_ind (i + Step, N);
              Data[i] =(0.5*((*W_1D)[Num_Plan][indi1]+(*W_1D)[Num_Plan][indi2])
                          + 2. * (*W_1D)[Num_Plan][i]) / 3.;
          }

          /* Calcul des coefficients d'ondelettes */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] -= Data [i];
    }
    
    /* copy the law resolution signal in the cube */
    for (i = 0; i < N; i++) (*W_1D)[*Nbr_Plan-1][i] = Data [i];

    free ((char *) Data);
}

/***************************************************************************/

wave_1d_spline3 (Signal, W_1D, N, Nbr_Plan)
float ***W_1D, *Signal;
int N, *Nbr_Plan;
/* wavelet transform with a b3-spline */
{
    int i,indi1,indi2,indi3,indi4,Num_Plan;
    int Step;
    float *Data;

    /* number of scales of the transform */
    /* The size of the wavelet for the i th scale is 2^{i+1}.
       We want to have a correct computation of the wavelet coefficients
       for at least a quater of the signal at the largest scale.
       Therefore we must have a number \verb|nx| of scales such that
       2^{nx+1} = 3 n / 4
    */
    *Nbr_Plan = log(3. * (float) N / 4.) / log(2.);
    sprintf (Send,"Nbr_Plan = %d\n", *Nbr_Plan);
    SCTPUT(Send);
    wave_1d_trou(Signal, W_1D, N, *Nbr_Plan);
}

/***************************************************************************/

wave_1d_robust (Signal, W_1D, N, Nbr_Plan)
float ***W_1D, *Signal;
int N, *Nbr_Plan;
/* wavelet transform with a b3-spline */
{
    int i,indi1,indi2,indi3,indi4,Num_Plan;
    int Step;
    float *Data;

    /* number of scales of the transform */
    /* The size of the wavelet for the i th scale is 2^{i+1}.
       We want to have a correct computation of the wavelet coefficients
       for at least a quater of the signal at the largest scale.
       Therefore we must have a number \verb|nx| of scales such that
       2^{nx+1} = 3 n / 4
    */
    *Nbr_Plan = log(3. * (float) N / 4.) / log(2.);
    sprintf (Send,"Nbr_Plan = %d\n", *Nbr_Plan);
    SCTPUT(Send);
    wave_1d_trou_med (Signal, W_1D, N, *Nbr_Plan);
}

/***************************************************************************/

wave_1d_trou (Signal, W_1D, N, Nbr_Plan)
float ***W_1D, *Signal;
int N, Nbr_Plan;
/* wavelet transform with a b3-spline */
{
    int i,indi1,indi2,indi3,indi4,Num_Plan;
    int Step;
    float *Data;

    /* memory allocation of the transform */
    *W_1D = f_matrix_alloc (Nbr_Plan, N);

    Data   = f_vector_alloc (N);
    for (i = 0; i < N; i++) Data[i] = Signal[i];

    for (Num_Plan = 0; Num_Plan < Nbr_Plan - 1; Num_Plan++)
    {
          /* copy the image in a cube */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] = Data [i];

          Step = pow(2., (float) Num_Plan) + 0.5;

          for (i = 0; i < N; i ++)
          {
              indi1 = test_ind (i - Step, N);
              indi2 = test_ind (i + Step, N);
              indi3 = test_ind (i - 2 * Step, N);
              indi4 = test_ind (i + 2 * Step, N);

              Data[i] =
                  0.0625 * ( (*W_1D)[Num_Plan][indi3]+(*W_1D)[Num_Plan][indi4])
                  + 0.25 * ( (*W_1D)[Num_Plan][indi1]+(*W_1D)[Num_Plan][indi2])
                          + 0.375 * (*W_1D)[Num_Plan][i];
          }

          /* Calcule the wavelet coefficients */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] -= Data [i];
    }
    
    /* copy the low resolution */
    for (i = 0; i < N; i++) (*W_1D)[Nbr_Plan-1][i] = Data [i];

    free ((char *) Data);
}

/***************************************************************************/

wave_1d_algo_trou_rec (W_1D, Signal, N, Nbr_Plan)
/* reconstruction from a a-trou algorithm */
float **W_1D, *Signal;
int N, Nbr_Plan;
{
    int i,j;

    for (i = 0; i < N; i++)
    {
        Signal[i] = 0.;
        for (j = 0; j < Nbr_Plan; j++) Signal[i] += W_1D[j][i];
    }
}

/***************************************************************************/

wave_1d_morlet (Signal, W_1D_re, W_1D_im, N, Nbr_Voie, Nbr_Plan, Nu_0, Scale0)
float ***W_1D_re, ***W_1D_im, *Signal, Nu_0, *Scale0;
int N, Nbr_Voie, *Nbr_Plan;
/* Morlet transform (Nu_0 > 0.8) */
{
    int i,j,k,l,jm,jp;
    float D_Scale, x;
    float Norm, Omega, Val, Coef;
    float Scale;

    *Scale0 = Scale = 2. * Nu_0;
    /* Number of scales */
    *Nbr_Plan = (float) Nbr_Voie * log((float) N / (12. * Scale)) / log(2.);

    D_Scale = pow (2., (1. / (float) Nbr_Voie));
    sprintf (Send,"Nbr_Plan = %d, Scale = %f\n", *Nbr_Plan, Scale);
    SCTPUT(Send);

    Norm = 1. / sqrt(2.*PI);
    Omega = 2. * PI * Nu_0;

    *W_1D_re = f_matrix_alloc (*Nbr_Plan, N);
    *W_1D_im = f_matrix_alloc (*Nbr_Plan, N);

    for (i = 0; i < *Nbr_Plan; i++)
    {
        k = 6. * Scale;
        for (j = 0; j < N; j++)
        {
            jm = ( (j-k)  >= 0) ? (j-k) : 0;
            jp = ( (j+k)  < N) ? (j+k) : (N-1);
            (*W_1D_re)[i][j] = 0.;
            (*W_1D_im)[i][j] = 0.;
            for (l = jm; l < jp; l++)
            {
                x = (float)(j - l)/ Scale;
                Coef = Norm * exp(-x*x/2);
                Val = Omega * x;
                (*W_1D_re)[i][j] += Coef * cos(Val) * Signal[l];
                (*W_1D_im)[i][j] -= Coef * sin(Val) * Signal[l];
            }
            (*W_1D_re)[i][j] /= Scale;
            (*W_1D_im)[i][j] /= Scale;
        }
        Scale *= D_Scale;
    }
}

/***************************************************************************/

static wave_mod (Wave_re, Wave_im, Wave_Mod, N, Nbr_Plan)
float **Wave_re, **Wave_im, **Wave_Mod;
int N, Nbr_Plan;
{
    int i,j;

    for (i = 0; i < Nbr_Plan; i++)
    for (j = 0; j < N; j++)
            Wave_Mod[i][j] = sqrt(Wave_re[i][j]*Wave_re[i][j] 
                                    +Wave_im[i][j]*Wave_im[i][j]);
}

/***************************************************************************/

static wave_phase (Wave_re, Wave_im, Wave_Phase, N, Nbr_Plan)
float **Wave_re, **Wave_im, **Wave_Phase;
int N, Nbr_Plan;
{
    int i,j;

    for (i = 0; i < Nbr_Plan; i++)
    {
        for (j = 0; j < N; j++)
        {
            ARG(Wave_re[i][j],Wave_im[i][j], Wave_Phase[Nbr_Plan+i][j]);
        }
    }
}

/***************************************************************************/

static wave_re (Wave_Mod, Wave_Ph, Wave_re, N, Nbr_Plan)
float **Wave_Mod, **Wave_Ph, **Wave_re;
int N, Nbr_Plan;
{
    int i,j;

    for (i = 0; i < Nbr_Plan; i++)
    for (j = 0; j < N; j++)
            Wave_re[i][j] = Wave_Mod[i][j] * cos(Wave_Ph[i][j]);
}

/***************************************************************************/

static wave_im (Wave_Mod, Wave_Ph, Wave_im, N, Nbr_Plan)
float **Wave_Mod, **Wave_Ph, **Wave_im;
int N, Nbr_Plan;
{
    int i,j;

    for (i = 0; i < Nbr_Plan; i++)
    for (j = 0; j < N; j++)
            Wave_im[i][j] = Wave_Mod[i][j] * sin(Wave_Ph[i][j]);
}

/***************************************************************************/

wave1d_transform (Signal, Nc, Type_Transform, Nbr_Voie, W_1D,
                 Nbr_Plan, Scale_0, Nu_0)
float *Signal, ***W_1D, *Scale_0, Nu_0;
int Nc, Nbr_Voie, Type_Transform, *Nbr_Plan;
{
    float **W_1D_re, **W_1D_im;
    float **Ptr;

    *Scale_0 = 0.;

    switch (Type_Transform)
    {
       case TO1_FRENCH: 
         wave_1d_french (Signal, W_1D, Nc, Nbr_Voie, Nbr_Plan, Scale_0);
         break;
       case TO1_MEX: 
         wave_1d_mex (Signal, W_1D, Nc, Nbr_Voie, Nbr_Plan, Scale_0);
         break;
       case TO1_LINEAR:
         wave_1d_linear (Signal, W_1D, Nc, Nbr_Plan); 
         break;
       case TO1_B1SPLINE :
         wave_1d_spline1 (Signal, W_1D, Nc, Nbr_Plan);
         break;
       case TO1_B3SPLINE :
         wave_1d_spline3 (Signal, W_1D, Nc, Nbr_Plan);
         break;
       case TO1_ROBUST :
         wave_1d_robust (Signal, W_1D, Nc, Nbr_Plan);
         break;
       case TO1_D1GAUS :
         wave_1d_d1gaus (Signal, W_1D, Nc, Nbr_Voie, Nbr_Plan, Scale_0);
         break;
       case TO1_MORLET:
         wave_1d_morlet (Signal, &W_1D_re, &W_1D_im, Nc, 
                         Nbr_Voie, Nbr_Plan, Nu_0, Scale_0);
         *W_1D = f_matrix_alloc (2*(*Nbr_Plan), Nc);
         wave_mod (W_1D_re, W_1D_im, *W_1D, Nc, *Nbr_Plan);
         wave_phase (W_1D_re, W_1D_im, *W_1D, Nc, *Nbr_Plan);
         f_matrix_free(W_1D_re,*Nbr_Plan);
         f_matrix_free(W_1D_im,*Nbr_Plan);
         break;
       }
} 

/*************************************************************************/

wave1d_recons (W_1D, Nc, Nbr_Plan, T_Transf, Nbr_Voie, Signal, Nu_0)
float **W_1D, **Signal, Nu_0;
int Nc, Nbr_Voie, T_Transf, Nbr_Plan;
{
    float **W_1D_re, **W_1D_im;
    float **Ptr;

    *Signal = f_vector_alloc (Nc);
    switch (T_Transf)
    {
       case TO1_FRENCH: 
         wave_1d_french_rec (W_1D, *Signal, Nc, Nbr_Voie, Nbr_Plan); 
         break;
       case TO1_MEX: 
         wave_1d_mex_rec (W_1D, *Signal, Nc, Nbr_Voie, Nbr_Plan);
         break;
       case TO1_LINEAR:
       case TO1_B1SPLINE :
       case TO1_B3SPLINE :
       case TO1_ROBUST :
         wave_1d_algo_trou_rec (W_1D, *Signal, Nc, Nbr_Plan);
         break;
       case TO1_D1GAUS :
         wave_1d_d1gaus_rec (Signal, W_1D, Nc, Nbr_Voie, Nbr_Plan);
         break;
       case TO1_MORLET:
/*
         W_1D_re = f_matrix_alloc (Nbr_Plan, Nc);
         W_1D_im = f_matrix_alloc (Nbr_Plan, Nc);
         Ptr = W_1D + Nbr_Plan * Nc;
         wave_re (W_1D, Ptr, W_1D_re, Nc, Nbr_Plan);
         wave_im (W_1D, Ptr, W_1D_im, Nc, Nbr_Plan);
         sprintf (Send,"Morlet's reconstruction not implemented \n");
         SCTPUT(Send);
         f_matrix_free(W_1D_re,Nbr_Plan);
         f_matrix_free(W_1D_im,Nbr_Plan);*/
         break;
       }
}

/*************************************************************************/
/*  Sorting algorithm:      Heapsort method                      */
/*  From: Numerical Recipes, Cambridge University Press, p. 226  */

static int sort(n,ra)

int n;
float ra[];
{
        int l,j,ir,i;
        float rra;

        l=(n >> 1)+1;
        ir=n;
        for (;;) {
                if (l > 1)
                        rra=ra[--l];
                else {
                        rra=ra[ir];
                        ra[ir]=ra[1];
                        if (--ir == 1) {
                                ra[1]=rra;
                                return (0);
                        }
                }
                i=l;
                j=l << 1;
                while (j <= ir) {
                        if (j < ir && ra[j] < ra[j+1]) ++j;
                        if (rra < ra[j]) {
                                ra[i]=ra[j];
                                j += (i=j);
                        }
                        else j=ir+1;
                }
                ra[i]=rra;
        }
}

/***********************************************************************/
static float median(List, N)

float *List;
int   N;

{
float *Histo, med;
int   pos, i;

pos = (N+1)/2 - 1;
Histo = f_vector_alloc(N);
for (i=0; i<N; i++) Histo[i] = List[i];
sort(N, &Histo[-1]);
med = Histo[pos];
free((char *) Histo);
return(med);

}

/***********************************************************************/
static int mirror(pos,N)

int pos,N;

{
if (pos < 0) return(-pos);
if (pos >= N) return(2*N-pos-2);
return(pos);

}

/***************************************************************************/

wave_1d_trou_med (Signal, W_1D, N, Nbr_Plan)
float ***W_1D, *Signal;
int N, Nbr_Plan;
/* wavelet transform with a b3-spline */
{
    int i,indi1,indi2,indi3,indi4,Num_Plan;
    int Step;
    float *Data;

    /* memory allocation of the transform */
    *W_1D = f_matrix_alloc (Nbr_Plan, N);

    Data   = f_vector_alloc (N);
    for (i = 0; i < N; i++) Data[i] = Signal[i];

    for (Num_Plan = 0; Num_Plan < Nbr_Plan - 1; Num_Plan++)
    {
          /* copy the image in a cube */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] = Data [i];

          Step = pow(2., (float) Num_Plan) + 0.5;

          for (i = 0; i < N; i ++)
          {
              indi1 = mirror (i - Step, N);
              indi2 = mirror (i + Step, N);
              indi3 = mirror (i - 2 * Step, N);
              indi4 = mirror (i + 2 * Step, N);

              Data[i] =
                  0.0625 * ( (*W_1D)[Num_Plan][indi3]+(*W_1D)[Num_Plan][indi4])
                  + 0.25 * ( (*W_1D)[Num_Plan][indi1]+(*W_1D)[Num_Plan][indi2])
                          + 0.375 * (*W_1D)[Num_Plan][i];
          }

          /* Calcule the wavelet coefficients */
          for (i = 0; i < N; i++) (*W_1D)[Num_Plan][i] -= Data [i];
    }
    
    /* copy the low resolution */
    for (i = 0; i < N; i++) (*W_1D)[Nbr_Plan-1][i] = Data [i];

    free ((char *) Data);
}

/***************************************************************************/

wave_1d_d1gaus (Signal, W_1D, N, Nbr_Voie, Nbr_Plan, Scale_0)
float ***W_1D, *Signal, *Scale_0;
int N, Nbr_Voie, *Nbr_Plan;
/* wavelet transform by the mexican hat */
{
    int i,j,k,l,jm,jp;
    float D_Scale, Scale, x;

   /* The minimum scale is such that the step size is just the place
       where the wavelet changes sign.
   */
    *Scale_0 = Scale = 1. / sqrt (3.);

    /*  *Nbr_Plan is the number of scales, i.e. the number of voices per
       octave multiplyed by the number of octaves. The number of octaves
       is the integral part of the binary logarithm of the ratio
       scale_max / scale_min. Since the extension of the
       wavelet is $8 scale$, we have $scale_{max} = n/8$.
    */
    *Nbr_Plan = (float) Nbr_Voie * log( (float) N / (8. * Scale))  / log(2.);
    sprintf (Send,"Nbr_Plan = %d\n", *Nbr_Plan);
    SCTPUT(Send);

    D_Scale = pow (2., (1. / (float) Nbr_Voie));

    /* Allocation memoire pour la transformee */
    *W_1D = f_matrix_alloc (*Nbr_Plan, N);

    for (i = 0; i < *Nbr_Plan; i++)
    {
        k = 4. * Scale;
        for (j = 0; j < N; j++)
        {
            jm = ( (j-k)  >= 0) ? (j-k) : 0;
            jp = ( (j+k)  < N) ? (j+k) : (N-1);
            (*W_1D)[i][j] = 0.;
            for (l = jm; l < jp; l++)
            {
                x = (float)(j - l)/ Scale;
                (*W_1D)[i][j] += (-1.)*x * exp(-x*x * .5) * Signal[l];
            }
            (*W_1D)[i][j] /= Scale;
        }
        Scale *= D_Scale;
    }
}

/****************************************************************************/

wave_1d_d1gaus_rec (W_1D, Signal, N, Nbr_Voie, Nbr_Plan)

float **W_1D, *Signal;
int N, Nbr_Voie, Nbr_Plan;
{
    int i,j,k,l,jm,jp;
    float D_Scale, Scale, x, Step;
    float C_psi_mex = PI;
    float tmp, Wave, Scale_0;

    Scale_0 = Scale = 1. / sqrt (3.);
    D_Scale = pow (2., (1. / (float) Nbr_Voie));
    Step = log(D_Scale);

    for (j = 0; j < N; j++) Signal[j] = 0.;

    for (i = 0; i < Nbr_Plan; i++)
    {
        k = 4. * Scale;
        for (j = 0; j < N; j++)
        {
            jm = ( (j-k)  >= 0) ? (j-k) : 0;
            jp = ( (j+k)  < N) ? (j+k) : (N-1);
            tmp = 0.;
            for (l = jm; l < jp; l++)
            {
                x = (float)(j - l)/ Scale;
                Wave = (-1.)*x * exp(-x*x * .5);
                tmp += Wave * W_1D[i][l];
            }
            Signal[j] += tmp / (Scale * C_psi_mex) * Step;
        }
        Scale *= D_Scale;
    }
}

/***************************************************************************/