/************************************************************************
 * The following is a collection of wavelet-analysis related routines
 * @ MAR-1994 EK (original code in fortran)
 * @ JUN-2001 - ported to C
 * void wavelet(double *array,
                int     size,
                int     filter,
                int     direction)
 * void wavelet_haar(double *array,
                     int     size,
                     int     direction)
 * void wavelet_haar_dec(double *array,
                         int     size)
 * int  wavelet_haar_dec_1level(double *array,
                                int     size)
 * void wavelet_haar_rec(double *array,
                         int     size)
 * int  wavelet_haar_rec_1level(double *array,
                                int     size)
 * void wavelet_D4(double *array,
                   int     size,
                   int     direction)
 * void wavelet_D4_dec(double *array,
                       int     size)
 * int  wavelet_D4_dec_1level(double *array,
                              int     size)
 * void wavelet_D4_rec(double *array,
                       int     size)
 * int  wavelet_D4_rec_1level(double *array,
                              int     size)
 * void wavelet_D8(double *array,
                   int     size,
                   int     direction)
 * void wavelet_D8_dec(double *array,
                       int     size)
 * int  wavelet_D8_dec_1level(double *array,
                              int     size)
 * void wavelet_D8_rec(double *array,
                       int     size)
 * int  wavelet_D8_rec_1level(double *array,
                              int     size)
 * void wavelet_D16(double *array,
                    int     size,
                    int     direction)
 * void wavelet_D16_dec(double *array,
                        int     size)
 * int  wavelet_D16_dec_1level(double *array,
                               int     size)
 * void wavelet_D16_rec(double *array,
                        int     size)
 * int  wavelete_D16_rec_1level(double *array,
                                int     size)
 ************************************************************************/

/************************************************************************
 * The following is a collection of wavelet-analysis related routines
 * @ MAR-1994 EK
 * @ JUL-2001 - ported to C
 ************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "mywv2lib.h"

/************************************************************************
 * Top-level call:
 * array : Holds the (input) data, and upon execution the wavelet (output) data
 * size  : size of the array
 * filter: Controls which filter to use
 *         1 -> Haar (or D-2) filter
 *         2 -> Daubechies  4-tap compact support wavelet filter
 *         3 -> Daubechies  8-tap compact support wavelet filter
 *         4 -> Daubechies 16-tap compact support wavelet filter
 *         Anything else -> Does nothing
 *         Rule of thumb : 2^filter is the filter support
 * direction: 1 -> Forward (or decomposition) transform
 *           -1 -> Inverse (or reconstruction) transform
 *           Anything else -> Does nothing
 ************************************************************************/
void wavelet(double *array,
             int     size,
             int     filter,
             int     direction)
{
   switch(filter)
   {
   case 1:
      wavelet_haar(array,size,direction);
      break;
   case 2:
      wavelet_D4(array,size,direction);
      break;
   case 3:
      wavelet_D8(array,size,direction);
      break;
   case 4:
      wavelet_D16(array,size,direction);
      break;
   default:
      // Do nothing
      break;
   }
}

/************************************************************************
 * Implements the case filter = 1 discussed above
 ************************************************************************/
void wavelet_haar(double *array,
                  int     size,
                  int     direction)
{
   // 2 cases: Forward (direction=1) or Inverse (direction=-1) transform
   if(direction==1)
   {
      wavelet_haar_dec(array,size);
   }
   else if(direction==(-1))
   {
      wavelet_haar_rec(array,size);
   }
   else
   {
      // Do nothing
   }
}

/************************************************************************
 * Implements the forward (or decomposition) Haar transform
 * discussed above.
 ************************************************************************/
void wavelet_haar_dec(double *array,
                      int     size)
{
   int lsize;

   lsize = size;
   while(lsize>1)
   {
      // The idea is that lsize returns the size of the low-pass band
      // For example, if size=512, then after 1-level of decomposition
      // (is the name of the function more clear now ?) lsize = 256
      // and so on.
      lsize = wavelet_haar_dec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of decomposition, used when doing forward Haar transform
 * Haar filter: [ 1  1 ] * 1/sqrt(2.0)
 *              [ 1 -1 ]
 ************************************************************************/
int  wavelet_haar_dec_1level(double *array,
                             int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int      size_low;
   int      size_high;
   double   sum;
   double   diff;
   double   temp1;
   double   temp2;
   double  *larray;
   double   sqrt2   = 1.41421356237309504880168872420969808;
   double   sqrt1_2 = 0.707106781186547524400844362104849039;

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      // sqrt1_2 = sqrt2/2.0;
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
      k = 0;
      l = 1;
      for(i=0,j=size_low;i<size_high;i++,j++)
      {
         // For example: i=0, k=0, l=1
         //              i=1, k=2, l=3
         //              i=2, k=4, l=5 ...
         temp1 = array[k];
         temp2 = array[l];
         sum   = temp1 + temp2;
         diff  = temp1 - temp2;
         larray[i] = sqrt1_2*sum;
         larray[j] = sqrt1_2*diff;
         k += 2;
         l += 2;
      }
      // The following case should only happen when size is odd
      // In this case, we want to keep the "extra" sample on the low-pass band
      // Assume that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      if(size_low!=size_high)
      {
         larray[size_low-1] = array[size-1]*sqrt2;
      }
      memcpy(array,larray,size*sizeof(double));
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}

/************************************************************************
 * Implements the inverse (or reconstruction) Haar transform
 * discussed above.
 ************************************************************************/
void wavelet_haar_rec(double *array,
                      int     size)
{
   register int i;
   register int j;
   int          lsize;
   int          count;

   // We first need to know how many levels of decomposition were performed
   // during the forward (decomposition) transform
   //   size<= 1 -> 0 levels
   // 1<size<= 2 -> 1 level
   // 2<size<= 4 -> 2 levels
   // 4<size<= 8 -> 3 levels
   // 8<size<=16 -> 4 levels ...
   count = 0;
   lsize = size;
   while(lsize>1)
   {
      lsize = lsize - lsize/2;
      count++;
   }
   // At this point, count = # of decomposition levels

   // For a given decomposition level, need to find out what is the
   // appropriate size to use.
   for(i=count-1;i>=0;i--)
   {
      lsize = size;
      for(j=0;j<i;j++)
      {
         lsize = lsize - lsize/2;
      }
      lsize = wavelet_haar_rec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of reconstruction, used when doing inverse Haar transform
 * Haar filter: [ 1  1 ] * 1/sqrt(2.0)
 *              [ 1 -1 ]
 ************************************************************************/
int  wavelet_haar_rec_1level(double *array,
                             int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int          size_low;
   int          size_high;
   double       sum;
   double       diff;
   double       temp1;
   double       temp2;
   double      *larray;
   double       sqrt2   = 1.41421356237309504880168872420969808;
   double       sqrt1_2 = 0.707106781186547524400844362104849039;

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      // sqrt1_2 = sqrt2/2.0;
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
      k = 0;
      l = 1;
      for(i=0,j=size_low;i<size_high;i++,j++)
      {
         // For example: i=0, k=0, l=1
         //              i=1, k=2, l=3
         //              i=2, k=4, l=5 ...
         temp1 = array[i];
         temp2 = array[j];
         sum   = temp1 + temp2;
         diff  = temp1 - temp2;
         larray[k] = sqrt1_2*sum;
         larray[l] = sqrt1_2*diff;
         k += 2;
         l += 2;
      }
      // The following case should only happen when size is odd
      // In this case, we have kept the "extra" sample on the low-pass band
      // assuming that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      // So, sample = sum/SQRT2
      if(size_low!=size_high)
      {
         larray[size-1] = sqrt1_2*array[size_low-1];
      }
      memcpy(array,larray,size*sizeof(double));
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}

/************************************************************************
 * Implements the case filter = 2 discussed above,
 * i.e. Daubechies-4 compact support tap-4 wavelet transform
 ************************************************************************/
void wavelet_D4(double *array,
                int     size,
                int     direction)
{
   // 2 cases: Forward (direction=1) or Inverse (direction=-1) transform
   if(direction==1)
   {
      wavelet_D4_dec(array,size);
   }
   else if(direction==(-1))
   {
      wavelet_D4_rec(array,size);
   }
   else
   {
      // Do nothing
   }
}

/************************************************************************
 * Implements the forward (or decomposition) Daubechies-4 compact support
 * 4-tap wavelet transform discussed above.
 ************************************************************************/
void wavelet_D4_dec(double *array,
                    int     size)
{
   int lsize;

   lsize = size;
   while(lsize>1)
   {
      // The idea is that lsize returns the size of the low-pass band
      // For example, if size=512, then after 1-level of decomposition
      // (is the name of the function more clear now ?) lsize = 256
      // and so on.
      lsize = wavelet_D4_dec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of decomposition, used when doing Daubechies-4
 * compact support 4-tap wavelet transform
 * Daubechies-4 filter:   
 * [ 0.482962913145  0.836516303738  0.224143868042 -0.129409522551 ] -LOW
 * [-0.129409522551 -0.224143868042  0.836516303738 -0.482962913145 ] -HIGH
 * "Allignment" : 
 *    x(2)   x(3)   x(4)    x(5)  x(6)   x(7)
 *               [ 0.483  0.837  0.224 -0.129 ]         -LOW
 * [-0.129 -0.224  0.837 -0.483 ]                       -HIGH
 ************************************************************************/
int wavelet_D4_dec_1level(double *array,
                          int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int          size_low;
   int          size_high;
   double       temp;
   double       sum;
   double       diff;
   double       sign;
   double      *larray;
   double       sqrt2   = 1.41421356237309504880168872420969808;
   double       sqrt1_2 = 0.707106781186547524400844362104849039;
   double       D4[4] = { 0.482962913145,
                          0.836516303738,
                          0.224143868042,
                         -0.129409522551 };

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
	  // adjust size to be even
	  size = 2*size_high;
      k = 0;
      l = size+1;
      for(i=0;i<size-1;i+=2)
      {
         sum  =  0.0;
         diff =  0.0;
         sign = -1.0;
         for(j=0;j<4;j++)
         {
            temp  = D4[j];
            // size=1024, i=0, k= 0, l=1025, j=0  ==>  sum=0  diff=1
            //                               j=1  ==>  sum=1  diff=0
            //                               j=2  ==>  sum=2  diff=1023
            //                               j=3  ==>  sum=3  diff=1022
            //            i=2, k= 2, l=1027, j=0  ==>  sum=2  diff=3
            //                               j=1  ==>  sum=3  diff=2
            //                               j=2  ==>  sum=4  diff=1
            //                               j=3  ==>  sum=5  diff=0
            sum  += temp*array[(k+j)%size];
            diff += sign*temp*array[(l-j)%size];
            sign  = -sign;
		 }
         larray[i]   = sum;
         larray[i+1] = diff;
         k += 2;
         l += 2;
      }
      // The following case should only happen when size is odd
      // In this case, we want to keep the "extra" sample on the low-pass band
      // Assume that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      if(size_low!=size_high)
      {
		 size++;
         array[size_low-1] = array[size-1]*sqrt2;
      }
      for(i=0;i<size-1;i+=2)
      {
         j = i/2;
         array[j]          = larray[i];
         array[j+size_low] = larray[i+1];
      }
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}

/************************************************************************
 * Implements the inverse (or reconstruction) Daubechies-4
 * compact support 4-tap wavelet transform discussed above.
 ************************************************************************/
void wavelet_D4_rec(double *array,
                    int     size)
{
   register int i;
   register int j;
   int          lsize;
   int          count;

   // We first need to know how many levels of decomposition were performed
   // during the forward (decomposition) transform
   //   size<= 1 -> 0 levels
   // 1<size<= 2 -> 1 level
   // 2<size<= 4 -> 2 levels
   // 4<size<= 8 -> 3 levels
   // 8<size<=16 -> 4 levels ...
   count = 0;
   lsize = size;
   while(lsize>1)
   {
      lsize = lsize - lsize/2;
      count++;
   }
   // At this point, count = # of decomposition levels

   // For a given decomposition level, need to find out what is the
   // appropriate size to use.
   for(i=count-1;i>=0;i--)
   {
      lsize = size;
      for(j=0;j<i;j++)
      {
         lsize = lsize - lsize/2;
      }
      lsize = wavelet_D4_rec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of reconstruction, used when doing inverse 
 * Daubechies-4 compact support 4-tap wavelet transform
 * Daubechies-4 filter:   
 * [ 0.482962913145  0.836516303738  0.224143868042 -0.129409522551 ] -LOW
 * [-0.129409522551 -0.224143868042  0.836516303738 -0.482962913145 ] -HIGH
 * "Allignment" : 
 *     c[512]  d[1]   c[1]   d[2]   c[2]   d[3]   c[3]    d[4]   c[4]
 *     Notice: "highs" go first and then "lows" i.e.
 *     d(512) c(512) | d(1) c(1) d(2) d(3) .... d(512) c(512) | d(1) c(1)
 *     d4[3]   d4[2]  d4[1]  d4[4]
 *     d4[4]  -d4[1]  d4[2] -d4[3]
 *                      3      2      1      4          -> low
 *                      4     -1      2     -3          -> high
 ************************************************************************/
int wavelet_D4_rec_1level(double *array,
                          int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int          size_low;
   int          size_high;
   double       temp1;
   double       temp2;
   double       temp3;
   double       temp4;
   double       sum;
   double       diff;
   double      *larray;
   double       sqrt2   = 1.41421356237309504880168872420969808;
   double       sqrt1_2 = 0.707106781186547524400844362104849039;
   double       D4[4] = { 0.482962913145,
                          0.836516303738,
                          0.224143868042,
                         -0.129409522551 };

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
	  // adjust size to be even
	  size = 2*size_high;
      for(i=0,k=size_high;i<size-1;i+=2,k++)
      {
         sum  = 0.0;
         diff = 0.0;
         for(j=0,l=0;j<4;j+=2,l++)
         {
            temp1 = D4[j];
            temp2 = D4[j+1];
            temp3 = array[(k-l)%size_high];
            temp4 = array[(k+l)%size_high+size_low];
            // For example : i=0, j=0, k=512, l=0, temp3=  0, temp4=512
            // (size=1024)   i=0, j=2, k=512, l=1, temp3=511, temp4=513
            //               i=2, j=0, k=513, l=0, temp3=  2, temp4=513
            //               i=2, j=2, k=513, l=1, temp3=  1, temp4=514
            //               i=4, j=0, k=514, l=0, temp3=  3, temp4=514
            //               i=4, j=2, k=514, l=1, temp3=  2, temp4=515
            //               i=6, j=0, k=515, l=0, temp3=  4, temp4=515
            //               i=6, j=2, k=515, l=1, temp3=  3, temp4=516
            sum  += temp1*temp3 + temp2*temp4;
            diff += temp2*temp3 - temp1*temp4;
		 }
         larray[i]   = sum;
         larray[i+1] = diff;
      }
      // The following case should only happen when size is odd
      // In this case, we have kept the "extra" sample on the low-pass band
      // assuming that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      // So, sample = sum/SQRT2
      if(size_low!=size_high)
      {
         size++;
         larray[i]   = array[size_low-1]*sqrt1_2;
      }
      memcpy(array,larray,size*sizeof(double));
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}

/************************************************************************
 * Implements the case filter = 3 discussed above,
 * i.e. Daubechies-8 compact support tap-8 wavelet transform
 ************************************************************************/
void wavelet_D8(double *array,
                int     size,
                int     direction)
{
   // 2 cases: Forward (direction=1) or Inverse (direction=-1) transform
   if(direction==1)
   {
      wavelet_D8_dec(array,size);
   }
   else if(direction==(-1))
   {
      wavelet_D8_rec(array,size);
   }
   else
   {
      // Do nothing
   }
}

/************************************************************************
 * Implements the forward (or decomposition) Daubechies-8 compact support
 * 8-tap wavelet transform discussed above.
 ************************************************************************/
void wavelet_D8_dec(double *array,
                    int     size)
{
   int lsize;

   lsize = size;
   while(lsize>1)
   {
      // The idea is that lsize returns the size of the low-pass band
      // For example, if size=512, then after 1-level of decomposition
      // (is the name of the function more clear now ?) lsize = 256
      // and so on.
      lsize = wavelet_D8_dec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of decomposition, used when doing Daubechies-8
 * compact support 8-tap wavelet transform
 * Daubechies-8 filter:   
 * [  0.230377813309, 0.714846570553, 0.630880767930, -0.027983769417,
 *   -0.187034811719, 0.030841381836, 0.032883011667, -0.010597401785 ] - LOW
 * [ -0.010597401785,-0.032883011667, 0.030841381836,  0.187034811719,
 *   -0.027983769417,-0.630880767930, 0.714846570553, -0.230377813309 ] - HIGH
 * "Allignment" : 
 *    0  1  2  3  4  5  6  7  8  9 10 11 12 13
 *                      .  .  .  .  .  .  .  .   -LOW
 *    .  .  .  .  .  .  .  .                     -HIGH
 ************************************************************************/
int wavelet_D8_dec_1level(double *array,
                          int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int          size_low;
   int          size_high;
   double       temp;
   double       sum;
   double       diff;
   double       sign;
   double      *larray;
   double       sqrt2   = 1.41421356237309504880168872420969808;
   double       sqrt1_2 = 0.707106781186547524400844362104849039;
   double       D8[8] = { 0.230377813309,
                          0.714846570553,
                          0.630880767930,
                         -0.027983769417,
                         -0.187034811719,
                          0.030841381836,
                          0.032883011667,
                         -0.010597401785 };

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
	  // adjust size to be even
	  size = 2*size_high;
      k = 0;
      l = 3*size+1;
      for(i=0;i<size-1;i+=2)
      {
         sum  =  0.0;
         diff =  0.0;
         sign = -1.0;
         for(j=0;j<8;j++)
		 {
            temp  = D8[j];
            // size=1024, i=0, k= 0, l=3075, j=0  ==>  sum=0  diff=1
            //                               j=1  ==>  sum=1  diff=0
            //                               j=2  ==>  sum=2  diff=1023
            //                               j=3  ==>  sum=3  diff=1022
            //                               j=4  ==>  sum=4  diff=1021
            //                               j=5  ==>  sum=5  diff=1020
            //                               j=6  ==>  sum=6  diff=1019
            //                               j=7  ==>  sum=7  diff=1018
            //            i=2, k= 2, l=3077, j=0  ==>  sum=2  diff=3
            //                               j=1  ==>  sum=3  diff=2
            //                               j=2  ==>  sum=4  diff=1
            //                               j=3  ==>  sum=5  diff=0
            //                               j=4  ==>  sum=6  diff=1023
            //                               j=5  ==>  sum=7  diff=1022
            //                               j=6  ==>  sum=8  diff=1021
            //                               j=7  ==>  sum=9  diff=1020
            sum  += temp*array[(k+j)%size];
            diff += sign*temp*array[(l-j)%size];
            sign  = -sign;
		 }
         larray[i]   = sum;
         larray[i+1] = diff;
         k += 2;
         l += 2;
      }
      // The following case should only happen when size is odd
      // In this case, we want to keep the "extra" sample on the low-pass band
      // Assume that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      if(size_low!=size_high)
      {
		 size++;
         array[size_low-1] = array[size-1]*sqrt2;
      }
      for(i=0;i<size-1;i+=2)
      {
         j = i/2;
         array[j]          = larray[i];
         array[j+size_low] = larray[i+1];
      }
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}

/************************************************************************
 * Implements the inverse (or reconstruction) Daubechies-8
 * compact support 8-tap wavelet transform discussed above.
 ************************************************************************/
void wavelet_D8_rec(double *array,
                    int     size)
{
   register int i;
   register int j;
   int          lsize;
   int          count;

   // We first need to know how many levels of decomposition were performed
   // during the forward (decomposition) transform
   //   size<= 1 -> 0 levels
   // 1<size<= 2 -> 1 level
   // 2<size<= 4 -> 2 levels
   // 4<size<= 8 -> 3 levels
   // 8<size<=16 -> 4 levels ...
   count = 0;
   lsize = size;
   while(lsize>1)
   {
      lsize = lsize - lsize/2;
      count++;
   }
   // At this point, count = # of decomposition levels

   // For a given decomposition level, need to find out what is the
   // appropriate size to use.
   for(i=count-1;i>=0;i--)
   {
      lsize = size;
      for(j=0;j<i;j++)
      {
         lsize = lsize - lsize/2;
      }
      lsize = wavelet_D8_rec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of reconstruction, used when doing inverse 
 * Daubechies-8 compact support 8-tap wavelet transform
 * Daubechies-8 filter:   
 * [  0.230377813309, 0.714846570553, 0.630880767930, -0.027983769417,
 *   -0.187034811719, 0.030841381836, 0.032883011667, -0.010597401785 ] - LOW
 * [ -0.010597401785,-0.032883011667, 0.030841381836,  0.187034811719,
 *   -0.027983769417,-0.630880767930, 0.714846570553, -0.230377813309 ] - HIGH
 ************************************************************************/
int wavelet_D8_rec_1level(double *array,
                          int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int          size_low;
   int          size_high;
   double       temp1;
   double       temp2;
   double       temp3;
   double       temp4;
   double       sum;
   double       diff;
   double      *larray;
   double       sqrt2   = 1.41421356237309504880168872420969808;
   double       sqrt1_2 = 0.707106781186547524400844362104849039;
   double       D8[8] = { 0.230377813309,
                          0.714846570553,
                          0.630880767930,
                         -0.027983769417,
                         -0.187034811719,
                          0.030841381836,
                          0.032883011667,
                         -0.010597401785 };

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
	  // adjust size to be even
	  size = 2*size_high;
      for(i=0,k=3*size_high;i<size-1;i+=2,k++)
      {
         sum  = 0.0;
         diff = 0.0;
         for(j=0,l=0;j<8;j+=2,l++)
		 {
            temp1 = D8[j];
            temp2 = D8[j+1];
            temp3 = array[(k-l)%size_high];
            temp4 = array[(k+l)%size_high+size_low];
            // For example : i=0, k=1536, j=0, l=0, temp3=  0, temp4=512
            // (size=1024)                j=2, l=1, temp3=512, temp4=513
            //                            j=4, l=2, temp3=511, temp4=514
            //                            j=6, l=3, temp3=510, temp4=515
            //               i=2, k=1538, j=0, l=0, temp3=  2, temp4=513
            //                            j=2, l=1, temp3=  1, temp4=514
            //                            j=4, l=2, temp3=512, temp4=515
            //                            j=6, l=3, temp3=511, temp4=516
            sum  += temp1*temp3 + temp2*temp4;
            diff += temp2*temp3 - temp1*temp4;
		 }
         larray[i]   = sum;
         larray[i+1] = diff;
      }
      // The following case should only happen when size is odd
      // In this case, we have kept the "extra" sample on the low-pass band
      // assuming that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      // So, sample = sum/SQRT2
      if(size_low!=size_high)
      {
         size++;
         larray[i]   = array[size_low-1]*sqrt1_2;
      }
      memcpy(array,larray,size*sizeof(double));
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}

/************************************************************************
 * Implements the case filter = 4 discussed above,
 * i.e. Daubechies-16 compact support tap-16 wavelet transform
 ************************************************************************/
void wavelet_D16(double *array,
                 int     size,
                 int     direction)
{
   // 2 cases: Forward (direction=1) or Inverse (direction=-1) transform
   if(direction==1)
   {
      wavelet_D16_dec(array,size);
   }
   else if(direction==(-1))
   {
      wavelet_D16_rec(array,size);
   }
   else
   {
      // Do nothing
   }
}

/************************************************************************
 * Implements the forward (or decomposition) Daubechies-16 compact support
 * 16-tap wavelet transform discussed above.
 ************************************************************************/
void wavelet_D16_dec(double *array,
                     int     size)
{
   int lsize;

   lsize = size;
   while(lsize>1)
   {
      // The idea is that lsize returns the size of the low-pass band
      // For example, if size=512, then after 1-level of decomposition
      // (is the name of the function more clear now ?) lsize = 256
      // and so on.
      lsize = wavelet_D16_dec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of decomposition, used when doing Daubechies-16
 * compact support 16-tap wavelet transform
 * Daubechies-16 filter:
 * [ 0.054415842243,  0.312871590914,  0.675630736297,  0.585354683654,
 *  -0.015829105256, -0.284015542962,  0.000472484574,  0.128747426620,
 *  -0.017369301002, -0.044088253931,  0.013981027917,  0.008746094047,
 *  -0.004870352993, -0.000391740373,  0.000675449409, -0.000117476784 ]
 * "Allignment" : 
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 21 23 24 25
 *                                           .  .  .  .  .  .  .  .  .  .  .  .
 * .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .
 ************************************************************************/
int wavelet_D16_dec_1level(double *array,
                           int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int          size_low;
   int          size_high;
   double       temp;
   double       sum;
   double       diff;
   double       sign;
   double      *larray;
   double       sqrt2   = 1.41421356237309504880168872420969808;
   double       sqrt1_2 = 0.707106781186547524400844362104849039;
   double       D16[16] = { 0.054415842243,
                            0.312871590914,
                            0.675630736297,
                            0.585354683654,
                           -0.015829105256,
                           -0.284015542962,
                            0.000472484574,
                            0.128747426620,
                           -0.017369301002,
                           -0.044088253931,
                            0.013981027917,
                            0.008746094047,
                           -0.004870352993,
                           -0.000391740373,
                            0.000675449409,
                           -0.000117476784};

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
	  // adjust size to be even
	  size = 2*size_high;
      k = 0;
      l = 7*size+1;
      for(i=0;i<size-1;i+=2)
      {
         sum  =  0.0;
         diff =  0.0;
         sign = -1.0;
         for(j=0;j<16;j++)
         {
            temp  = D16[j];
            sum  += temp*array[(k+j)%size];
            diff += sign*temp*array[(l-j)%size];
            sign  = -sign;
		 }
         larray[i]   = sum;
         larray[i+1] = diff;
         k += 2;
         l += 2;
      }
      // The following case should only happen when size is odd
      // In this case, we want to keep the "extra" sample on the low-pass band
      // Assume that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      if(size_low!=size_high)
      {
		 size++;
         array[size_low-1] = array[size-1]*sqrt2;
      }
      for(i=0;i<size-1;i+=2)
      {
         j = i/2;
         array[j]          = larray[i];
         array[j+size_low] = larray[i+1];
      }
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}

/************************************************************************
 * Implements the inverse (or reconstruction) Daubechies-16
 * compact support 16-tap wavelet transform discussed above.
 ************************************************************************/
void wavelet_D16_rec(double *array,
                     int     size)
{
   register int i;
   register int j;
   int          lsize;
   int          count;

   // We first need to know how many levels of decomposition were performed
   // during the forward (decomposition) transform
   //   size<= 1 -> 0 levels
   // 1<size<= 2 -> 1 level
   // 2<size<= 4 -> 2 levels
   // 4<size<= 8 -> 3 levels
   // 8<size<=16 -> 4 levels ...
   count = 0;
   lsize = size;
   while(lsize>1)
   {
      lsize = lsize - lsize/2;
      count++;
   }
   // At this point, count = # of decomposition levels

   // For a given decomposition level, need to find out what is the
   // appropriate size to use.
   for(i=count-1;i>=0;i--)
   {
      lsize = size;
      for(j=0;j<i;j++)
      {
         lsize = lsize - lsize/2;
      }
      lsize = wavelet_D16_rec_1level(array,lsize);
   }
}

/************************************************************************
 * Implements 1-level of reconstruction, used when doing inverse 
 * Daubechies-16 compact support 16-tap wavelet transform
 * Daubechies-8 filter:   
 * [ 0.054415842243,  0.312871590914,  0.675630736297,  0.585354683654,
 *  -0.015829105256, -0.284015542962,  0.000472484574,  0.128747426620,
 *  -0.017369301002, -0.044088253931,  0.013981027917,  0.008746094047,
 *  -0.004870352993, -0.000391740373,  0.000675449409, -0.000117476784 ]
 ************************************************************************/
int wavelet_D16_rec_1level(double *array,
                           int     size)
{
   register int i;
   register int j;
   register int k;
   register int l;
   int          size_low;
   int          size_high;
   double       temp1;
   double       temp2;
   double       temp3;
   double       temp4;
   double       sum;
   double       diff;
   double      *larray;
   double       sqrt2   = 1.41421356237309504880168872420969808;
   double       sqrt1_2 = 0.707106781186547524400844362104849039;
   double       D16[16] = { 0.054415842243,
                            0.312871590914,
                            0.675630736297,
                            0.585354683654,
                           -0.015829105256,
                           -0.284015542962,
                            0.000472484574,
                            0.128747426620,
                           -0.017369301002,
                           -0.044088253931,
                            0.013981027917,
                            0.008746094047,
                           -0.004870352993,
                           -0.000391740373,
                            0.000675449409,
                           -0.000117476784};

   if(size>1)
   {
      larray = (double *) malloc(size*sizeof(double));
      if(larray==NULL)
      {
         printf("Not enough memory (%d bytes) for %s - exiting\n",
                size*sizeof(double),"larray");
         exit(-1);
      }
      size_high = size/2;
      // size_high = size/2 (integer arithmetic) should give us the size of the
      // HIGH-PASS filtered band, i.e. the differences. 
      size_low  = size - size_high;
      // size of the LOW-PASS filtered band, i.e. the sums.
	  // adjust size to be even
	  size = 2*size_high;
      for(i=0,k=7*size_high;i<size-1;i+=2,k++)
      {
         sum  = 0.0;
         diff = 0.0;
         for(j=0,l=0;j<16;j+=2,l++)
		 {
            temp1 = D16[j];
            temp2 = D16[j+1];
            temp3 = array[(k-l)%size_high];
            temp4 = array[(k+l)%size_high+size_low];
            sum  += temp1*temp3 + temp2*temp4;
            diff += temp2*temp3 - temp1*temp4;
		 }
         larray[i]   = sum;
         larray[i+1] = diff;
      }
      // The following case should only happen when size is odd
      // In this case, we have kept the "extra" sample on the low-pass band
      // assuming that the "missing" sample is = to the last. In this case,
      // sum = 2*(the last)/sqrt2 = sqrt2*(the last)
      // diff = 0 (no need to store !!!)
      // So, sample = sum/SQRT2
      if(size_low!=size_high)
      {
         size++;
         larray[i]   = array[size_low-1]*sqrt1_2;
      }
      memcpy(array,larray,size*sizeof(double));
      free(larray);
      return(size_low);
   }
   else
   {
      return(size);
   }
}