Main Page   Namespace List   Class Hierarchy   Compound List   File List   Compound Members   File Members  

daub.h

Go to the documentation of this file.
00001 
00002 #include <math.h>
00003 
00131 template <class T>
00132 class Daubechies : public liftbase<T, double> {
00133   
00134 protected:
00135   void predict( T& vec, int N, transDirection direction )
00136   {
00137     assert( false );
00138   } // predict
00139   
00140   virtual void update( T& vec, int N, transDirection direction )
00141   {
00142     assert( false );
00143   } // update
00144   
00145 private:
00147   double h0, h1, h2, h3;
00149   double g0, g1, g2, g3;
00150   
00151   double Ih0, Ih1, Ih2, Ih3;
00152   double Ig0, Ig1, Ig2, Ig3;
00153 
00154 public:
00155   
00160   void forwardStep( T& a, const int n )
00161   {
00162     if (n >= 4) {
00163       int i, j;
00164       const int half = n >> 1;
00165       
00166       double* tmp = new double[n];
00167       
00168       for (i = 0, j = 0; j < n-3; j += 2, i++) {
00169         tmp[i]      = a[j]*h0 + a[j+1]*h1 + a[j+2]*h2 + a[j+3]*h3;
00170         tmp[i+half] = a[j]*g0 + a[j+1]*g1 + a[j+2]*g2 + a[j+3]*g3;
00171       }
00172       
00173       tmp[i]      = a[n-2]*h0 + a[n-1]*h1 + a[0]*h2 + a[1]*h3;
00174       tmp[i+half] = a[n-2]*g0 + a[n-1]*g1 + a[0]*g2 + a[1]*g3;
00175       
00176       for (i = 0; i < n; i++) {
00177         a[i] = tmp[i];
00178       }
00179       delete [] tmp;
00180     }
00181   }
00182   
00188   void forwardStepRev( T& a, const int n )
00189   {
00190     if (n >= 4) {
00191       int i, j;
00192       const int half = n >> 1;
00193       
00194       double* tmp = new double[n];
00195       
00196       for (i = 0, j = 0; j < n-3; j += 2, i++) {
00197         tmp[i+half] = a[j]*h0 + a[j+1]*h1 + a[j+2]*h2 + a[j+3]*h3;
00198         tmp[i]      = a[j]*g0 + a[j+1]*g1 + a[j+2]*g2 + a[j+3]*g3;
00199       }
00200       
00201       tmp[i+half] = a[n-2]*h0 + a[n-1]*h1 + a[0]*h2 + a[1]*h3;
00202       tmp[i]      = a[n-2]*g0 + a[n-1]*g1 + a[0]*g2 + a[1]*g3;
00203       
00204       for (i = 0; i < n; i++) {
00205         a[i] = tmp[i];
00206       }
00207       delete [] tmp;
00208     }
00209   }
00210   
00214   void inverseStep( T& a, const int n )
00215   {
00216     if (n >= 4) {
00217       int i, j;
00218       const int half = n >> 1;
00219       const int halfPls1 = half + 1;
00220       
00221       double* tmp = new double[n];
00222       
00223       //      last smooth val  last coef.  first smooth  first coef
00224       tmp[0] = a[half-1]*Ih0 + a[n-1]*Ih1 + a[0]*Ih2 + a[half]*Ih3;
00225       tmp[1] = a[half-1]*Ig0 + a[n-1]*Ig1 + a[0]*Ig2 + a[half]*Ig3;
00226       for (i = 0, j = 2; i < half-1; i++) {
00227         //     smooth val     coef. val       smooth val    coef. val
00228         tmp[j++] = a[i]*Ih0 + a[i+half]*Ih1 + a[i+1]*Ih2 + a[i+halfPls1]*Ih3;
00229         tmp[j++] = a[i]*Ig0 + a[i+half]*Ig1 + a[i+1]*Ig2 + a[i+halfPls1]*Ig3;
00230       }
00231       for (i = 0; i < n; i++) {
00232         a[i] = tmp[i];
00233       }
00234       delete [] tmp;
00235     }
00236   } // inverseStep
00237   
00238 
00243   Daubechies() 
00244   {
00245     const double sqrt_3 = sqrt( 3 );
00246     const double denom = 4 * sqrt( 2 );
00247     
00248     //
00249     // forward transform scaling (smoothing) coefficients
00250     //
00251     h0 = (1 + sqrt_3)/denom;
00252     h1 = (3 + sqrt_3)/denom;
00253     h2 = (3 - sqrt_3)/denom;
00254     h3 = (1 - sqrt_3)/denom;
00255     //
00256     // forward transform wavelet coefficients (a.k.a. high
00257     // pass filter coefficients)
00258     //
00259     g0 =  h3;
00260     g1 = -h2;
00261     g2 =  h1;
00262     g3 = -h0;
00263     
00264     Ih0 = h2;
00265     Ih1 = g2;  // h1
00266     Ih2 = h0;
00267     Ih3 = g0;  // h3
00268     
00269     Ig0 = h3;
00270     Ig1 = g3;  // -h0
00271     Ig2 = h1;
00272     Ig3 = g1;  // -h2
00273   }
00274 
00278   void forwardTrans( T& ts, int N )
00279   {
00280     int n;
00281     for (n = N; n >= 4; n >>= 1) {
00282       forwardStep( ts, n );
00283     }
00284   } // forwardTrans
00285   
00286   
00290   void inverseTrans( T& coef, int N )
00291   {
00292     int n;
00293     for (n = 4; n <= N; n <<= 1) {
00294       inverseStep( coef, n );
00295     }
00296   } // inverseTrans
00297   
00298 }; // Daubechies
00299 

Generated at Thu May 22 21:12:35 2003 for Hurst Exponent Calculation and Supporting Statistics by doxygen1.2.8.1 written by Dimitri van Heesch, © 1997-2001