Wavelet.cc

00001 #include "Wavelet.hpp"
00002 #include <string.h>
00003 
00004 //extern int debugRefine ;
00005 
00006 void GetGeneralBoundaryConditions(int BC[2] , int *BG) {
00007   BG[0]=BG[1]=-1 ;
00008   if (BC[1]==PERIODIC) BG[1]=PERIODIC ;
00009 }
00010 
00011 void out_of_memory() {
00012   cerr << "out of memory error!" << endl;
00013   exit(-1) ;
00014 }
00015 
00016 void BoundaryCorrection(int *BC, double *d, int l) {
00017  if (BC[0]==0             ) d[0   ]=0 ;
00018  if (BC[1]==0             ) d[1<<l]=0 ;
00019 
00020  if (BC[0]==1             ) d[1       ]=0 ;
00021  if (BC[1]==1             ) d[(1<<l)-1]=0 ;
00022 
00023  if (BC[1]==PERIODIC) d[1<<l]=0 ;
00024 }
00025 
00026 void FD2stencil::Set(int which, int wdth, int *js, double *ws) {
00027   width[which]=wdth; 
00028   for (int i=0; i<wdth; i++) {
00029     j[which][i]=js[i] ;
00030     w[which][i]=ws[i] ;
00031   }
00032 }
00033 
00034 void FD2stencil::Mirror(int which, int from, double flip) {
00035   assert(which != from) ;
00036   width[which]=width[from] ;
00037   for (int i=0; i<width[from]; i++) {
00038     j[which][i]=-j[from][width[from]-1-i] ;
00039     w[which][i]= w[from][width[from]-1-i]*flip ;
00040   }
00041 }
00042 
00043 /* READLINE reads a complete line from FILE *fid to char *line.
00044    If (necessary) white spaces and empty lines are skipped until a  
00045    non-empty lines is found. This one is read. */
00046 void READLINE(FILE *fid,char *line) {
00047   int i=0,f=0,j=0 ;
00048   do{ 
00049     if (fscanf(fid,"%c",&(line[i]))==EOF) break ;
00050     if ((line[i]!=10)&&(line[i]!=32)) f=1 ;
00051     j=i ; i=i+f ;
00052   } while((line[j]!=10)||(f==0));
00053   line[i]=0 ;
00054 }
00055 
00056 void Wavelets::ReadInnerEntries(FILE *file,double *a,int *U,int *O) {
00057   fscanf(file,"%d  %d",U,O) ;
00058   if ((*U)+(*O)>=0) 
00059     for (int i=0 ;i<=(*U)+(*O) ;i++) { fscanf(file,"%le ",&a[i]) ;}
00060 }
00061 
00062 
00063 Wavelets::Wavelets(const char *typ,const int N) {
00064 
00065   std::set_new_handler(out_of_memory) ;
00066 
00067   int flag=1 ;
00068   InterpoletFlag=LiftingFlag=false ;
00069   Buffer=NULL ;
00070 
00071   flag =GetCoefficients(typ,"left" ,N) ;
00072   flag&=GetCoefficients(typ,"right",N) ;
00073   if (!flag) { std::cout<< "Could not get wavelet coefficients\n" ; exit(-1) ;}
00074 }
00075 
00076 Wavelets::~Wavelets() {;}
00077 
00078 int Wavelets::GetCoefficients(const char *typ,const int NN)  {
00079   return GetCoefficients(typ,"left" ,NN) & GetCoefficients(typ,"right",NN) ;
00080 }
00081 
00082 int Wavelets::GetCoefficients(const char *typ,const char *leftright,const int NN) {
00083   FILE *file       ;
00084   int i,j,n,left=0 ;
00085   bool found=false ;
00086   char line[10000] ; 
00087   char typs[100]   ;
00088   char lfrg[100]   ;
00089   
00090   int KK,LL,RR,KWW ;
00091   
00092   int op,bcT,bcF,iswT,iswF ;
00093   
00094   if (strcmp(typ,"Interpolet")==0) { return GetICoefficients(typ,leftright,NN) ; }
00095   
00096   if (strncmp(leftright,"left",4)==0) left=1 ;
00097   
00098   char name[1000] ;
00099   sprintf(name,"%s/MATLAB/Data/coeff.dat",_WAVELET_ROOT_) ;
00100   file=fopen(name,"r") ;/*assert(file) ;*/
00101   if (file == NULL) {std::cout<< "Coeff-file not found: "<<name<<'\n'; assert(0);} 
00102   
00103   do {
00104     READLINE(file,line) ;
00105     if (line[0]==0) {std::cout<< "Filter coefficients for "<<typ<<'('<<NN<<") not found\n"; assert(0);}
00106     
00107     if (strncmp(line,"begin",5)==0) {
00108       fscanf(file,"%s",typs) ;                 // Wavelet Family
00109       fscanf(file,"%s",lfrg) ;                 // left/right
00110       READLINE(file,line);  // comments
00111       fscanf(file,"%d %d %d %d %d",&n,&LL,&RR,&KK,&KWW) ;
00112       
00113       if ((strcmp(typs,typ)==0) && (n==NN) && (strcmp(leftright,lfrg)==0)) {
00114         
00115         found=true ;
00116         strcpy(type,&typs[0]) ;
00117         N=NN                  ;
00118         
00119         if (left) {
00120           L=LL ;
00121           R=RR ;
00122           K0=KK   ;
00123           KW0=KWW ;
00124           MX=K0+L+R-1 ;
00125           
00126           /*********************/
00127           /*                   */
00128           /*      HJ / GJ      */
00129           READLINE(file,line) ;
00130           ReadInnerEntries(file,HJ.a,&HJ.U,&HJ.O) ;
00131           for (i=0;i<2;i++) {
00132             HJ.IsWT=0 ;
00133             READLINE(file,line) ;
00134             HJ.AL[i].Read(file) ;
00135           }
00136           
00137           READLINE(file,line) ;
00138           ReadInnerEntries(file,GJ.a,&GJ.U,&GJ.O) ;
00139           for (i=0;i<2;i++) {
00140             GJ.IsWT=1 ;
00141             READLINE(file,line) ;
00142             GJ.AL[i].Read(file) ;
00143           }
00144           /*****************************/
00145           /*   Quadratur-Matrizen      */
00146           /*                           */
00147           READLINE(file,line) ;
00148           ReadInnerEntries(file,W.a,&W.U,&W.O) ;
00149           for (i=0;i<2;i++) {
00150             READLINE(file,line) ;
00151             W.AL[i].Read(file) ;
00152           }
00153           
00154           READLINE(file,line) ;
00155           ReadInnerEntries(file,NULL,&i,&j) ;READLINE(file,line) ;
00156           for (i=0;i<2;i++) {
00157             READLINE(file,line) ;
00158             IW.AL[i].Read(file) ;
00159           }
00160           
00161           READLINE(file,line) ;
00162           ReadInnerEntries(file,Phi.a,&Phi.U,&Phi.O) ;
00163           for (i=0;i<2;i++) {
00164             READLINE(file,line) ;
00165             Phi.AL[i].Read(file) ;
00166           }
00167           /********************************/
00168           /*     Projektions-Matrizen     */
00169           /*                              */
00170           for (i=0;i<1;i++) {
00171             READLINE(file,line) ;
00172             OP0[i].Read(file) ;
00173           }         
00174           
00175           /********************************/
00176           /*     Operator-Matrizen        */
00177           /*                              */
00178           for (op=0;op<2;op++)
00179             for (iswT=0;iswT<2;iswT++)
00180               for (iswF=0;iswF<2;iswF++) {   
00181                 Operators[op][iswT][iswF].IsWT =iswT ;
00182                 Operators[op][iswT][iswF].IsWF =iswF ;
00183                 Operators[op][iswT][iswF].power=op+1 ;
00184                 
00185                 READLINE(file,line) ;
00186                 ReadInnerEntries(file,Operators[op][iswT][iswF].a,&Operators[op][iswT][iswF].U,&Operators[op][iswT][iswF].O) ;
00187                 for (bcT=0;bcT<2;bcT++)
00188                   for (bcF=0;bcF<2;bcF++) {
00189                     READLINE(file,line) ;
00190                     Operators[op][iswT][iswF].AL [bcT][bcF].Read(file) ;
00191                     Operators[op][iswT][iswF].ALT[bcT][bcF].Read(file) ;
00192                   }
00193                 
00194               } 
00195           /****************************/
00196           /*     ATau-Matrizen        */
00197           /*                          */
00198           for (i=0;i<2;i++) {
00199             READLINE(file,line) ;
00200             ATau0[i].Read(file) ;
00201           }
00202         } else {  // right
00203           
00204           K1 =KK  ;
00205           KW1=KWW ;
00206           MX =MX+K1 ;
00207           
00208           /*********************/
00209           /*                   */
00210           /*      HJ / GJ      */
00211           READLINE(file,line) ;
00212           ReadInnerEntries(file,NULL,&i,&j) ;READLINE(file,line) ;
00213           for (i=0;i<2;i++) {
00214             READLINE(file,line) ;
00215             HJ.AR[i].Read(file) ;
00216           }
00217           
00218           READLINE(file,line) ;
00219           ReadInnerEntries(file,NULL,&i,&j) ;READLINE(file,line) ;
00220           for (i=0;i<2;i++) {
00221             READLINE(file,line) ;
00222             GJ.AR[i].Read(file) ;
00223           }
00224           
00225           /*****************************/
00226           /*   Quadratur-Matrizen      */
00227           /*                           */
00228           READLINE(file,line) ;
00229           ReadInnerEntries(file,NULL,&i,&j) ;READLINE(file,line) ;
00230           for (i=0;i<2;i++) {
00231             READLINE(file,line) ;
00232             W.AR[i].Read(file) ;
00233           }
00234           
00235           READLINE(file,line) ;
00236           ReadInnerEntries(file,NULL,&i,&j) ;READLINE(file,line) ;
00237           for (i=0;i<2;i++) {
00238             READLINE(file,line) ;
00239             IW.AR[i].Read(file) ;
00240           }
00241           
00242           READLINE(file,line) ;
00243           ReadInnerEntries(file,NULL,&i,&j) ;READLINE(file,line) ;
00244           for (i=0;i<2;i++) {
00245             READLINE(file,line) ;
00246             Phi.AR[i].Read(file) ;
00247           }
00248           /********************************/
00249           /*     Projektions-Matrizen     */
00250           /*                              */
00251           for (i=0;i<1;i++) {
00252             READLINE(file,line) ;
00253             OP1[i].Read(file) ;
00254           }         
00255           /********************************/
00256           /*     Operator-Matrizen        */
00257           /*                              */
00258           for (op=0;op<2;op++)
00259             for (iswT=0;iswT<2;iswT++)
00260               for (iswF=0;iswF<2;iswF++) {
00261                 
00262                 READLINE(file,line) ;
00263                 ReadInnerEntries(file,NULL,&i,&j) ;READLINE(file,line) ;
00264                 for (bcT=0;bcT<2;bcT++)
00265                   for (bcF=0;bcF<2;bcF++) {
00266                     READLINE(file,line) ;
00267                     Operators[op][iswT][iswF].AR [bcT][bcF].Read(file) ;
00268                     Operators[op][iswT][iswF].ART[bcT][bcF].Read(file) ;
00269                   }
00270               } 
00271           
00272           /****************************/
00273           /*     ATau-Matrizen        */
00274           /*                          */
00275           for (i=0;i<2;i++) {
00276             READLINE(file,line) ;
00277             ATau1[i].Read(file) ;
00278           }
00279           
00280           /****************************/
00281           /*     duale Trafo Matrizen */
00282           /*                          */
00283           HtJ.Copy(&HJ) ;
00284           GtJ.Copy(&GJ) ;
00285           
00286         }        // left right 
00287       }          // end if typ/quad ...
00288     }            // end if found begin
00289   } while ((strncmp(&line[0],"stop",4)!=0) && (!found) ) ;
00290   
00291   if (!found) {
00292     L=R=N=0 ;
00293   }
00294   
00295   MX1=MX ;
00296   Level0=0 ;
00297   i=MX1 ;
00298   while ((i=(i>>1)) > 0) { Level0++ ;}
00299   if (N==3) Level0=4 ;
00300   if (N==1) Level0=1 ;
00301   
00302   fclose(file) ;
00303   
00304   if (left) ReadFDMatrices() ; 
00305   
00306   return(found ? 1 : 0) ;
00307 }
00308 
00309 int Wavelets::GetICoefficients(const char *typ,const char *leftright,const int NN)
00310 {
00311   FILE *file         ;
00312   int i,n,left=0 ;
00313   int found=0    ;
00314   char line[10000]    ; 
00315   char typs[100]     ;
00316   char lfrg[100]     ;
00317   
00318   int KK,LL,RR,KWW, op,bcF,bcT,iswT,iswF ;
00319   
00320   InterpoletFlag=true ;
00321   
00322   if (strncmp(leftright,"left",4)==0) left=1 ;
00323   
00324   char name[1000] ;
00325   sprintf(name,"%s/MATLAB/Data/Icoeff.dat",_WAVELET_ROOT_) ;
00326   file=fopen(name,"r") ; 
00327   if (file == NULL) {std::cout<< "Coeff-file not found: "<<name<<'\n'; return 0;} 
00328   
00329   do {
00330     READLINE(file,line) ;
00331     if (line[0]==0) {std::cout<< "Filter coefficients for "<<typ<<'('<<NN<<") not found\n"; assert(0);}
00332     
00333     if (strncmp(line,"begin",5)==0) {
00334       fscanf(file,"%s",typs) ;  /* Wavelet Family  */
00335       fscanf(file,"%s",lfrg) ;  /* left/right      */
00336       READLINE(file,line)    ;  /* comments         */
00337       fscanf(file,"%d %d %d %d %d %d",&n,&LL,&RR,&KK,&KWW , &NBC) ;
00338       if ((strcmp(typs,typ)==0) && (n==NN) && (strcmp(leftright,lfrg)==0)) {
00339         found=1 ;
00340         strcpy(type,&typs[0]) ;
00341         N=NN                  ;
00342         
00343         if (left) {
00344           
00345           L=LL ;
00346           R=RR ;
00347           K0=KK   ;
00348           KW0=KWW ;
00349           MX=K0+L+R /* -1 */ ;
00350           
00351           /*********************/
00352           /*                   */
00353           /*      HJ / GJ      */
00354           READLINE(file,line) ;
00355           ReadInnerEntries(file,HJ.a,&HJ.U,&HJ.O) ;
00356           for (i=0;i<NBC;i++) {
00357             HJ.IsWT=0 ;
00358             READLINE(file,line) ;
00359             HJ.AL[i].Read(file) ;
00360           }
00361           
00362           READLINE(file,line) ;
00363           ReadInnerEntries(file,GtJ.a,&GtJ.U,&GtJ.O) ;
00364           for (i=0;i<NBC;i++) {
00365             GtJ.IsWT=1 ;
00366             READLINE(file,line) ;
00367             GtJ.AL[i].Read(file) ;
00368           }
00369           
00370           /************************/
00371           /*  Lifting Matrizen    */
00372           /*                      */
00373           READLINE(file,line) ;
00374           ReadInnerEntries(file,Q.a,&Q.U,&Q.O) ;
00375           for (i=0;i<NBC;i++) {
00376             READLINE(file,line) ;
00377             Q.QL[i].Read(file)  ;
00378             READLINE(file,line) ;
00379             Q.PL[i].Read(file)  ;
00380           }
00381           
00382           /********************************/
00383           /*     Projektions-Matrizen     */
00384           /*                              */
00385           for (i=0;i<NBC-1;i++) {
00386             READLINE(file,line) ;
00387             OP0[i].Read(file) ;
00388           }         
00389           
00390           /********************************/
00391           /*       Masse-Matrizen         */
00392           /*                              */
00393           READLINE(file,line) ;
00394           ReadInnerEntries(file,MassMatrix.a,&MassMatrix.U,&MassMatrix.O) ;
00395           MassMatrix.power=-1.0 ;
00396           
00397           /********************************/
00398           /*    Steifigkeits-Matrizen     */
00399           /*                              */
00400           READLINE(file,line) ;
00401           ReadInnerEntries(file,StiffnessMatrix.a,&StiffnessMatrix.U,&StiffnessMatrix.O) ;
00402           StiffnessMatrix.power=1.0 ;
00403           
00404           /********************************/
00405           /*    RandIntegral-Matrizen     */
00406           /*                              */
00407           READLINE(file,line)   ;
00408           IntegralsL.Read(file) ;
00409           
00410           /********************************/
00411           /*     Operator-Matrizen        */
00412           /*                              */
00413           
00414           //int cnt=0 ;
00415           for (op=0;op<2;op++)
00416             for (iswT=0;iswT<2;iswT++)
00417               for (iswF=0;iswF<2;iswF++) {
00418                 Operators[op][iswT][iswF].IsWT =iswT ;
00419                 Operators[op][iswT][iswF].IsWF =iswF ;
00420                 Operators[op][iswT][iswF].power=op+1 ;
00421                 
00422                 READLINE(file,line) ;
00423                 ReadInnerEntries(file,Operators[op][iswT][iswF].a,&Operators[op][iswT][iswF].U,&Operators[op][iswT][iswF].O) ;
00424                 for (bcT=0;bcT<NBC;bcT++)
00425                   for (bcF=0;bcF<NBC;bcF++) {
00426                     READLINE(file,line) ;
00427                     
00428                     Operators[op][iswT][iswF].AL [bcT][bcF].Read(file) ;
00429                     Operators[op][iswT][iswF].ALT[bcT][bcF].Read(file) ;
00430                   }
00431               }
00432           
00433           /******************************/
00434           /*   FD-Operator-Matrizen:    */
00435           /*                            */
00436           for (op=0;op<1;op++) {
00437             FDOperators[op].power=1 ;
00438             FDOperators[op].IsWT =FDOperators[op].IsWF=false ;
00439             
00440             READLINE(file,line) ;
00441             ReadInnerEntries(file,FDOperators[op].a,&FDOperators[op].U,&FDOperators[op].O) ;
00442             for (bcT=0;bcT<NBC;bcT++)
00443               for (bcF=0;bcF<NBC;bcF++) {
00444                 READLINE(file,line) ;
00445                 FDOperators[op].AL[bcT][bcF].Read(file) ;
00446               }
00447           }
00448           
00449           
00450         } else { /* right */
00451           
00452           K1 =KK  ;
00453           KW1=KWW ;
00454           MX =(MX+K1)/2 ;
00455           
00456           /*********************/
00457           /*                   */
00458           /*      HJ / GJ      */
00459           for (i=0;i<NBC;i++) {
00460             READLINE(file,line) ;
00461             HJ.AR[i].Read(file) ;
00462           }
00463           
00464           for (i=0;i<NBC;i++) {
00465             READLINE(file,line) ;
00466             GtJ.AR[i].Read(file) ;
00467           }
00468           
00469           /************************/
00470           /*                      */
00471           /* Lifting Matrizen     */
00472           for (i=0;i<NBC;i++) {
00473             READLINE(file,line) ;
00474             Q.QR[i].Read(file)  ;
00475             READLINE(file,line) ;
00476             Q.PR[i].Read(file)  ;
00477           }
00478           
00479           /********************************/
00480           /*     Projektions-Matrizen     */
00481           /*                              */
00482           for (i=0;i<NBC-1;i++) {
00483             READLINE(file,line) ;
00484             OP1[i].Read(file) ;
00485           }         
00486           /********************************/
00487           /*     Operator-Matrizen        */
00488           /*                              */
00489           for (op=0;op<2;op++)
00490             for (iswT=0;iswT<2;iswT++)
00491               for (iswF=0;iswF<2;iswF++) {
00492                 
00493                 for (bcT=0;bcT<NBC;bcT++)
00494                   for (bcF=0;bcF<NBC;bcF++) {
00495                     READLINE(file,line) ;
00496                     Operators[op][iswT][iswF].AR [bcT][bcF].Read(file) ;
00497                     Operators[op][iswT][iswF].ART[bcT][bcF].Read(file) ;
00498                   }
00499                 
00500               } 
00501           
00502           /****************************/
00503           /*   FD-Operator-Matrizen   */
00504           /*                          */
00505           for (op=0;op<1;op++) {
00506             for (bcT=0;bcT<NBC;bcT++)
00507               for (bcF=0;bcF<NBC;bcF++) {
00508                 READLINE(file,line) ;     
00509                 FDOperators[op].AR[bcT][bcF].Read(file) ;
00510               }
00511           }
00512           
00513           
00514           /**********************************/
00515           /* restlichen Transform-Matrizen  */
00516           /* aufbauen                       */
00517           /*                                */
00518           
00519           HtJ.U=HtJ.O=0 ;
00520           HtJ.a[0]   =1.0 ;
00521           HtJ.IsWT   =0   ;
00522           
00523           GJ.U=-1 ;
00524           GJ.O=1  ;
00525           GJ.a[0]=1.0 ;
00526           GJ.IsWT=1   ;
00527           
00528           // Sonderbehandlung fuer RandWavelet bei homogene Neumann-RB
00529           GJ.AL[2].Init(1,3) ;
00530           GJ.AR[2].Init(1,3) ;
00531           GJ.AL[2].a[0][0]=GJ.AL[2].a[0][1]=GJ.AR[2].a[0][1]=GJ.AR[2].a[0][2]=0 ;
00532           GJ.AL[2].a[0][2]=GJ.AR[2].a[0][0]=1 ;
00533           
00534           HtJ.AL[2].Init(2,3) ;
00535           HtJ.AR[2].Init(2,3) ;
00536           HtJ.AL[2].a[0][0]=1 ;
00537           HtJ.AL[2].a[0][1]=0 ;
00538           HtJ.AL[2].a[0][2]=0 ;
00539           HtJ.AL[2].a[1][0]=0 ;
00540           HtJ.AL[2].a[1][1]=0 ;
00541           HtJ.AL[2].a[1][2]=0 ;
00542           
00543           HtJ.AR[2].a[1][0]=0 ;
00544           HtJ.AR[2].a[1][1]=0 ;
00545           HtJ.AR[2].a[1][2]=1 ;
00546           HtJ.AR[2].a[0][0]=0 ;
00547           HtJ.AR[2].a[0][1]=0 ;
00548           HtJ.AR[2].a[0][2]=0 ;
00549           
00550         } 
00551       } /* end if typ/quad ...*/
00552     } /* end if found begin */
00553   } while ((strncmp(&line[0],"stop",4)!=0)&&(found==0)) ;
00554   
00555   if (found==0) {
00556     L=R=N=0 ;
00557   }
00558   
00559   MX1=MX  ;
00560   
00561   Level0=0  ;
00562   i=2*MX1-1 ;
00563   while ((i=(i>>1))>0) { Level0++ ;}
00564   
00565   if (N==2) { Level0=2; }// required for FD24 and such stuff
00566   if (N==4) { Level0=3; }// required for FD24 and such stuff
00567   
00568   fclose(file) ;
00569   
00570   if (left) ReadFDMatrices() ;
00571   
00572   return(found) ;
00573 }
00574 
00575 
00576 //
00577 // FD Matrices
00578 //
00579 
00580 void Wavelets::ReadFDMatrices() {
00581   int  fd,bcF,bcT ;
00582   
00583   FILE *file ;
00584   char name[1000],line[1000] ;
00585   sprintf(name,"%s/MATLAB/Data/FD.dat",_WAVELET_ROOT_) ;
00586   file=fopen(name,"r") ; 
00587   if (file == NULL) {std::cout<<"FD-file not found: "<<name<<'\n';} 
00588   
00589   for (fd=1; fd<MAXINT; fd++) {
00590     // Read comment 
00591     READLINE(file,line) ;
00592     if (line[0]==0) break ;
00593     sscanf(line,"%s %le",name, &FDOperators[fd].power ) ;
00594 
00595     // Read c,U,O   
00596     ReadInnerEntries(file,FDOperators[fd].a,&FDOperators[fd].U,&FDOperators[fd].O) ;
00597     
00598     // Read for all BC the boundary modifications
00599     for (bcF=-1; bcF<=-1; bcF++) 
00600       for (bcT=-1; bcT<=-1; bcT++) {
00601         FDOperators[fd].AL[bcT+1][bcF+1].Read(file) ;
00602         FDOperators[fd].AR[bcT+1][bcF+1].Read(file) ;
00603       }
00604   }
00605   fclose(file) ; 
00606 }
00607 
00608 
00609 int Wavelets::FDOp(int op) {
00610  switch (op) {
00611   case FD12:    return  1;
00612   case FD22:    return  2;
00613   case FD24:    return  3;
00614   case FD40:    return  4;
00615   case FD60:    return  5;
00616   case FD80:    return  6;
00617   case FD14:    return  7;
00618   case FD16:    return  8;
00619   case DIVFD4:  return  9;
00620   case GRADFD4: return 10;
00621   case FD11:    return 11;
00622  default: assert(0) ;
00623  }
00624  return -MAXINT ;
00625 } 
00626 
00627 int Wavelets::TypeID(char *typ) {
00628   if (strcmp(typ,"Interpolet")==0) return 0 ;
00629   if (strcmp(typ,"Daubechies"   )==0) return 1 ;
00630   return -1 ;
00631 }
00632 
00633 
00634 //
00635 double Wavelets::WaveletIntegral(int l, int s, int *BC) {
00636   assert( InterpoletFlag) ;
00637   assert(!LiftingFlag   ) ;
00638   
00639   // special cases
00640   if ((l>Level0) && (s==0            ) && (BC[0]==1)) return ScalingFunctionIntegral(l,       2 ,BC) ;
00641   if ((l>Level0) && (s==(1<<(l-1))-1 ) && (BC[1]==1)) return ScalingFunctionIntegral(l,(1<<l)-2 ,BC) ;
00642   if (l>Level0)                                       return ScalingFunctionIntegral(l,2*s+1,BC) ;
00643   return ScalingFunctionIntegral(l,s,BC) ; // l==Level0 
00644 }
00645   
00646 double Wavelets::ScalingFunctionIntegral(int l, int s, int *BC) {
00647   int n=1<<l ;
00648   if (BC[1]==PERIODIC) return 1./n ;
00649   // boundaries
00650   if (s<N  ) { int in[2]={BC[0]+1, s  }; return IntegralsL[in]/n ; }
00651   if (s>n-N) { int in[2]={BC[1]+1, n-s}; return IntegralsL[in]/n ; }   
00652     
00653   return 1./n ;
00654 }
00655 
00656 
00657 void Wavelets::Print(int allflag) {
00658   std::cout<< "Wavelets ("<<type<<") Interpolets "<< ((InterpoletFlag) ? "yes" : "no") << this <<'\n' ;
00659   std::cout<< "N="<<N<<" L="<<L<<" R="<<R<<" K0="<<K0<<" K1="<<K1<<" KW0="<<KW0<<" KW1="<<KW1<<" Level0="<<Level0<<" NBC="<<NBC<<'\n';
00660   
00661   if (allflag & 1) {
00662     std::cout<< "primal low  pass filters H^l \n"; HJ.Print() ;
00663     std::cout<< "primal high pass filters G^l \n"; GJ.Print() ;
00664     std::cout<< "dual   low  pass filters H^tl\n"; HtJ.Print() ;
00665     std::cout<< "dual   high pass filters G^tl\n"; GtJ.Print() ;
00666   }
00667  
00668   if (allflag & 2) {
00669     ;
00670   } 
00671  
00672 }
Generated at Mon Aug 19 10:02:32 2002 for AWFD by