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

   1 /* fft/hc_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 void
  21 FUNCTION(fft_halfcomplex,pass_2) (const BASE in[],
  22                                   const size_t istride,
  23                                   BASE out[],
  24                                   const size_t ostride,
  25                                   const size_t product,
  26                                   const size_t n,
  27                                   const TYPE(gsl_complex) twiddle[])
  28 {
  29   size_t i, j, k, k1, jump;
  30   size_t factor, q, m, product_1;
  31   i = 0;
  32   j = 0;
  33
  34   factor = 2;
  35   m = n / factor;
  36   q = n / product;
  37   product_1 = product / factor;
  38   jump = (factor - 1) * q;
  39
  40   for (k1 = 0; k1 < product_1; k1++)
  41     {
  42       const ATOMIC r0 = VECTOR(in,istride,2 * k1 * q);
  43       const ATOMIC r1 = VECTOR(in,istride,2 * k1 * q + 2 * q - 1);
  44
  45       const ATOMIC s0 = r0 + r1;
  46       const ATOMIC s1 = r0 - r1;
  47
  48       VECTOR(out,ostride,q * k1) = s0;
  49       VECTOR(out,ostride,q * k1 + m) = s1;
  50     }
  51
  52   if (q == 1)
  53     return;
  54
  55   for (k = 1; k < (q + 1) / 2; k++)
  56     {
  57       const ATOMIC w_real = GSL_REAL(twiddle[k - 1]);
  58       const ATOMIC w_imag = GSL_IMAG(twiddle[k - 1]);
  59
  60       for (k1 = 0; k1 < product_1; k1++)
  61         {
  62           const size_t from0 = 2 * k1 * q + 2 * k - 1;
  63           const size_t from1 = 2 * k1 * q - 2 * k + 2 * q - 1;
  64
  65           const ATOMIC z0_real = VECTOR(in,istride,from0);
  66           const ATOMIC z0_imag = VECTOR(in,istride,from0 + 1);
  67
  68           const ATOMIC z1_real = VECTOR(in,istride,from1);
  69           const ATOMIC z1_imag = VECTOR(in,istride,from1 + 1);
  70
  71           /* compute x = W(2) z */
  72
  73           /* x0 = z0 + z1 */
  74           const ATOMIC x0_real = z0_real + z1_real;
  75           const ATOMIC x0_imag = z0_imag - z1_imag;
  76
  77           /* x1 = z0 - z1 */
  78           const ATOMIC x1_real = z0_real - z1_real;
  79           const ATOMIC x1_imag = z0_imag + z1_imag;
  80
  81           const size_t to0 = k1 * q + 2 * k - 1;
  82           const size_t to1 = to0 + m;
  83
  84           VECTOR(out,ostride,to0) = x0_real;
  85           VECTOR(out,ostride,to0 + 1) = x0_imag;
  86
  87           VECTOR(out,ostride,to1) = w_real * x1_real - w_imag * x1_imag;
  88           VECTOR(out,ostride,to1 + 1) = w_imag * x1_real + w_real * x1_imag;
  89
  90         }
  91     }
  92
  93   if (q % 2 == 1)
  94     return;
  95
  96   for (k1 = 0; k1 < product_1; k1++)
  97     {
  98       const size_t from0 = 2 * k1 * q + q - 1;
  99       const size_t to0 = k1 * q + q - 1;
 100       const size_t to1 = to0 + m;
 101
 102       VECTOR(out,ostride,to0) = 2 * VECTOR(in,istride,from0);
 103       VECTOR(out,ostride,to1) = -2 * VECTOR(in,istride,from0 + 1);
 104     }
 105   return;
 106 }