31.4 Examples

The following program demonstrates the use of the one-dimensional wavelet transform functions. It computes an approximation to an input signal (of length 256) using the 20 largest components of the wavelet transform, while setting the others to zero.

     #include <stdio.h>
     #include <math.h>
     #include <gsl/gsl_sort.h>
     #include <gsl/gsl_wavelet.h>
     
     int
     main (int argc, char **argv)
     {
       int i, n = 256, nc = 20;
       double *data = malloc (n * sizeof (double));
       double *abscoeff = malloc (n * sizeof (double));
       size_t *p = malloc (n * sizeof (size_t));
     
       FILE * f;
       gsl_wavelet *w;
       gsl_wavelet_workspace *work;
     
       w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
       work = gsl_wavelet_workspace_alloc (n);
     
       f = fopen (argv[1], "r");
       for (i = 0; i < n; i++)
         {
           fscanf (f, "%lg", &data[i]);
         }
       fclose (f);
     
       gsl_wavelet_transform_forward (w, data, 1, n, work);
     
       for (i = 0; i < n; i++)
         {
           abscoeff[i] = fabs (data[i]);
         }
       
       gsl_sort_index (p, abscoeff, 1, n);
       
       for (i = 0; (i + nc) < n; i++)
         data[p[i]] = 0;
       
       gsl_wavelet_transform_inverse (w, data, 1, n, work);
       
       for (i = 0; i < n; i++)
         {
           printf ("%g\n", data[i]);
         }
       
       gsl_wavelet_free (w);
       gsl_wavelet_workspace_free (work);
     
       free (data);
       free (abscoeff);
       free (p);
       return 0;
     }

The output can be used with the gnu plotutils graph program,

     $ ./a.out ecg.dat > dwt.dat
     $ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat > dwt.ps