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

   1 /* fft/hc_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_halfcomplex.h>
  27
  28 #include "hc_pass.h"
  29
  30 int
  31 FUNCTION(gsl_fft_halfcomplex,backward) (BASE data[], const size_t stride,
  32                                         const size_t n,
  33                                         const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
  34                                         TYPE(gsl_fft_real_workspace) * work)
  35 {
  36   int status = FUNCTION(gsl_fft_halfcomplex,transform) (data, stride, n, wavetable, work) ;
  37   return status ;
  38 }
  39
  40 int
  41 FUNCTION(gsl_fft_halfcomplex,inverse) (BASE data[], const size_t stride,
  42                                        const size_t n,
  43                                        const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
  44                                        TYPE(gsl_fft_real_workspace) * work)
  45 {
  46   int status = FUNCTION(gsl_fft_halfcomplex,transform) (data, stride, n, wavetable, work);
  47
  48   if (status)
  49     {
  50       return status;
  51     }
  52
  53   /* normalize inverse fft with 1/n */
  54
  55   {
  56     const double norm = 1.0 / n;
  57     size_t i;
  58     for (i = 0; i < n; i++)
  59       {
  60         data[stride*i] *= norm;
  61       }
  62   }
  63   return status;
  64 }
  65
  66 int
  67 FUNCTION(gsl_fft_halfcomplex,transform) (BASE data[], const size_t stride, const size_t n,
  68                                          const TYPE(gsl_fft_halfcomplex_wavetable) * wavetable,
  69                                          TYPE(gsl_fft_real_workspace) * work)
  70 {
  71   BASE * const scratch = work->scratch;
  72
  73   BASE * in;
  74   BASE * out;
  75   size_t istride, ostride ;
  76
  77
  78   size_t factor, product, q;
  79   size_t i;
  80   size_t nf;
  81   int state;
  82   int product_1;
  83   int tskip;
  84   TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4;
  85
  86   if (n == 0)
  87     {
  88       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  89     }
  90
  91   if (n == 1)
  92     {                           /* FFT of one data point is the identity */
  93       return 0;
  94     }
  95
  96   if (n != wavetable->n)
  97     {
  98       GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
  99     }
 100
 101   if (n != work->n)
 102     {
 103       GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
 104     }
 105
 106   nf = wavetable->nf;
 107   product = 1;
 108   state = 0;
 109
 110   for (i = 0; i < nf; i++)
 111     {
 112       factor = wavetable->factor[i];
 113       product_1 = product;
 114       product *= factor;
 115       q = n / product;
 116
 117       tskip = (q + 1) / 2 - 1;
 118
 119       if (state == 0)
 120         {
 121           in = data;
 122           istride = stride;
 123           out = scratch;
 124           ostride = 1;
 125           state = 1;
 126         }
 127       else
 128         {
 129           in = scratch;
 130           istride = 1;
 131           out = data;
 132           ostride = stride;
 133           state = 0;
 134         }
 135
 136       if (factor == 2)
 137         {
 138           twiddle1 = wavetable->twiddle[i];
 139           FUNCTION(fft_halfcomplex,pass_2) (in, istride, out, ostride,
 140                                             product, n, twiddle1);
 141         }
 142       else if (factor == 3)
 143         {
 144           twiddle1 = wavetable->twiddle[i];
 145           twiddle2 = twiddle1 + tskip;
 146           FUNCTION(fft_halfcomplex,pass_3) (in, istride, out, ostride,
 147                                             product, n, twiddle1, twiddle2);
 148         }
 149       else if (factor == 4)
 150         {
 151           twiddle1 = wavetable->twiddle[i];
 152           twiddle2 = twiddle1 + tskip;
 153           twiddle3 = twiddle2 + tskip;
 154           FUNCTION(fft_halfcomplex,pass_4) (in, istride, out, ostride,
 155                                             product, n, twiddle1, twiddle2,
 156                                             twiddle3);
 157         }
 158       else if (factor == 5)
 159         {
 160           twiddle1 = wavetable->twiddle[i];
 161           twiddle2 = twiddle1 + tskip;
 162           twiddle3 = twiddle2 + tskip;
 163           twiddle4 = twiddle3 + tskip;
 164           FUNCTION(fft_halfcomplex,pass_5) (in, istride, out, ostride,
 165                                             product, n, twiddle1, twiddle2,
 166                                             twiddle3, twiddle4);
 167         }
 168       else
 169         {
 170           twiddle1 = wavetable->twiddle[i];
 171           FUNCTION(fft_halfcomplex,pass_n) (in, istride, out, ostride,
 172                                             factor, product, n, twiddle1);
 173         }
 174     }
 175
 176   if (state == 1)               /* copy results back from scratch to data */
 177     {
 178       for (i = 0; i < n; i++)
 179         {
 180           data[stride*i] = scratch[i] ;
 181         }
 182     }
 183
 184   return 0;
 185
 186 }
 187
 188