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

   1 /* fft/real_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_real,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 k, k1;
  30
  31   const size_t factor = 2;
  32   const size_t m = n / factor;
  33   const size_t q = n / product;
  34   const size_t product_1 = product / factor;
  35
  36   for (k1 = 0; k1 < q; k1++)
  37     {
  38       const size_t from0 = k1 * product_1;
  39       const size_t from1 = from0 + m;
  40
  41       const ATOMIC r0 = VECTOR(in,istride,from0);
  42       const ATOMIC r1 = VECTOR(in,istride,from1);
  43
  44       const ATOMIC s0 = r0 + r1;
  45       const ATOMIC s1 = r0 - r1;
  46
  47       const size_t to0 = product * k1;
  48       const size_t to1 = to0 + product - 1;
  49
  50       VECTOR(out,ostride,to0) = s0;
  51       VECTOR(out,ostride,to1) = s1;
  52     }
  53
  54   if (product_1 == 1)
  55     return;
  56
  57   for (k = 1; k < (product_1 + 1) / 2; k++)
  58     {
  59
  60       /* forward real transform: w -> conjugate(w) */
  61       const ATOMIC w_real = GSL_REAL(twiddle[k - 1]);
  62       const ATOMIC w_imag = -GSL_IMAG(twiddle[k - 1]);
  63
  64       for (k1 = 0; k1 < q; k1++)
  65         {
  66           const size_t from0 = k1 * product_1 + 2 * k - 1;
  67           const size_t from1 = from0 + m;
  68
  69           const ATOMIC f0_real = VECTOR(in,istride,from0);
  70           const ATOMIC f0_imag = VECTOR(in,istride,from0 + 1);
  71
  72           const ATOMIC f1_real = VECTOR(in,istride,from1);
  73           const ATOMIC f1_imag = VECTOR(in,istride,from1 + 1);
  74
  75           const ATOMIC z0_real = f0_real;
  76           const ATOMIC z0_imag = f0_imag;
  77
  78           const ATOMIC z1_real = w_real * f1_real - w_imag * f1_imag;
  79           const ATOMIC z1_imag = w_real * f1_imag + w_imag * f1_real;
  80
  81           /* compute x = W(2) z */
  82
  83           /* x0 = z0 + z1 */
  84           const ATOMIC x0_real = z0_real + z1_real;
  85           const ATOMIC x0_imag = z0_imag + z1_imag;
  86
  87           /* x1 = z0 - z1 */
  88           const ATOMIC x1_real = z0_real - z1_real;
  89           const ATOMIC x1_imag = z0_imag - z1_imag;
  90
  91           const size_t to0 = k1 * product + 2 * k - 1;
  92           const size_t to1 = k1 * product + product - 2 * k - 1;
  93
  94           VECTOR(out,ostride,to0) = x0_real;
  95           VECTOR(out,ostride,to0 + 1) = x0_imag;
  96
  97           /* stored in conjugate location */
  98           VECTOR(out,ostride,to1) = x1_real;
  99           VECTOR(out,ostride,to1 + 1) = -x1_imag;
 100         }
 101     }
 102
 103   if (product_1 % 2 == 1)
 104     return;
 105
 106   for (k1 = 0; k1 < q; k1++)
 107     {
 108       const size_t from0 = k1 * product_1 + product_1 - 1;
 109       const size_t from1 = from0 + m;
 110       const size_t to0 = k1 * product + product_1 - 1;
 111
 112       VECTOR(out,ostride,to0) = VECTOR(in,istride,from0);
 113       VECTOR(out,ostride,to0 + 1) = -VECTOR(in,istride,from1);
 114     }
 115   return;
 116 }