/************************************************************************
 * The following is a collection of wavelet-analysis related routines
 * @ MAR-1994 EK (original code in fortran)
 * @ JUN-2001 - ported to C
 ************************************************************************/
#include <stdio.h>
#include "mywv2lib.h"
#include <math.h>
#include <stdlib.h>

int main(int   argc,
         char *argv[])
{
   double *array1;
   double *array2;
   double *array3;
   double temp;
   double mse;
   double rms;
   register int i;
   int filter;
   int size   = 16384;

   if(argc>1 && atoi(argv[1])>0)
   {
	   size = atoi(argv[1]);
   }
   array1 = (double *) malloc(size*sizeof(double));
   if(array1==NULL)
   {
	   printf("Not enough memory (%d bytes) for %s - exiting\n",
		      size*sizeof(double),"array1");
	   exit(-1);
   }
   array2 = (double *) malloc(size*sizeof(double));
   if(array2==NULL)
   {
	   free(array1);
	   printf("Not enough memory (%d bytes) for %s - exiting\n",
		      size*sizeof(double),"array2");
	   exit(-1);
   }
   array3 = (double *) malloc(size*sizeof(double));
   if(array3==NULL)
   {
	   free(array2);
	   free(array1);
	   printf("Not enough memory (%d bytes) for %s - exiting\n",
		      size*sizeof(double),"array3");
	   exit(-1);
   }
   for(i=0;i<size;i++)
   {
      array3[i] = (double) rand()/(double) RAND_MAX;
   }
   for(filter=1;filter<=4;filter++)
   {
	   for(i=0;i<size;i++)
	   {
		  array1[i] = array3[i];
	   }
	   wavelet(array1,size,filter,1);
	   for(i=0;i<size;i++)
	   {
		  array2[i] = array1[i];
	   }
	   wavelet(array1,size,filter,-1);
	   mse = 0.0;
	   for(i=0;i<size;i++)
	   {
		  temp = array1[i] - array3[i];
		  mse += temp*temp;
	   }
	   mse/=(double) size;
	   rms = sqrt(mse);
	   printf("Filter = %d, RMS = %e\n",filter,rms);
   }
   free(array3);
   free(array2);
   free(array1);
   return(-1);
}