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

   1 /* fft/c_pass_2.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 static int
  21 FUNCTION(fft_complex,pass_2) (const BASE in[],
  22                               const size_t istride,
  23                               BASE out[],
  24                               const size_t ostride,
  25                               const gsl_fft_direction sign,
  26                               const size_t product,
  27                               const size_t n,
  28                               const TYPE(gsl_complex) twiddle[])
  29 {
  30   size_t i = 0, j = 0;
  31   size_t k, k1;
  32
  33   const size_t factor = 2;
  34   const size_t m = n / factor;
  35   const size_t q = n / product;
  36   const size_t product_1 = product / factor;
  37   const size_t jump = (factor - 1) * product_1;
  38
  39   for (k = 0; k < q; k++)
  40     {
  41       ATOMIC w_real, w_imag;
  42
  43       if (k == 0)
  44         {
  45           w_real = 1.0;
  46           w_imag = 0.0;
  47         }
  48       else
  49         {
  50           if (sign == gsl_fft_forward)
  51             {
  52               /* forward tranform */
  53               w_real = GSL_REAL(twiddle[k - 1]);
  54               w_imag = GSL_IMAG(twiddle[k - 1]);
  55             }
  56           else
  57             {
  58               /* backward tranform: w -> conjugate(w) */
  59               w_real = GSL_REAL(twiddle[k - 1]);
  60               w_imag = -GSL_IMAG(twiddle[k - 1]);
  61             }
  62         }
  63
  64       for (k1 = 0; k1 < product_1; k1++)
  65         {
  66           const ATOMIC z0_real = REAL(in,istride,i);
  67           const ATOMIC z0_imag = IMAG(in,istride,i);
  68
  69           const ATOMIC z1_real = REAL(in,istride,i+m);
  70           const ATOMIC z1_imag = IMAG(in,istride,i+m);
  71
  72           /* compute x = W(2) z */
  73
  74           /* x0 = z0 + z1 */
  75           const ATOMIC x0_real = z0_real + z1_real;
  76           const ATOMIC x0_imag = z0_imag + z1_imag;
  77
  78           /* x1 = z0 - z1 */
  79           const ATOMIC x1_real = z0_real - z1_real;
  80           const ATOMIC x1_imag = z0_imag - z1_imag;
  81
  82           /* apply twiddle factors */
  83
  84           /* out0 = 1 * x0 */
  85           REAL(out,ostride,j) = x0_real;
  86           IMAG(out,ostride,j) = x0_imag;
  87
  88           /* out1 = w * x1 */
  89           REAL(out,ostride,j+product_1) = w_real * x1_real - w_imag * x1_imag;
  90           IMAG(out,ostride,j+product_1) = w_real * x1_imag + w_imag * x1_real;
  91
  92           i++;
  93           j++;
  94         }
  95       j += jump;
  96     }
  97   return 0;
  98 }