var/external/wavelet/dwtai2_w.c

   1 // Author: Aime' Fournier
   2 // E-mail: fournier@ucar.edu
   3
   4 #include "qf.h"                 // for pqf ...
   5
   6 // See 2009/8/4 e-mail from john@michalakes.us
   7 #ifndef CRAY
   8 #   ifdef NOUNDERSCORE
   9 #      define DWTAI2_W dwtai2_w
  10 #   else
  11 #      ifdef F2CSTYLE
  12 #         define DWTAI2_W dwtai2_w__
  13 #      else
  14 #         define DWTAI2_W dwtai2_w_
  15 #      endif
  16 #   endif
  17 #endif
  18
  19 void DWTAI2_W (char *nam,       // in: filter name
  20            int *ran,            // in: filter length
  21            real *u,             // in/out: data & transform
  22            real *s,             // scratch space
  23            int *iv,             // v-space lengths
  24            int *ia,             // w-space starts
  25            int *lm) {           // maximum level
  26   extern real l2norm();         // l2 norm
  27   int *ja[2],*jv[2];            // reshapes of ia & iv
  28   pqf *g,*h;                    // wavelet hpf & lpf
  29 //real r;                       // only for printf(l2norm())
  30
  31 //printf("dwtai2_w::dwtai2_w(%s,%2d,{",nam,*ran);
  32 //for( int i=0; i<4*ia[*lm+1]; i+=ia[*lm+1] ) {
  33 //   printf("{");
  34 //   for( int j=i; j<i+4; j++ ) printf("%11.4e,",u[j]);
  35 //   printf("\b...},");
  36 //}
  37 //printf("\b...},{%6.3f,...},{%2d,...},{%2d,...},%2d)\n",*s,*iv,*ia,*lm);
  38
  39   for( int i=0; i<2; i++ ) {
  40      ja[i]=ia+(*lm+1)*i;
  41      jv[i]=iv+(*lm+1)*i;
  42 //   for( int j=0; j<=*lm; j++ )
  43 //      printf("dwtai2_w: {ia,iv}[%d][%d]={%d,%d}\n",i,j,ja[i][j],jv[i][j]);
  44   }
  45
  46   h=qf(nam,*ran,0);             // assign h
  47   g=qf(nam,*ran,1);             // assign g
  48 //PrintFilter(nam,*ran,0,h);
  49 //PrintFilter(nam,*ran,1,g);
  50
  51 //r=l2norm(u,ja[1][0],jv[0][0],jv[1][0]);
  52 //printf("dwtai2_w:          ||u||=%15.9E\n",r);
  53   dwtai2(u,s,jv,ja,*lm,h,g);
  54 //printf("dwtai2_w: 1-||Wu||/||u||=%15.9E\n",1.-l2norm(u,ja[1][0],ja[0][0],ja[1][0])/r);
  55
  56   free(g->fp);free(g);free(h->fp);free(h);
  57 }