#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_wavelet.h>

int
main (int argc, char **argv)
{
  (void)(argc); /* avoid unused parameter warning */
  int i, n = 256, nc = 20;
  double *orig_data = malloc (n * sizeof (double));
  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", &orig_data[i]);
      data[i] = orig_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 %g\n", orig_data[i], data[i]);
    }
  
  gsl_wavelet_free (w);
  gsl_wavelet_workspace_free (work);

  free (data);
  free (orig_data);
  free (abscoeff);
  free (p);
  return 0;
}