htswanalysis/MACS/lib/gsl/gsl-1.11/fft/test_real_source.c

   1 /* fft/test_real.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include "bitreverse.h"
  21 #include "signals.h"
  22 #include "compare.h"
  23
  24 void FUNCTION(test_real,func) (size_t stride, size_t n);
  25 void FUNCTION(test_real,bitreverse_order) (size_t stride, size_t n);
  26 void FUNCTION(test_real,radix2) (size_t stride, size_t n);
  27
  28 void FUNCTION(test_real,func) (size_t stride, size_t n)
  29 {
  30   size_t i ;
  31   int status ;
  32
  33   TYPE(gsl_fft_real_wavetable) * rw ;
  34   TYPE(gsl_fft_halfcomplex_wavetable) * hcw ;
  35   TYPE(gsl_fft_real_workspace) * rwork ;
  36
  37   BASE * real_data = (BASE *) malloc (n * stride * sizeof (BASE));
  38   BASE * complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
  39   BASE * complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
  40   BASE * fft_complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
  41
  42   for (i = 0 ; i < n * stride ; i++)
  43     {
  44       real_data[i] = (BASE)i ;
  45     }
  46
  47   for (i = 0 ; i < 2 * n * stride ; i++)
  48     {
  49       complex_data[i] = (BASE)(i + 1000.0) ;
  50       complex_tmp[i] = (BASE)(i + 2000.0) ;
  51       fft_complex_data[i] = (BASE)(i + 3000.0) ;
  52     }
  53
  54   gsl_set_error_handler (NULL); /* abort on any errors */
  55
  56   /* mixed radix real fft */
  57
  58   rw = FUNCTION(gsl_fft_real_wavetable,alloc) (n);
  59   gsl_test (rw == 0, NAME(gsl_fft_real_wavetable)
  60             "_alloc, n = %d, stride = %d", n, stride);
  61
  62   rwork = FUNCTION(gsl_fft_real_workspace,alloc) (n);
  63   gsl_test (rwork == 0, NAME(gsl_fft_real_workspace)
  64             "_alloc, n = %d", n);
  65
  66   FUNCTION(fft_signal,real_noise) (n, stride, complex_data, fft_complex_data);
  67   memcpy (complex_tmp, complex_data, 2 * n * stride * sizeof (BASE));
  68
  69   for (i = 0; i < n; i++)
  70     {
  71       real_data[i*stride] = REAL(complex_data,stride,i);
  72     }
  73
  74   FUNCTION(gsl_fft_real,transform) (real_data, stride, n, rw, rwork);
  75   FUNCTION(gsl_fft_halfcomplex,unpack) (real_data, complex_data, stride, n);
  76
  77   status = FUNCTION(compare_complex,results) ("dft", fft_complex_data,
  78                                               "fft of noise", complex_data,
  79                                               stride, n, 1e6);
  80   gsl_test (status, NAME(gsl_fft_real)
  81             " with signal_real_noise, n = %d, stride = %d", n, stride);
  82
  83   /* compute the inverse fft */
  84
  85   hcw = FUNCTION(gsl_fft_halfcomplex_wavetable,alloc) (n);
  86   gsl_test (hcw == 0, NAME(gsl_fft_halfcomplex_wavetable)
  87             "_alloc, n = %d, stride = %d", n, stride);
  88
  89   status = FUNCTION(gsl_fft_halfcomplex,transform) (real_data, stride, n, hcw, rwork);
  90
  91   for (i = 0; i < n; i++)
  92     {
  93       real_data[i*stride] /= n;
  94     }
  95
  96   FUNCTION(gsl_fft_real,unpack) (real_data, complex_data, stride, n);
  97
  98   status = FUNCTION(compare_complex,results) ("orig", complex_tmp,
  99                                               "fft inverse", complex_data,
 100                                               stride, n, 1e6);
 101
 102   gsl_test (status, NAME(gsl_fft_halfcomplex)
 103             " with data from signal_noise, n = %d, stride = %d", n, stride);
 104
 105   FUNCTION(gsl_fft_real_workspace,free) (rwork);
 106   FUNCTION(gsl_fft_real_wavetable,free) (rw);
 107   FUNCTION(gsl_fft_halfcomplex_wavetable,free) (hcw);
 108
 109   free(real_data) ;
 110   free(complex_data) ;
 111   free(complex_tmp) ;
 112   free(fft_complex_data) ;
 113 }
 114
 115
 116 void
 117 FUNCTION(test_real,bitreverse_order) (size_t stride, size_t n)
 118 {
 119   int status ;
 120   size_t logn, i ;
 121
 122   BASE * tmp = (BASE *) malloc (n * stride * sizeof (BASE));
 123   BASE * data = (BASE *) malloc (n * stride * sizeof (BASE));
 124   BASE * reversed_data = (BASE *) malloc (n * stride * sizeof (BASE));
 125
 126   for (i = 0; i <  stride * n; i++)
 127     {
 128       data[i] = (BASE)i ;
 129     }
 130
 131   memcpy (tmp, data, n * stride * sizeof(BASE)) ;
 132
 133   logn = 0 ;  while (n > (1U <<logn)) {logn++;} ;
 134
 135   /* do a naive bit reversal as a baseline for testing the other routines */
 136
 137   for (i = 0; i < n; i++)
 138     {
 139       size_t i_tmp = i ;
 140       size_t j = 0 ;
 141       size_t bit ;
 142
 143       for (bit = 0; bit < logn; bit++)
 144         {
 145           j <<= 1;              /* reverse shift i into j */
 146           j |= i_tmp & 1;
 147           i_tmp >>= 1;
 148         }
 149
 150       reversed_data[j*stride] = data[i*stride] ;
 151     }
 152
 153   FUNCTION(fft_real,bitreverse_order) (data, stride, n, logn);
 154
 155   status = FUNCTION(compare_real,results) ("naive bit reverse",
 156                                            reversed_data,
 157                                            "fft_complex_bitreverse_order",
 158                                            data,
 159                                            stride, n, 1e6);
 160
 161   gsl_test (status, NAME(fft_real) "_bitreverse_order, n = %d", n);
 162
 163   free (reversed_data) ;
 164   free (data) ;
 165   free (tmp) ;
 166 }
 167
 168
 169 void FUNCTION(test_real,radix2) (size_t stride, size_t n)
 170 {
 171   size_t i ;
 172   int status ;
 173
 174   BASE * real_data = (BASE *) malloc (n * stride * sizeof (BASE));
 175   BASE * complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 176   BASE * complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 177   BASE * fft_complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 178
 179   for (i = 0 ; i < n * stride ; i++)
 180     {
 181       real_data[i] = (BASE)i ;
 182     }
 183
 184   for (i = 0 ; i < 2 * n * stride ; i++)
 185     {
 186       complex_data[i] = (BASE)(i + 1000.0) ;
 187       complex_tmp[i] = (BASE)(i + 2000.0) ;
 188       fft_complex_data[i] = (BASE)(i + 3000.0) ;
 189     }
 190
 191   gsl_set_error_handler (NULL); /* abort on any errors */
 192
 193   /* radix-2 real fft */
 194
 195   FUNCTION(fft_signal,real_noise) (n, stride, complex_data, fft_complex_data);
 196   memcpy (complex_tmp, complex_data, 2 * n * stride * sizeof (BASE));
 197
 198   for (i = 0; i < n; i++)
 199     {
 200       real_data[i*stride] = REAL(complex_data,stride,i);
 201     }
 202
 203   FUNCTION(gsl_fft_real,radix2_transform) (real_data, stride, n);
 204   FUNCTION(gsl_fft_halfcomplex,radix2_unpack) (real_data, complex_data, stride, n);
 205
 206   status = FUNCTION(compare_complex,results) ("dft", fft_complex_data,
 207                                               "fft of noise", complex_data,
 208                                               stride, n, 1e6);
 209   gsl_test (status, NAME(gsl_fft_real)
 210             "_radix2 with signal_real_noise, n = %d, stride = %d", n, stride);
 211
 212   /* compute the inverse fft */
 213
 214   status = FUNCTION(gsl_fft_halfcomplex,radix2_transform) (real_data, stride, n);
 215
 216   for (i = 0; i < n; i++)
 217     {
 218       real_data[i*stride] /= n;
 219     }
 220
 221   FUNCTION(gsl_fft_real,unpack) (real_data, complex_data, stride, n);
 222
 223   status = FUNCTION(compare_complex,results) ("orig", complex_tmp,
 224                                               "fft inverse", complex_data,
 225                                               stride, n, 1e6);
 226
 227   gsl_test (status, NAME(gsl_fft_halfcomplex)
 228             "_radix2 with data from signal_noise, n = %d, stride = %d", n, stride);
 229
 230
 231   free(real_data) ;
 232   free(complex_data) ;
 233   free(complex_tmp) ;
 234   free(fft_complex_data) ;
 235 }