daubtst.c

Dwt离散db小波 daubechi.c daubtst.c

/**********************************************************************/
/*	tested daubechies.c					      */
/**********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int decomp(double *sig,int length_sig,int length_filt,double *h,double *g);
int reconst(double *spec,int length_spec,int length_filt,double *h,double *g);

#define N 2048

void main(int argc, char *argv[])
{
 FILE *inout, *out;
 double g[8],h[8],sig[N];
 int order,filtlength,i,j,dir,n;
 char infilename[20],outfilename[20],numva[3];

/*********************************************************************/
/*	analysis of command line				     */
/*********************************************************************/

if(argc<5)	/*	missing options				*/
  {fprintf(stderr,"Usage:\n");
   fprintf(stderr,"daub -Iinputfile -Ooutputfile -D{1/-1} -L{1/2/3/4}\n");
   fprintf(stderr,"\tinputfile\tname of the inputfile\n");
   fprintf(stderr,"\toutputfile\tname of the outputfile\n");
   fprintf(stderr,"\tD\t\tdirection :\t1=wavelet transform\n");
   fprintf(stderr,"\t\t\t \t\t-1=inverse wavelet transform\n");
   fprintf(stderr,"\tL\t\torder of wavelet 1,2,3,4\n");
   exit(1);
   }  
 

 for(i=1;i<argc;i++)
    switch(argv[i][1])
       {case 'I' : for(j=2;j<strlen(argv[i]);j++) 
			infilename[j-2]=argv[i][j];
                    infilename[j-2]='\0';
                    break;
        case 'O' :for(j=2;j<strlen(argv[i]);j++) 
			outfilename[j-2]=argv[i][j];
                     outfilename[j-2]='\0';
                    break;
        case 'D' : for(j=2;j<strlen(argv[i]);j++) numva[j-2]=argv[i][j];
                    numva[j-2]='\0';
                   dir=atoi(numva);
                   strcpy(numva,"");;
                   break;
        case 'L' :for(j=2;j<strlen(argv[i]);j++) numva[j-2]=argv[i][j];
                    numva[j-2]='\0';
                   order=atoi(numva);
                   strcpy(numva,"");
                   break;
       }

/*****************************************************************************/
/*		Data input						     */
/*****************************************************************************/
 
inout=fopen(infilename,"r");
 if(inout==NULL) {fprintf(stderr,"Signal file not found\n");exit(1);}
 
    i=0; 
    while(fscanf(inout,"%le",sig+i)!=EOF) i++;
    n=i;

/****************************************************************************/
/*	direct or inverse wavelet transformation			    */
/****************************************************************************/

if(daub_init(order,&filtlength,h,g))
            {fprintf(stderr,"Filter of order %d not implemented\n");exit(1);}

switch(dir)
  {
   case 1: 

             if(decomp(sig,n,filtlength,h,g)) 
                {fprintf(stderr,"not enough memory\n");exit(1);}
            break;


    case -1: if(reconst(sig,n,filtlength,h,g))
                {fprintf(stderr,"not enough memory\n");exit(1);}
             break;
    }

 

 
   out=fopen(outfilename,"w");
   if(out==NULL) {fprintf(stderr,"error while opening output file\n");
                       exit(1);
		      }
		       
   for(i=0;i<n;i++) fprintf(out,"%e\n",sig[i]);
   fclose(out);
}