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

   1 /* fft/c_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_complex,radix2_forward) (TYPE(gsl_complex_packed_array) data,
  22                                           const size_t stride, const size_t n)
  23 {
  24   gsl_fft_direction sign = gsl_fft_forward;
  25   int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
  26   return status;
  27 }
  28
  29 int
  30 FUNCTION(gsl_fft_complex,radix2_backward) (TYPE(gsl_complex_packed_array) data,
  31                                            const size_t stride, const size_t n)
  32 {
  33   gsl_fft_direction sign = gsl_fft_backward;
  34   int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
  35   return status;
  36 }
  37
  38 int
  39 FUNCTION(gsl_fft_complex,radix2_inverse) (TYPE(gsl_complex_packed_array) data,
  40                                           const size_t stride, const size_t n)
  41 {
  42   gsl_fft_direction sign = gsl_fft_backward;
  43   int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
  44
  45   if (status)
  46     {
  47       return status;
  48     }
  49
  50   /* normalize inverse fft with 1/n */
  51
  52   {
  53     const ATOMIC norm = 1.0 / n;
  54     size_t i;
  55     for (i = 0; i < n; i++)
  56       {
  57         REAL(data,stride,i) *= norm;
  58         IMAG(data,stride,i) *= norm;
  59       }
  60   }
  61
  62   return status;
  63 }
  64
  65
  66
  67 int
  68 FUNCTION(gsl_fft_complex,radix2_transform) (TYPE(gsl_complex_packed_array) data,
  69                                             const size_t stride,
  70                                             const size_t n,
  71                                             const gsl_fft_direction sign)
  72 {
  73   int result ;
  74   size_t dual;
  75   size_t bit;
  76   size_t logn = 0;
  77   int status;
  78
  79   if (n == 1) /* identity operation */
  80     {
  81       return 0 ;
  82     }
  83
  84   /* make sure that n is a power of 2 */
  85
  86   result = fft_binary_logn(n) ;
  87
  88   if (result == -1)
  89     {
  90       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
  91     }
  92   else
  93     {
  94       logn = result ;
  95     }
  96
  97   /* bit reverse the ordering of input data for decimation in time algorithm */
  98
  99   status = FUNCTION(fft_complex,bitreverse_order) (data, stride, n, logn) ;
 100
 101   /* apply fft recursion */
 102
 103   dual = 1;
 104
 105   for (bit = 0; bit < logn; bit++)
 106     {
 107       ATOMIC w_real = 1.0;
 108       ATOMIC w_imag = 0.0;
 109
 110       const double theta = 2.0 * ((int) sign) * M_PI / (2.0 * (double) dual);
 111
 112       const ATOMIC s = sin (theta);
 113       const ATOMIC t = sin (theta / 2.0);
 114       const ATOMIC s2 = 2.0 * t * t;
 115
 116       size_t a, b;
 117
 118       /* a = 0 */
 119
 120       for (b = 0; b < n; b += 2 * dual)
 121         {
 122           const size_t i = b ;
 123           const size_t j = b + dual;
 124
 125           const ATOMIC z1_real = REAL(data,stride,j) ;
 126           const ATOMIC z1_imag = IMAG(data,stride,j) ;
 127
 128           const ATOMIC wd_real = z1_real ;
 129           const ATOMIC wd_imag = z1_imag ;
 130
 131           REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
 132           IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
 133           REAL(data,stride,i) += wd_real;
 134           IMAG(data,stride,i) += wd_imag;
 135         }
 136
 137       /* a = 1 .. (dual-1) */
 138
 139       for (a = 1; a < dual; a++)
 140         {
 141
 142           /* trignometric recurrence for w-> exp(i theta) w */
 143
 144           {
 145             const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
 146             const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
 147             w_real = tmp_real;
 148             w_imag = tmp_imag;
 149           }
 150
 151           for (b = 0; b < n; b += 2 * dual)
 152             {
 153               const size_t i = b + a;
 154               const size_t j = b + a + dual;
 155
 156               const ATOMIC z1_real = REAL(data,stride,j) ;
 157               const ATOMIC z1_imag = IMAG(data,stride,j) ;
 158
 159               const ATOMIC wd_real = w_real * z1_real - w_imag * z1_imag;
 160               const ATOMIC wd_imag = w_real * z1_imag + w_imag * z1_real;
 161
 162               REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
 163               IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
 164               REAL(data,stride,i) += wd_real;
 165               IMAG(data,stride,i) += wd_imag;
 166             }
 167         }
 168       dual *= 2;
 169     }
 170
 171   return 0;
 172
 173 }
 174
 175
 176 int
 177 FUNCTION(gsl_fft_complex,radix2_dif_forward) (TYPE(gsl_complex_packed_array) data,
 178                                               const size_t stride,
 179                                               const size_t n)
 180 {
 181   gsl_fft_direction sign = gsl_fft_forward;
 182   int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
 183   return status;
 184 }
 185
 186 int
 187 FUNCTION(gsl_fft_complex,radix2_dif_backward) (TYPE(gsl_complex_packed_array) data,
 188                                                const size_t stride,
 189                                                const size_t n)
 190 {
 191   gsl_fft_direction sign = gsl_fft_backward;
 192   int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
 193   return status;
 194 }
 195
 196 int
 197 FUNCTION(gsl_fft_complex,radix2_dif_inverse) (TYPE(gsl_complex_packed_array) data,
 198                                               const size_t stride,
 199                                               const size_t n)
 200 {
 201   gsl_fft_direction sign = gsl_fft_backward;
 202   int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
 203
 204   if (status)
 205     {
 206       return status;
 207     }
 208
 209   /* normalize inverse fft with 1/n */
 210
 211   {
 212     const ATOMIC norm = 1.0 / n;
 213     size_t i;
 214     for (i = 0; i < n; i++)
 215       {
 216         REAL(data,stride,i) *= norm;
 217         IMAG(data,stride,i) *= norm;
 218       }
 219   }
 220
 221   return status;
 222 }
 223
 224 int
 225 FUNCTION(gsl_fft_complex,radix2_dif_transform) (TYPE(gsl_complex_packed_array) data,
 226                                       const size_t stride,
 227                                       const size_t n,
 228                                       const gsl_fft_direction sign)
 229 {
 230   int result ;
 231   size_t dual;
 232   size_t bit;
 233   size_t logn = 0;
 234   int status;
 235
 236   if (n == 1) /* identity operation */
 237     {
 238       return 0 ;
 239     }
 240
 241   /* make sure that n is a power of 2 */
 242
 243   result = fft_binary_logn(n) ;
 244
 245   if (result == -1)
 246     {
 247       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
 248     }
 249   else
 250     {
 251       logn = result ;
 252     }
 253
 254   /* apply fft recursion */
 255
 256   dual = n / 2;
 257
 258   for (bit = 0; bit < logn; bit++)
 259     {
 260       ATOMIC w_real = 1.0;
 261       ATOMIC w_imag = 0.0;
 262
 263       const double theta = 2.0 * ((int) sign) * M_PI / ((double) (2 * dual));
 264
 265       const ATOMIC s = sin (theta);
 266       const ATOMIC t = sin (theta / 2.0);
 267       const ATOMIC s2 = 2.0 * t * t;
 268
 269       size_t a, b;
 270
 271       for (b = 0; b < dual; b++)
 272         {
 273           for (a = 0; a < n; a+= 2 * dual)
 274             {
 275               const size_t i = b + a;
 276               const size_t j = b + a + dual;
 277
 278               const ATOMIC t1_real = REAL(data,stride,i) + REAL(data,stride,j);
 279               const ATOMIC t1_imag = IMAG(data,stride,i) + IMAG(data,stride,j);
 280               const ATOMIC t2_real = REAL(data,stride,i) - REAL(data,stride,j);
 281               const ATOMIC t2_imag = IMAG(data,stride,i) - IMAG(data,stride,j);
 282
 283               REAL(data,stride,i) = t1_real;
 284               IMAG(data,stride,i) = t1_imag;
 285               REAL(data,stride,j) = w_real*t2_real - w_imag * t2_imag;
 286               IMAG(data,stride,j) = w_real*t2_imag + w_imag * t2_real;
 287             }
 288
 289           /* trignometric recurrence for w-> exp(i theta) w */
 290
 291           {
 292             const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
 293             const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
 294             w_real = tmp_real;
 295             w_imag = tmp_imag;
 296           }
 297         }
 298       dual /= 2;
 299     }
 300
 301   /* bit reverse the ordering of output data for decimation in
 302      frequency algorithm */
 303
 304   status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
 305
 306   return 0;
 307
 308 }
 309
 310
 311
 312
 313
 314
 315
 316