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 }
00139
00140 virtual void update( T& vec, int N, transDirection direction )
00141 {
00142 assert( false );
00143 }
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
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
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 }
00237
00238
00243 Daubechies()
00244 {
00245 const double sqrt_3 = sqrt( 3 );
00246 const double denom = 4 * sqrt( 2 );
00247
00248
00249
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
00257
00258
00259 g0 = h3;
00260 g1 = -h2;
00261 g2 = h1;
00262 g3 = -h0;
00263
00264 Ih0 = h2;
00265 Ih1 = g2;
00266 Ih2 = h0;
00267 Ih3 = g0;
00268
00269 Ig0 = h3;
00270 Ig1 = g3;
00271 Ig2 = h1;
00272 Ig3 = g1;
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 }
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 }
00297
00298 };
00299