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

   1 /* fft/real_main.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 <config.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_complex.h>
  26 #include <gsl/gsl_fft_real.h>
  27
  28 #include "real_pass.h"
  29
  30 int
  31 FUNCTION(gsl_fft_real,transform) (BASE data[], const size_t stride, const size_t n,
  32                                   const TYPE(gsl_fft_real_wavetable) * wavetable,
  33                                   TYPE(gsl_fft_real_workspace) * work)
  34 {
  35   const size_t nf = wavetable->nf;
  36
  37   size_t i;
  38
  39   size_t q, product = 1;
  40   size_t tskip;
  41   size_t product_1;
  42
  43   BASE *const scratch = work->scratch;
  44   TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4;
  45
  46   size_t state = 0;
  47   BASE *in = data;
  48   size_t istride = stride ;
  49   BASE *out = scratch;
  50   size_t ostride = 1 ;
  51
  52   if (n == 0)
  53     {
  54       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  55     }
  56
  57   if (n == 1)
  58     {                           /* FFT of one data point is the identity */
  59       return 0;
  60     }
  61
  62   if (n != wavetable->n)
  63     {
  64       GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
  65     }
  66
  67   if (n != work->n)
  68     {
  69       GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
  70     }
  71
  72   for (i = 0; i < nf; i++)
  73     {
  74       const size_t factor = wavetable->factor[i];
  75       product_1 = product;
  76       product *= factor;
  77       q = n / product;
  78
  79       tskip = (product_1 + 1) / 2 - 1;
  80
  81       if (state == 0)
  82         {
  83           in = data;
  84           istride = stride;
  85           out = scratch;
  86           ostride = 1;
  87           state = 1;
  88         }
  89       else
  90         {
  91           in = scratch;
  92           istride = 1;
  93           out = data;
  94           ostride = stride;
  95           state = 0;
  96         }
  97
  98       if (factor == 2)
  99         {
 100           twiddle1 = wavetable->twiddle[i];
 101           FUNCTION(fft_real,pass_2) (in, istride, out, ostride, product, n, twiddle1);
 102         }
 103       else if (factor == 3)
 104         {
 105           twiddle1 = wavetable->twiddle[i];
 106           twiddle2 = twiddle1 + tskip;
 107           FUNCTION(fft_real,pass_3) (in, istride, out, ostride, product, n, twiddle1,
 108                                twiddle2);
 109         }
 110       else if (factor == 4)
 111         {
 112           twiddle1 = wavetable->twiddle[i];
 113           twiddle2 = twiddle1 + tskip;
 114           twiddle3 = twiddle2 + tskip;
 115           FUNCTION(fft_real,pass_4) (in, istride, out, ostride, product, n, twiddle1,
 116                                      twiddle2, twiddle3);
 117         }
 118       else if (factor == 5)
 119         {
 120           twiddle1 = wavetable->twiddle[i];
 121           twiddle2 = twiddle1 + tskip;
 122           twiddle3 = twiddle2 + tskip;
 123           twiddle4 = twiddle3 + tskip;
 124           FUNCTION(fft_real,pass_5) (in, istride, out, ostride, product, n, twiddle1,
 125                                      twiddle2, twiddle3, twiddle4);
 126         }
 127       else
 128         {
 129           twiddle1 = wavetable->twiddle[i];
 130           FUNCTION(fft_real,pass_n) (in, istride, out, ostride, factor, product, n,
 131                                      twiddle1);
 132         }
 133     }
 134
 135   if (state == 1)               /* copy results back from scratch to data */
 136     {
 137       for (i = 0; i < n; i++)
 138         {
 139           data[stride*i] = scratch[i] ;
 140         }
 141     }
 142
 143   return 0;
 144
 145 }