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

   1 /* fft/real_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_real,radix2_transform) (BASE data[], const size_t stride,  const size_t n)
  22 {
  23   int result ;
  24   size_t p, p_1, q;
  25   size_t i;
  26   size_t logn = 0;
  27   int status;
  28
  29   if (n == 1) /* identity operation */
  30     {
  31       return 0 ;
  32     }
  33
  34   /* make sure that n is a power of 2 */
  35
  36   result = fft_binary_logn(n) ;
  37
  38   if (result == -1)
  39     {
  40       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
  41     }
  42   else
  43     {
  44       logn = result ;
  45     }
  46
  47   /* bit reverse the ordering of input data for decimation in time algorithm */
  48
  49   status = FUNCTION(fft_real,bitreverse_order)(data, stride, n, logn) ;
  50
  51   /* apply fft recursion */
  52
  53   p = 1; q = n ;
  54
  55   for (i = 1; i <= logn; i++)
  56     {
  57       size_t a, b;
  58
  59       p_1 = p ;
  60       p = 2 * p ;
  61       q = q / 2 ;
  62
  63       /* a = 0 */
  64
  65       for (b = 0; b < q; b++)
  66         {
  67           ATOMIC t0_real = VECTOR(data,stride,b*p) + VECTOR(data,stride,b*p + p_1) ;
  68           ATOMIC t1_real = VECTOR(data,stride,b*p) - VECTOR(data,stride,b*p + p_1) ;
  69
  70           VECTOR(data,stride,b*p) = t0_real ;
  71           VECTOR(data,stride,b*p + p_1) = t1_real ;
  72         }
  73
  74       /* a = 1 ... p_{i-1}/2 - 1 */
  75
  76       {
  77         ATOMIC w_real = 1.0;
  78         ATOMIC w_imag = 0.0;
  79
  80         const double theta = - 2.0 * M_PI / p;
  81
  82         const ATOMIC s = sin (theta);
  83         const ATOMIC t = sin (theta / 2.0);
  84         const ATOMIC s2 = 2.0 * t * t;
  85
  86         for (a = 1; a < (p_1)/2; a++)
  87           {
  88             /* trignometric recurrence for w-> exp(i theta) w */
  89
  90             {
  91               const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
  92               const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
  93               w_real = tmp_real;
  94               w_imag = tmp_imag;
  95             }
  96
  97             for (b = 0; b < q; b++)
  98               {
  99                 ATOMIC z0_real = VECTOR(data,stride,b*p + a) ;
 100                 ATOMIC z0_imag = VECTOR(data,stride,b*p + p_1 - a) ;
 101                 ATOMIC z1_real = VECTOR(data,stride,b*p + p_1 + a) ;
 102                 ATOMIC z1_imag = VECTOR(data,stride,b*p + p - a) ;
 103
 104                 /* t0 = z0 + w * z1 */
 105
 106                 ATOMIC t0_real = z0_real + w_real * z1_real - w_imag * z1_imag;
 107                 ATOMIC t0_imag = z0_imag + w_real * z1_imag + w_imag * z1_real;
 108
 109                 /* t1 = z0 - w * z1 */
 110
 111                 ATOMIC t1_real = z0_real - w_real * z1_real + w_imag * z1_imag;
 112                 ATOMIC t1_imag = z0_imag - w_real * z1_imag - w_imag * z1_real;
 113
 114                 VECTOR(data,stride,b*p + a) = t0_real ;
 115                 VECTOR(data,stride,b*p + p - a) = t0_imag ;
 116
 117                 VECTOR(data,stride,b*p + p_1 - a) = t1_real ;
 118                 VECTOR(data,stride,b*p + p_1 + a) = -t1_imag ;
 119               }
 120           }
 121       }
 122
 123       if (p_1 >  1)
 124         {
 125           for (b = 0; b < q; b++)
 126             {
 127               /* a = p_{i-1}/2 */
 128
 129               VECTOR(data,stride,b*p + p - p_1/2) *= -1 ;
 130             }
 131         }
 132     }
 133   return 0;
 134 }