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

   1 /* fft/test_complex.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_complex,func) (size_t stride, size_t n);
  25 int FUNCTION(test,offset) (const BASE data[], size_t stride,
  26                            size_t n, size_t offset);
  27 void FUNCTION(test_complex,bitreverse_order) (size_t stride, size_t n) ;
  28 void FUNCTION(test_complex,radix2) (size_t stride, size_t n);
  29
  30 int FUNCTION(test,offset) (const BASE data[], size_t stride,
  31                            size_t n, size_t offset)
  32 {
  33   int status = 0 ;
  34   size_t i, j, k = 0 ;
  35
  36   for (i = 0; i < n; i++)
  37     {
  38       k += 2 ;
  39
  40       for (j = 1; j < stride; j++)
  41         {
  42           status |= data[k] != k + offset ;
  43           k++ ;
  44           status |= data[k] != k + offset ;
  45           k++ ;
  46         }
  47     }
  48   return status ;
  49 }
  50
  51
  52 void FUNCTION(test_complex,func) (size_t stride, size_t n)
  53 {
  54   size_t i ;
  55   int status ;
  56
  57   TYPE(gsl_fft_complex_wavetable) * cw ;
  58   TYPE(gsl_fft_complex_workspace) * cwork ;
  59
  60   BASE * complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
  61   BASE * complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
  62   BASE * fft_complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
  63   BASE * fft_complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
  64
  65   for (i = 0 ; i < 2 * n * stride ; i++)
  66     {
  67       complex_data[i] = (BASE)i ;
  68       complex_tmp[i] = (BASE)(i + 1000.0) ;
  69       fft_complex_data[i] = (BASE)(i + 2000.0) ;
  70       fft_complex_tmp[i] = (BASE)(i + 3000.0) ;
  71     }
  72
  73   gsl_set_error_handler (NULL); /* abort on any errors */
  74
  75   /* Test allocation */
  76
  77   {
  78     cw = FUNCTION(gsl_fft_complex_wavetable,alloc) (n);
  79     gsl_test (cw == 0, NAME(gsl_fft_complex_wavetable)
  80               "_alloc, n = %d, stride = %d", n, stride);
  81   }
  82
  83   {
  84     cwork = FUNCTION(gsl_fft_complex_workspace,alloc) (n);
  85     gsl_test (cwork == 0, NAME(gsl_fft_complex_workspace)
  86               "_alloc, n = %d", n);
  87   }
  88
  89
  90   /* Test mixed radix fft with noise */
  91
  92   {
  93     FUNCTION(fft_signal,complex_noise) (n, stride, complex_data, fft_complex_data);
  94     for (i = 0 ; i < n ; i++)
  95       {
  96         REAL(complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
  97         IMAG(complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
  98       }
  99
 100     FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
 101
 102     for (i = 0 ; i < n ; i++)
 103       {
 104         REAL(fft_complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
 105         IMAG(fft_complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
 106       }
 107
 108     status = FUNCTION(compare_complex,results) ("dft", fft_complex_data,
 109                                                 "fft of noise", complex_data,
 110                                                 stride, n, 1e6);
 111     gsl_test (status, NAME(gsl_fft_complex)
 112               "_forward with signal_noise, n = %d, stride = %d",  n, stride);
 113
 114     if (stride > 1)
 115       {
 116         status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
 117
 118         gsl_test (status, NAME(gsl_fft_complex)
 119                   "_forward avoids unstrided data, n = %d, stride = %d",
 120                   n, stride);
 121       }
 122   }
 123
 124   /* Test the inverse fft */
 125
 126   {
 127     status = FUNCTION(gsl_fft_complex,inverse) (complex_data, stride, n, cw, cwork);
 128     status = FUNCTION(compare_complex,results) ("orig", complex_tmp,
 129                                                 "fft inverse", complex_data,
 130                                                 stride, n, 1e6);
 131     gsl_test (status, NAME(gsl_fft_complex)
 132               "_inverse with signal_noise, n = %d, stride = %d", n, stride);
 133
 134     if (stride > 1)
 135       {
 136         status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
 137
 138         gsl_test (status, NAME(gsl_fft_complex)
 139                   "_inverse other data untouched, n = %d, stride = %d",
 140                   n, stride);
 141       }
 142
 143   }
 144
 145   /* Test the backward fft */
 146
 147   {
 148     status = FUNCTION(gsl_fft_complex,backward) (fft_complex_tmp, stride, n, cw, cwork);
 149
 150     for (i = 0; i < n; i++)
 151       {
 152         REAL(complex_tmp,stride,i) *= n;
 153         IMAG(complex_tmp,stride,i) *= n;
 154       }
 155     status = FUNCTION(compare_complex,results) ("orig",
 156                                                 complex_tmp,
 157                                                 "fft backward",
 158                                                 fft_complex_tmp,
 159                                                 stride, n, 1e6);
 160
 161     gsl_test (status, NAME(gsl_fft_complex)
 162               "_backward with signal_noise, n = %d, stride = %d", n, stride);
 163
 164     if (stride > 1)
 165       {
 166         status = FUNCTION(test, offset) (fft_complex_tmp, stride, n, 3000) ;
 167
 168         gsl_test (status, NAME(gsl_fft_complex)
 169                   "_backward avoids unstrided data, n = %d, stride = %d",
 170                   n, stride);
 171       }
 172
 173   }
 174
 175   /* Test a pulse signal */
 176
 177   {
 178     FUNCTION(fft_signal,complex_pulse) (1, n, stride, 1.0, 0.0, complex_data,
 179                                         fft_complex_data);
 180     FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
 181     status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
 182                                                 "fft of pulse", complex_data,
 183                                                 stride, n, 1e6);
 184     gsl_test (status, NAME(gsl_fft_complex)
 185               "_forward with signal_pulse, n = %d, stride = %d", n, stride);
 186
 187   }
 188
 189
 190   /* Test a constant signal */
 191
 192   {
 193     FUNCTION(fft_signal,complex_constant) (n, stride, 1.0, 0.0, complex_data,
 194                                            fft_complex_data);
 195     FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
 196     status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
 197                                                 "fft of constant",
 198                                                 complex_data,
 199                                                 stride, n, 1e6);
 200     gsl_test (status, NAME(gsl_fft_complex)
 201               "_forward with signal_constant, n = %d, stride = %d", n, stride);
 202   }
 203
 204   /* Test an exponential (cos/sin) signal */
 205
 206   {
 207     status = 0;
 208     for (i = 0; i < n; i++)
 209       {
 210         FUNCTION(fft_signal,complex_exp) ((int)i, n, stride, 1.0, 0.0, complex_data,
 211                                           fft_complex_data);
 212         FUNCTION(gsl_fft_complex,forward) (complex_data, stride, n, cw, cwork);
 213         status |= FUNCTION(compare_complex,results) ("analytic",
 214                                                      fft_complex_data,
 215                                                      "fft of exp",
 216                                                      complex_data,
 217                                                      stride, n, 1e6);
 218       }
 219     gsl_test (status, NAME(gsl_fft_complex)
 220               "_forward with signal_exp, n = %d, stride = %d", n, stride);
 221   }
 222
 223   FUNCTION(gsl_fft_complex_wavetable,free) (cw);
 224   FUNCTION(gsl_fft_complex_workspace,free) (cwork);
 225
 226   free (complex_data);
 227   free (complex_tmp);
 228   free (fft_complex_data);
 229   free (fft_complex_tmp);
 230 }
 231
 232
 233 void
 234 FUNCTION(test_complex,bitreverse_order) (size_t stride, size_t n)
 235 {
 236   int status ;
 237   size_t logn, i ;
 238
 239   BASE * tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 240   BASE * data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 241   BASE * reversed_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 242
 243   for (i = 0; i <  2 * stride * n; i++)
 244     {
 245       data[i] = (BASE)i ;
 246     }
 247
 248   memcpy (tmp, data, 2 * n * stride * sizeof(BASE)) ;
 249
 250   logn = 0 ; while (n > (1U<<logn)) {logn++ ; } ;
 251
 252   /* do a naive bit reversal as a baseline for testing the other routines */
 253
 254   for (i = 0; i < n; i++)
 255     {
 256       size_t i_tmp = i ;
 257       size_t j = 0 ;
 258       size_t bit ;
 259
 260       for (bit = 0; bit < logn; bit++)
 261         {
 262           j <<= 1;              /* reverse shift i into j */
 263           j |= i_tmp & 1;
 264           i_tmp >>= 1;
 265         }
 266
 267       reversed_data[2*j*stride] = data[2*i*stride] ;
 268       reversed_data[2*j*stride+1] = data[2*i*stride+1] ;
 269     }
 270
 271   FUNCTION(fft_complex,bitreverse_order) (data, stride, n, logn);
 272
 273   status = FUNCTION(compare_complex,results) ("naive bit reverse",
 274                                               reversed_data,
 275                                               "fft_complex_bitreverse_order",
 276                                               data,
 277                                               stride, n, 1e6);
 278
 279   gsl_test (status, "fft_complex_bitreverse_order, n = %d", n);
 280
 281   free (reversed_data) ;
 282   free (data) ;
 283   free (tmp) ;
 284 }
 285
 286 void FUNCTION(test_complex,radix2) (size_t stride, size_t n)
 287 {
 288   size_t i ;
 289   int status ;
 290
 291   BASE * complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 292   BASE * complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 293   BASE * fft_complex_data = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 294   BASE * fft_complex_tmp = (BASE *) malloc (2 * n * stride * sizeof (BASE));
 295
 296   for (i = 0 ; i < 2 * n * stride ; i++)
 297     {
 298       complex_data[i] = (BASE)i ;
 299       complex_tmp[i] = (BASE)(i + 1000.0) ;
 300       fft_complex_data[i] = (BASE)(i + 2000.0) ;
 301       fft_complex_tmp[i] = (BASE)(i + 3000.0) ;
 302     }
 303
 304   gsl_set_error_handler (NULL); /* abort on any errors */
 305
 306   /* Test radix-2 fft with noise */
 307
 308   {
 309     FUNCTION(fft_signal,complex_noise) (n, stride, complex_data,
 310                                         fft_complex_data);
 311     for (i = 0 ; i < n ; i++)
 312       {
 313         REAL(complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
 314         IMAG(complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
 315       }
 316
 317     FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
 318
 319     for (i = 0 ; i < n ; i++)
 320       {
 321         REAL(fft_complex_tmp,stride,i) = REAL(complex_data,stride,i) ;
 322         IMAG(fft_complex_tmp,stride,i) = IMAG(complex_data,stride,i) ;
 323       }
 324
 325     status = FUNCTION(compare_complex,results) ("dft", fft_complex_data,
 326                                                 "fft of noise", complex_data,
 327                                                 stride, n, 1e6);
 328     gsl_test (status, NAME(gsl_fft_complex)
 329               "_radix2_forward with signal_noise, n = %d, stride = %d",
 330               n, stride);
 331
 332     if (stride > 1)
 333       {
 334         status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
 335
 336         gsl_test (status, NAME(gsl_fft_complex)
 337                   "_radix2_forward avoids unstrided data, n = %d, stride = %d",
 338                   n, stride);
 339       }
 340   }
 341
 342   /* Test the inverse fft */
 343
 344   {
 345     status = FUNCTION(gsl_fft_complex,radix2_inverse) (complex_data, stride, n);
 346     status = FUNCTION(compare_complex,results) ("orig", complex_tmp,
 347                                                 "fft inverse", complex_data,
 348                                                 stride, n, 1e6);
 349     gsl_test (status, NAME(gsl_fft_complex)
 350               "_radix2_inverse with signal_noise, n = %d, stride = %d", n, stride);
 351
 352     if (stride > 1)
 353       {
 354         status = FUNCTION(test, offset) (complex_data, stride, n, 0) ;
 355
 356         gsl_test (status, NAME(gsl_fft_complex)
 357                   "_radix2_inverse other data untouched, n = %d, stride = %d",
 358                   n, stride);
 359       }
 360
 361   }
 362
 363   /* Test the backward fft */
 364
 365   {
 366     status = FUNCTION(gsl_fft_complex,radix2_backward) (fft_complex_tmp, stride, n);
 367
 368     for (i = 0; i < n; i++)
 369       {
 370         REAL(complex_tmp,stride,i) *= n;
 371         IMAG(complex_tmp,stride,i) *= n;
 372       }
 373     status = FUNCTION(compare_complex,results) ("orig",
 374                                                 complex_tmp,
 375                                                 "fft backward",
 376                                                 fft_complex_tmp,
 377                                                 stride, n, 1e6);
 378
 379     gsl_test (status, NAME(gsl_fft_complex)
 380               "_radix2_backward with signal_noise, n = %d, stride = %d", n, stride);
 381
 382     if (stride > 1)
 383       {
 384         status = FUNCTION(test, offset) (fft_complex_tmp, stride, n, 3000) ;
 385
 386         gsl_test (status, NAME(gsl_fft_complex)
 387                   "_radix2_backward avoids unstrided data, n = %d, stride = %d",
 388                   n, stride);
 389       }
 390
 391   }
 392
 393   /* Test a pulse signal */
 394
 395   {
 396     FUNCTION(fft_signal,complex_pulse) (1, n, stride, 1.0, 0.0, complex_data,
 397                                         fft_complex_data);
 398     FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
 399     status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
 400                                                 "fft of pulse", complex_data,
 401                                                 stride, n, 1e6);
 402     gsl_test (status, NAME(gsl_fft_complex)
 403               "_radix2_forward with signal_pulse, n = %d, stride = %d", n, stride);
 404
 405   }
 406
 407
 408   /* Test a constant signal */
 409
 410   {
 411     FUNCTION(fft_signal,complex_constant) (n, stride, 1.0, 0.0, complex_data,
 412                                            fft_complex_data);
 413     FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
 414     status = FUNCTION(compare_complex,results) ("analytic", fft_complex_data,
 415                                                 "fft of constant",
 416                                                 complex_data,
 417                                                 stride, n, 1e6);
 418     gsl_test (status, NAME(gsl_fft_complex)
 419               "_radix2_forward with signal_constant, n = %d, stride = %d",
 420               n, stride);
 421   }
 422
 423   /* Test an exponential (cos/sin) signal */
 424
 425   {
 426     status = 0;
 427     for (i = 0; i < n; i++)
 428       {
 429         FUNCTION(fft_signal,complex_exp) ((int)i, n, stride, 1.0, 0.0, complex_data,
 430                                           fft_complex_data);
 431         FUNCTION(gsl_fft_complex,radix2_forward) (complex_data, stride, n);
 432         status |= FUNCTION(compare_complex,results) ("analytic",
 433                                                      fft_complex_data,
 434                                                      "fft of exp",
 435                                                      complex_data,
 436                                                      stride, n, 1e6);
 437       }
 438     gsl_test (status, NAME(gsl_fft_complex)
 439               "_radix2_forward with signal_exp, n = %d, stride = %d", n, stride);
 440   }
 441
 442
 443   free (complex_data);
 444   free (complex_tmp);
 445   free (fft_complex_data);
 446   free (fft_complex_tmp);
 447 }
 448