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

   1 /* fft/c_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 "c_pass.h"
  21
  22 int
  23 FUNCTION(gsl_fft_complex,forward) (TYPE(gsl_complex_packed_array) data,
  24                                    const size_t stride,
  25                                    const size_t n,
  26                                    const TYPE(gsl_fft_complex_wavetable) * wavetable,
  27                                    TYPE(gsl_fft_complex_workspace) * work)
  28 {
  29   gsl_fft_direction sign = gsl_fft_forward;
  30   int status = FUNCTION(gsl_fft_complex,transform) (data, stride, n,
  31                                                     wavetable, work, sign);
  32   return status;
  33 }
  34
  35 int
  36 FUNCTION(gsl_fft_complex,backward) (TYPE(gsl_complex_packed_array) data,
  37                                     const size_t stride,
  38                                     const size_t n,
  39                                     const TYPE(gsl_fft_complex_wavetable) * wavetable,
  40                                     TYPE(gsl_fft_complex_workspace) * work)
  41 {
  42   gsl_fft_direction sign = gsl_fft_backward;
  43   int status = FUNCTION(gsl_fft_complex,transform) (data, stride, n,
  44                                                     wavetable, work, sign);
  45   return status;
  46 }
  47
  48 int
  49 FUNCTION(gsl_fft_complex,inverse) (TYPE(gsl_complex_packed_array) data,
  50                                    const size_t stride,
  51                                    const size_t n,
  52                                    const TYPE(gsl_fft_complex_wavetable) * wavetable,
  53                                    TYPE(gsl_fft_complex_workspace) * work)
  54 {
  55   gsl_fft_direction sign = gsl_fft_backward;
  56   int status = FUNCTION(gsl_fft_complex,transform) (data, stride, n,
  57                                                     wavetable, work, sign);
  58
  59   if (status)
  60     {
  61       return status;
  62     }
  63
  64   /* normalize inverse fft with 1/n */
  65
  66   {
  67     const ATOMIC norm = ONE / (ATOMIC)n;
  68     size_t i;
  69     for (i = 0; i < n; i++)
  70       {
  71         REAL(data,stride,i) *= norm;
  72         IMAG(data,stride,i) *= norm;
  73       }
  74   }
  75   return status;
  76 }
  77
  78 int
  79 FUNCTION(gsl_fft_complex,transform) (TYPE(gsl_complex_packed_array) data,
  80                                      const size_t stride,
  81                                      const size_t n,
  82                                      const TYPE(gsl_fft_complex_wavetable) * wavetable,
  83                                      TYPE(gsl_fft_complex_workspace) * work,
  84                                      const gsl_fft_direction sign)
  85 {
  86   const size_t nf = wavetable->nf;
  87
  88   size_t i;
  89
  90   size_t q, product = 1;
  91
  92   TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4,
  93     *twiddle5, *twiddle6;
  94
  95   size_t state = 0;
  96
  97   BASE * const scratch = work->scratch;
  98
  99   BASE * in = data;
 100   size_t istride = stride;
 101
 102   BASE * out = scratch;
 103   size_t ostride = 1;
 104
 105   if (n == 0)
 106     {
 107       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
 108     }
 109
 110   if (n == 1)
 111     {                           /* FFT of 1 data point is the identity */
 112       return 0;
 113     }
 114
 115   if (n != wavetable->n)
 116     {
 117       GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
 118     }
 119
 120   if (n != work->n)
 121     {
 122       GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
 123     }
 124
 125   for (i = 0; i < nf; i++)
 126     {
 127       const size_t factor = wavetable->factor[i];
 128       product *= factor;
 129       q = n / product;
 130
 131       if (state == 0)
 132         {
 133           in = data;
 134           istride = stride;
 135           out = scratch;
 136           ostride = 1;
 137           state = 1;
 138         }
 139       else
 140         {
 141           in = scratch;
 142           istride = 1;
 143           out = data;
 144           ostride = stride;
 145           state = 0;
 146         }
 147
 148       if (factor == 2)
 149         {
 150           twiddle1 = wavetable->twiddle[i];
 151           FUNCTION(fft_complex,pass_2) (in, istride, out, ostride, sign,
 152                                         product, n, twiddle1);
 153         }
 154       else if (factor == 3)
 155         {
 156           twiddle1 = wavetable->twiddle[i];
 157           twiddle2 = twiddle1 + q;
 158           FUNCTION(fft_complex,pass_3) (in, istride, out, ostride, sign,
 159                                         product, n, twiddle1, twiddle2);
 160         }
 161       else if (factor == 4)
 162         {
 163           twiddle1 = wavetable->twiddle[i];
 164           twiddle2 = twiddle1 + q;
 165           twiddle3 = twiddle2 + q;
 166           FUNCTION(fft_complex,pass_4) (in, istride, out, ostride, sign,
 167                                         product, n, twiddle1, twiddle2,
 168                                         twiddle3);
 169         }
 170       else if (factor == 5)
 171         {
 172           twiddle1 = wavetable->twiddle[i];
 173           twiddle2 = twiddle1 + q;
 174           twiddle3 = twiddle2 + q;
 175           twiddle4 = twiddle3 + q;
 176           FUNCTION(fft_complex,pass_5) (in, istride, out, ostride, sign,
 177                                         product, n, twiddle1, twiddle2,
 178                                         twiddle3, twiddle4);
 179         }
 180       else if (factor == 6)
 181         {
 182           twiddle1 = wavetable->twiddle[i];
 183           twiddle2 = twiddle1 + q;
 184           twiddle3 = twiddle2 + q;
 185           twiddle4 = twiddle3 + q;
 186           twiddle5 = twiddle4 + q;
 187           FUNCTION(fft_complex,pass_6) (in, istride, out, ostride, sign,
 188                                         product, n, twiddle1, twiddle2,
 189                                         twiddle3, twiddle4, twiddle5);
 190         }
 191       else if (factor == 7)
 192         {
 193           twiddle1 = wavetable->twiddle[i];
 194           twiddle2 = twiddle1 + q;
 195           twiddle3 = twiddle2 + q;
 196           twiddle4 = twiddle3 + q;
 197           twiddle5 = twiddle4 + q;
 198           twiddle6 = twiddle5 + q;
 199           FUNCTION(fft_complex,pass_7) (in, istride, out, ostride, sign,
 200                                         product, n, twiddle1, twiddle2,
 201                                         twiddle3, twiddle4, twiddle5,
 202                                         twiddle6);
 203         }
 204       else
 205         {
 206           twiddle1 = wavetable->twiddle[i];
 207           FUNCTION(fft_complex,pass_n) (in, istride, out, ostride, sign,
 208                                         factor, product, n, twiddle1);
 209         }
 210     }
 211
 212   if (state == 1)               /* copy results back from scratch to data */
 213     {
 214       for (i = 0; i < n; i++)
 215         {
 216           REAL(data,stride,i) = REAL(scratch,1,i) ;
 217           IMAG(data,stride,i) = IMAG(scratch,1,i) ;
 218         }
 219     }
 220
 221   return 0;
 222
 223 }