#include <stdio.h>

/**
   \file

   This file contains code to implement Lifting Scheme wavelets.
   Lifting Scheme wavelets are described in Wim Sweldens' tutorial
   paper <i>Building Your Own Wavelets at Home</i> which is available
   on the Web.

   Lifting Scheme wavelets are a conceptual extension of Haar wavelets.
   One of the disadvantages of Haar wavelets is the high frequency
   (largest) coefficient spectrum can miss detail (even to odd
   transitions, for example).  Lifting Scheme wavelets properly
   represent change in all coefficient spectrum.  This makes lifting
   scheme wavelets a better choice for some algorithms which do
   filtering based on the absolute standard deviation calculated on the
   high frequency coefficients.

 */


/**
   This class contains functions that are useful in debugging
   wavelet code.
 */
class debug {
   public:
   
   /**
      Print the result of a wavelet transform so that each coefficient
      spectrum is printed in its own block, showing the structure of
      the result.
   */
   void pr_ordered( const DoubleVec& coef)
   {
      const char *format = "%7.4f ";
      printf("{");
      int len = coef.size();
      int cnt = 0;
      int num_in_freq = 0;
      for (int i = 0; i < len; i++) {
         printf(format, coef[i]);
         cnt++;
         if (num_in_freq == 0) {
            printf("\n");
            cnt = 0;
            num_in_freq = 1;
         }
         else if (cnt == num_in_freq) {
            printf("\n");
            cnt = 0;
            num_in_freq = num_in_freq * 2;
         }
      }
      printf("}\n");
   } // pr_ordered

   /**
      Print wavelet transform intermediate results.

      The wavelet transform calculates a set of averages and
      coefficients (differences).  These are calculated recursively.
      For example, for a time series with size 32, the first iteration
      of the recursive wavelet transform will calculate a set of 16
      averages and 16 coefficients (differences).

      The next set of differences and avarages will be calculated at
      the start of the array (0..15) and will consist of 8 averages
      and 8 differences.

      This function is passed a value for <tt>N</tt> which is the size
      of the average and difference spectrum.  It prints each of these.

   */
   void pr_intermediate( DoubleVec& ts, int N)
   {
      const char* fmt = "%7.4f ";
      int i;
      int half = N >> 1;
      
      printf(" ave = ");
      for (i = 0; i < half; i++) {
         printf(fmt, ts[i] );
      }
      printf("\n");
      printf("coef = ");
      for (i = half; i < N; i++) {
         printf(fmt, ts[i] );
      }
      printf("\n");
   } // pr_intermediate


   /**
      Print the first N elements of <tt>ts</tt>.  If N is less than
      the length of <tt>ts</tt> the print a dividing mark (e.g., "|")
      and print the rest of the array.

      This function is useful for displaying the flow of the
      wavelet computation.
   */
   void pr_vector( DoubleVec& ts)
   {
     const int N = ts.size();
     const char* fmt = "%7.4f ";

     int i;
     for (i = 0; i < N; i++) {
       printf(fmt, ts[i] );
     }
     int len = ts.size();
     if (len != N) {
       printf(" | ");
       
       for (i = N; i < len; i++) {
	 printf(fmt, ts[i] );
       }
     }
     printf("\n");
   } // pr_vector

};  // debug





/**
   Support for four point polynomial interpolation, using the Lagrange
   formula.

 */
class interpolation {

   private:
   double fourPointTable[4][4];
   double twoPointTable[4][4];

   /**
     <p>
     The polynomial interpolation algorithm assumes that the known
     points are located at x-coordinates 0, 1,.. N-1.  An interpolated
     point is calculated at <b><i>x</i></b>, using N coefficients.  The
     polynomial coefficients for the point <b><i>x</i></b> can be 
     calculated staticly, using the Lagrange method.
     </p>

     @param x the x-coordinate of the interpolated point
     @param N the number of polynomial points.
     @param c[] an array for returning the coefficients

     */
  void lagrange( double x, int N, double c[] )
  {
    double num, denom;

    for (int i = 0; i < N; i++) {
      num = 1;
      denom = 1;
      for (int k = 0; k < N; k++) {
	if (i != k) {
	  num = num * (x - k);
	  denom = denom * (i - k);
	}
      } // for k
      c[i] = num / denom;
    } // for i
  } // lagrange


  /**
    <p>
    For a given N-point polynomial interpolation, fill the coefficient
    table, for points 0.5 ... (N-0.5).
    </p>

   */
  void fillTable( const int N, double table[4][4] )
  {
    double x;
    double n = N;
    int i = 0;

    for (x = 0.5; x < n; x = x + 1.0) {
      lagrange( x, N, table[i] );
      i++;
    }
  } // fillTable


  /**
    Print an N x N table polynomial coefficient table
   */
  void printTable( double table[4][4], int N)
  {
    printf("%d-point interpolation table:\n", N);
    double x = 0.5;
    for (int i = 0; i < N; i++) {
      printf("%4.2f: ", x);
      for (int j = 0; j < N; j++) {
	printf("%6.4f", table[i][j] );
	if (j < N-1)
	  printf(", ");
      }
      printf("\n");
      x = x + 1.0;
    }    
  }


  /**
    Print the 4-point and 2-point polynomial coefficient
    tables.
   */
  void printTables()
  {
    printTable( fourPointTable, 4 );
    printTable( twoPointTable, 2 );
    printf("\n");
  } // printTables


    /**
    <p>
    For the polynomial interpolation point x-coordinate 
    <b><i>x</i></b>, return the associated polynomial
    interpolation coefficients.
    </p>

    @param x the x-coordinate for the interpolated pont
    @param n the number of polynomial interpolation points
    @param c[] an array to return the polynomial coefficients
   */
  void getCoef( double x, int n, double c[])
  {
    int j = (int)x;

    if (j < 0 || j >= n) {
      printf("poly::getCoef: n = %d, bad x value = %f\n", n, x);
    }

    if (n == 2) {
      c[2] = 0.0;
      c[3] = 0.0;
    }
    else if (n != 4) {
      printf("poly::getCoef: bad value for N\n");
      return;
    }

    for (int i = 0; i < n; i++) {
      if (n == 4) {
	c[i] = fourPointTable[j][i];
      }
      else { // n == 2
	c[i] = twoPointTable[j][i];
      }
    }

  } // getCoef


   public:

   /**
      The interpolation class constructor calculates the coefficients for
      a four point polynomial interpolated at -0.5, 0.5, 1.5, 2.5 and 3.5.
   */
   interpolation() 
   {
    // Fill in the 4-point polynomial interplation table
    // for the points 0.5, 1.5, 2.5, 3.5
    fillTable( 4, fourPointTable );
    
    // Fill in the 2-point polynomial interpolation table
    // for 0.5 and 1.5
    fillTable( 2, twoPointTable );
   } // interpolation class constructor


  /**
    <p>
    Given four points at the x,y coordinates {0,d<sub>0</sub>},
    {1,d<sub>1</sub>}, {2,d<sub>2</sub>}, {3,d<sub>3</sub>}
    return the y-coordinate value for the polynomial interpolated
    point at <b><i>x</i></b>.
    </p>
    
    @param x the x-coordinate for the point to be interpolated
    @param N the number of interpolation points
    @param d[] an array containing the y-coordinate values for
           the known points (which are located at x-coordinates
	   0..N-1).
   */
  double interpPoint( double x, int N, double d[])
  {
    double c[4];
    double point = 0;
    
    int n = 4;
    if (N < 4)
      n = N;

    getCoef( x, n, c );

    if (n == 4) {
      point = c[0]*d[0] + c[1]*d[1] + c[2]*d[2] + c[3]*d[3]; 
    }
    else if (n == 2) {
      point = c[0]*d[0] + c[1]*d[1];
    }

    return point;
  } // interpPoint

};  // class interpolation






/**
   This class implements the wavelet Lifting Scheme with polynomial
   interpolation.

   This version of wavelets is based on Wim Sweldens' tutorial paper
   <i>Building Your Own Wavelets at Home</i>.  This tutorial was
   presented at SIGGraph.  The tutorial is available on the 
   Web.

   One of the attractive features of Lifting Scheme wavelets is that
   the transform and the inverse transform are exact mirrors of each other.
   To arrive at the inverse transform the order of the operations is
   reversed and subtraction and addition operations are exchanged.
   This allows a generic Lifting Scheme base class to be constructed,