htswanalysis/MACS/lib/gsl/gsl-1.11/doc/examples/dwt.c

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <gsl/gsl_sort.h>
   4 #include <gsl/gsl_wavelet.h>
   5
   6 int
   7 main (int argc, char **argv)
   8 {
   9   int i, n = 256, nc = 20;
  10   double *data = malloc (n * sizeof (double));
  11   double *abscoeff = malloc (n * sizeof (double));
  12   size_t *p = malloc (n * sizeof (size_t));
  13
  14   FILE * f;
  15   gsl_wavelet *w;
  16   gsl_wavelet_workspace *work;
  17
  18   w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
  19   work = gsl_wavelet_workspace_alloc (n);
  20
  21   f = fopen (argv[1], "r");
  22   for (i = 0; i < n; i++)
  23     {
  24       fscanf (f, "%lg", &data[i]);
  25     }
  26   fclose (f);
  27
  28   gsl_wavelet_transform_forward (w, data, 1, n, work);
  29
  30   for (i = 0; i < n; i++)
  31     {
  32       abscoeff[i] = fabs (data[i]);
  33     }
  34
  35   gsl_sort_index (p, abscoeff, 1, n);
  36
  37   for (i = 0; (i + nc) < n; i++)
  38     data[p[i]] = 0;
  39
  40   gsl_wavelet_transform_inverse (w, data, 1, n, work);
  41
  42   for (i = 0; i < n; i++)
  43     {
  44       printf ("%g\n", data[i]);
  45     }
  46
  47   gsl_wavelet_free (w);
  48   gsl_wavelet_workspace_free (work);
  49
  50   free (data);
  51   free (abscoeff);
  52   free (p);
  53   return 0;
  54 }