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

   1 /* fft/c_pass_3.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_3) (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) * twiddle1,
  29                               const TYPE(gsl_complex) * twiddle2)
  30 {
  31   size_t i = 0, j = 0;
  32   size_t k, k1;
  33
  34   const size_t factor = 3;
  35   const size_t m = n / factor;
  36   const size_t q = n / product;
  37   const size_t product_1 = product / factor;
  38   const size_t jump = (factor - 1) * product_1;
  39
  40   const ATOMIC tau = sqrt (3.0) / 2.0;
  41
  42   for (k = 0; k < q; k++)
  43     {
  44       ATOMIC w1_real, w1_imag, w2_real, w2_imag;
  45
  46       if (k == 0)
  47         {
  48           w1_real = 1.0;
  49           w1_imag = 0.0;
  50           w2_real = 1.0;
  51           w2_imag = 0.0;
  52         }
  53       else
  54         {
  55           if (sign == gsl_fft_forward)
  56             {
  57               /* forward tranform */
  58               w1_real = GSL_REAL(twiddle1[k - 1]);
  59               w1_imag = GSL_IMAG(twiddle1[k - 1]);
  60               w2_real = GSL_REAL(twiddle2[k - 1]);
  61               w2_imag = GSL_IMAG(twiddle2[k - 1]);
  62             }
  63           else
  64             {
  65               /* backward tranform: w -> conjugate(w) */
  66               w1_real = GSL_REAL(twiddle1[k - 1]);
  67               w1_imag = -GSL_IMAG(twiddle1[k - 1]);
  68               w2_real = GSL_REAL(twiddle2[k - 1]);
  69               w2_imag = -GSL_IMAG(twiddle2[k - 1]);
  70             }
  71         }
  72
  73       for (k1 = 0; k1 < product_1; k1++)
  74         {
  75           const ATOMIC z0_real = REAL(in,istride,i);
  76           const ATOMIC z0_imag = IMAG(in,istride,i);
  77           const ATOMIC z1_real = REAL(in,istride,i+m);
  78           const ATOMIC z1_imag = IMAG(in,istride,i+m);
  79           const ATOMIC z2_real = REAL(in,istride,i+2*m);
  80           const ATOMIC z2_imag = IMAG(in,istride,i+2*m);
  81
  82           /* compute x = W(3) z */
  83
  84           /* t1 = z1 + z2 */
  85           const ATOMIC t1_real = z1_real + z2_real;
  86           const ATOMIC t1_imag = z1_imag + z2_imag;
  87
  88           /* t2 = z0 - t1/2 */
  89           const ATOMIC t2_real = z0_real - t1_real / 2.0;
  90           const ATOMIC t2_imag = z0_imag - t1_imag / 2.0;
  91
  92           /* t3 = (+/-) sin(pi/3)*(z1 - z2) */
  93           const ATOMIC t3_real = ((int) sign) * tau * (z1_real - z2_real);
  94           const ATOMIC t3_imag = ((int) sign) * tau * (z1_imag - z2_imag);
  95
  96           /* x0 = z0 + t1 */
  97           const ATOMIC x0_real = z0_real + t1_real;
  98           const ATOMIC x0_imag = z0_imag + t1_imag;
  99
 100           /* x1 = t2 + i t3 */
 101           const ATOMIC x1_real = t2_real - t3_imag;
 102           const ATOMIC x1_imag = t2_imag + t3_real;
 103
 104           /* x2 = t2 - i t3 */
 105           const ATOMIC x2_real = t2_real + t3_imag;
 106           const ATOMIC x2_imag = t2_imag - t3_real;
 107
 108           /* apply twiddle factors */
 109
 110           /* to0 = 1 * x0 */
 111           REAL(out,ostride,j) = x0_real;
 112           IMAG(out,ostride,j) = x0_imag;
 113
 114           /* to1 = w1 * x1 */
 115           REAL(out,ostride,j+product_1) = w1_real * x1_real - w1_imag * x1_imag;
 116           IMAG(out,ostride,j+product_1) = w1_real * x1_imag + w1_imag * x1_real;
 117
 118           /* to2 = w2 * x2 */
 119           REAL(out,ostride,j+2*product_1) = w2_real * x2_real - w2_imag * x2_imag;
 120           IMAG(out,ostride,j+2*product_1) = w2_real * x2_imag + w2_imag * x2_real;
 121
 122           i++; j++;
 123         }
 124       j += jump;
 125     }
 126   return 0;
 127 }