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

   1 /* fft/hc_radix2.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 int
  21 FUNCTION(gsl_fft_halfcomplex,radix2_backward) (BASE data[],
  22                                                const size_t stride,
  23                                                const size_t n)
  24 {
  25   int status = FUNCTION(gsl_fft_halfcomplex,radix2_transform) (data, stride, n) ;
  26   return status ;
  27 }
  28
  29 int
  30 FUNCTION(gsl_fft_halfcomplex,radix2_inverse) (BASE data[],
  31                                               const size_t stride,
  32                                               const size_t n)
  33 {
  34   int status = FUNCTION(gsl_fft_halfcomplex,radix2_transform) (data, stride, n);
  35
  36   if (status)
  37     {
  38       return status;
  39     }
  40
  41   /* normalize inverse fft with 1/n */
  42
  43   {
  44     const ATOMIC norm = 1.0 / n;
  45     size_t i;
  46     for (i = 0; i < n; i++)
  47       {
  48         data[stride*i] *= norm;
  49       }
  50   }
  51   return status;
  52 }
  53
  54 int
  55 FUNCTION(gsl_fft_halfcomplex,radix2_transform) (BASE data[],
  56                                                 const size_t stride,
  57                                                 const size_t n)
  58 {
  59   int result ;
  60   size_t p, p_1, q;
  61   size_t i;
  62   size_t logn = 0;
  63   int status;
  64
  65   if (n == 1) /* identity operation */
  66     {
  67       return 0 ;
  68     }
  69
  70   /* make sure that n is a power of 2 */
  71
  72   result = fft_binary_logn(n) ;
  73
  74   if (result == -1)
  75     {
  76       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
  77     }
  78   else
  79     {
  80       logn = result ;
  81     }
  82
  83   /* apply fft recursion */
  84
  85   p = n; q = 1 ; p_1 = n/2 ;
  86
  87   for (i = 1; i <= logn; i++)
  88     {
  89       size_t a, b;
  90
  91       /* a = 0 */
  92
  93       for (b = 0; b < q; b++)
  94         {
  95           const ATOMIC z0 = VECTOR(data,stride,b*p);
  96           const ATOMIC z1 = VECTOR(data,stride,b*p + p_1);
  97
  98           const ATOMIC t0_real = z0 + z1 ;
  99           const ATOMIC t1_real = z0 - z1 ;
 100
 101           VECTOR(data,stride,b*p) = t0_real;
 102           VECTOR(data,stride,b*p + p_1) = t1_real ;
 103         }
 104
 105       /* a = 1 ... p_{i-1}/2 - 1 */
 106
 107       {
 108         ATOMIC w_real = 1.0;
 109         ATOMIC w_imag = 0.0;
 110
 111         const ATOMIC theta = 2.0 * M_PI / p;
 112
 113         const ATOMIC s = sin (theta);
 114         const ATOMIC t = sin (theta / 2.0);
 115         const ATOMIC s2 = 2.0 * t * t;
 116
 117         for (a = 1; a < (p_1)/2; a++)
 118           {
 119             /* trignometric recurrence for w-> exp(i theta) w */
 120
 121             {
 122               const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
 123               const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
 124               w_real = tmp_real;
 125               w_imag = tmp_imag;
 126             }
 127
 128             for (b = 0; b < q; b++)
 129               {
 130                 ATOMIC z0_real = VECTOR(data,stride,b*p + a) ;
 131                 ATOMIC z0_imag = VECTOR(data,stride,b*p + p - a) ;
 132                 ATOMIC z1_real = VECTOR(data,stride,b*p + p_1 - a) ;
 133                 ATOMIC z1_imag = -VECTOR(data,stride,b*p + p_1 + a) ;
 134
 135                 /* t0 = z0 + z1 */
 136
 137                 ATOMIC t0_real = z0_real + z1_real;
 138                 ATOMIC t0_imag = z0_imag + z1_imag;
 139
 140                 /* t1 = (z0 - z1) */
 141
 142                 ATOMIC t1_real = z0_real -  z1_real;
 143                 ATOMIC t1_imag = z0_imag -  z1_imag;
 144
 145                 VECTOR(data,stride,b*p + a) = t0_real ;
 146                 VECTOR(data,stride,b*p + p_1 - a) = t0_imag ;
 147
 148                 VECTOR(data,stride,b*p + p_1 + a) = (w_real * t1_real - w_imag * t1_imag) ;
 149                 VECTOR(data,stride,b*p + p - a) = (w_real * t1_imag + w_imag * t1_real) ;
 150               }
 151           }
 152       }
 153
 154       if (p_1 >  1) {
 155         for (b = 0; b < q; b++) {
 156           VECTOR(data,stride,b*p + p_1/2) *= 2 ;
 157           VECTOR(data,stride,b*p + p_1 + p_1/2) *= -2 ;
 158         }
 159       }
 160
 161       p_1 = p_1 / 2 ;
 162       p = p / 2 ;
 163       q = q * 2 ;
 164     }
 165
 166   /* bit reverse the ordering of output data for decimation in
 167      frequency algorithm */
 168
 169   status = FUNCTION(fft_real,bitreverse_order)(data, stride, n, logn) ;
 170
 171   return 0;
 172
 173 }