// $Id: Wavelet.cc,v 0.1 2001/08/10 19:10:00 klimenko Exp $

#define WAVEDWT_CC

#include "Wavelet.hh"

namespace datacondAPIwat {

// constructors

Wavelet::Wavelet(int mH=1, int mL=1, int tree=0, enum BORDER border=B_CYCLE) : 
   m_WaveType(HAAR)
{ 
   m_H = mH;
   m_L = mL;
   m_Border = border;
   m_Level = 0;
   m_TreeType = tree;
}

Wavelet::Wavelet(const Wavelet &w) 
{ 
   m_H = w.m_H;
   m_L = w.m_L;
   m_Border = w.m_Border;
   m_Level = w.m_Level;
   m_TreeType = w.m_TreeType;
   m_WaveType = w.m_WaveType;
}

// destructor
Wavelet::~Wavelet()
{ }

Wavelet* Wavelet::Clone() const
{
  return new Wavelet(*this);
}

//*******************************
//*  wavedata acess functions   *
//*******************************

int Wavelet::getOffset(int level, int layer)
{       
   int n=0;

   for(int i=0; i<level; i++)
      if((layer>>i)&1) n += 1<<(level-1-i);

   return n;
}

int Wavelet::convertL2F(int level, int layer)
{       
   int n = layer;
   int j;
   for(int i=1; i<level; i++) {
      j = (1<<i) & (n);
      if(j) n = ((1<<i)-1) ^ (n);
   }
   
   return n;
}

int Wavelet::convertF2L(int level, int index)
{       
   int n = index;
   int j;
   for(int i=level-1; i>=1; i--) {
      j = ((1<<i) & (n));
      if(j) n = ((1<<i)-1) ^ (n);
   }
   return n;
}

//*******************************
//*  WAT functions              *
//*******************************

// Calculates polynomial interpolation coefficient using Lagrange formula.
inline double Lagrange(const int n, const int i, const double x)
{
    double c = 1.;
    double xn = x+n/2.-0.5;	// shift to the center of interpolation interval

    for(int j=0; j<n; j++) 
       if(j!=i) c *= (xn-j)/(i-j);

    return c;
}

// Polynomial interpolation using Neville's alghorithm.
double Nevill(const int n, 
	      const double* pX, 
	      const double* pY, 
	      const double x)
{
   if(n <= 1){
      return *pY;
   }
   int m = n-1;
   double PL = Nevill(m, pX, pY, x);
   double PR = Nevill(m, pX+1, pY+1, x);
   return ((x - *(pX+m))*PL + (*pX - x)*PR)/(*pX - *(pX+m));
}

} // namespace datacondAPIwat