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

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <gsl/gsl_errno.h>
   4 #include <gsl/gsl_fft_complex.h>
   5
   6 #define REAL(z,i) ((z)[2*(i)])
   7 #define IMAG(z,i) ((z)[2*(i)+1])
   8
   9 int
  10 main (void)
  11 {
  12   int i;
  13   const int n = 630;
  14   double data[2*n];
  15
  16   gsl_fft_complex_wavetable * wavetable;
  17   gsl_fft_complex_workspace * workspace;
  18
  19   for (i = 0; i < n; i++)
  20     {
  21       REAL(data,i) = 0.0;
  22       IMAG(data,i) = 0.0;
  23     }
  24
  25   data[0] = 1.0;
  26
  27   for (i = 1; i <= 10; i++)
  28     {
  29       REAL(data,i) = REAL(data,n-i) = 1.0;
  30     }
  31
  32   for (i = 0; i < n; i++)
  33     {
  34       printf ("%d: %e %e\n", i, REAL(data,i),
  35                                 IMAG(data,i));
  36     }
  37   printf ("\n");
  38
  39   wavetable = gsl_fft_complex_wavetable_alloc (n);
  40   workspace = gsl_fft_complex_workspace_alloc (n);
  41
  42   for (i = 0; i < wavetable->nf; i++)
  43     {
  44        printf ("# factor %d: %d\n", i,
  45                wavetable->factor[i]);
  46     }
  47
  48   gsl_fft_complex_forward (data, 1, n,
  49                            wavetable, workspace);
  50
  51   for (i = 0; i < n; i++)
  52     {
  53       printf ("%d: %e %e\n", i, REAL(data,i),
  54                                 IMAG(data,i));
  55     }
  56
  57   gsl_fft_complex_wavetable_free (wavetable);
  58   gsl_fft_complex_workspace_free (workspace);
  59   return 0;
  60 }