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

   1 /* fft/c_pass_6.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_6) (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                               const TYPE(gsl_complex) twiddle3[],
  31                               const TYPE(gsl_complex) twiddle4[],
  32                               const TYPE(gsl_complex) twiddle5[])
  33 {
  34
  35   size_t i = 0, j = 0;
  36   size_t k, k1;
  37
  38   const size_t factor = 6;
  39   const size_t m = n / factor;
  40   const size_t q = n / product;
  41   const size_t p_1 = product / factor;
  42   const size_t jump = (factor - 1) * p_1;
  43
  44   const ATOMIC tau = sqrt (3.0) / 2.0;
  45
  46   for (k = 0; k < q; k++)
  47     {
  48       ATOMIC w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag, w4_real,
  49         w4_imag, w5_real, w5_imag;
  50
  51       if (k == 0)
  52         {
  53           w1_real = 1.0;
  54           w1_imag = 0.0;
  55           w2_real = 1.0;
  56           w2_imag = 0.0;
  57           w3_real = 1.0;
  58           w3_imag = 0.0;
  59           w4_real = 1.0;
  60           w4_imag = 0.0;
  61           w5_real = 1.0;
  62           w5_imag = 0.0;
  63         }
  64       else
  65         {
  66           if (sign == gsl_fft_forward)
  67             {
  68               /* forward tranform */
  69               w1_real = GSL_REAL(twiddle1[k - 1]);
  70               w1_imag = GSL_IMAG(twiddle1[k - 1]);
  71               w2_real = GSL_REAL(twiddle2[k - 1]);
  72               w2_imag = GSL_IMAG(twiddle2[k - 1]);
  73               w3_real = GSL_REAL(twiddle3[k - 1]);
  74               w3_imag = GSL_IMAG(twiddle3[k - 1]);
  75               w4_real = GSL_REAL(twiddle4[k - 1]);
  76               w4_imag = GSL_IMAG(twiddle4[k - 1]);
  77               w5_real = GSL_REAL(twiddle5[k - 1]);
  78               w5_imag = GSL_IMAG(twiddle5[k - 1]);
  79             }
  80           else
  81             {
  82               /* backward tranform: w -> conjugate(w) */
  83               w1_real = GSL_REAL(twiddle1[k - 1]);
  84               w1_imag = -GSL_IMAG(twiddle1[k - 1]);
  85               w2_real = GSL_REAL(twiddle2[k - 1]);
  86               w2_imag = -GSL_IMAG(twiddle2[k - 1]);
  87               w3_real = GSL_REAL(twiddle3[k - 1]);
  88               w3_imag = -GSL_IMAG(twiddle3[k - 1]);
  89               w4_real = GSL_REAL(twiddle4[k - 1]);
  90               w4_imag = -GSL_IMAG(twiddle4[k - 1]);
  91               w5_real = GSL_REAL(twiddle5[k - 1]);
  92               w5_imag = -GSL_IMAG(twiddle5[k - 1]);
  93             }
  94         }
  95
  96       for (k1 = 0; k1 < p_1; k1++)
  97         {
  98           const ATOMIC z0_real = REAL(in,istride,i);
  99           const ATOMIC z0_imag = IMAG(in,istride,i);
 100           const ATOMIC z1_real = REAL(in,istride,i+m);
 101           const ATOMIC z1_imag = IMAG(in,istride,i+m);
 102           const ATOMIC z2_real = REAL(in,istride,i+2*m);
 103           const ATOMIC z2_imag = IMAG(in,istride,i+2*m);
 104           const ATOMIC z3_real = REAL(in,istride,i+3*m);
 105           const ATOMIC z3_imag = IMAG(in,istride,i+3*m);
 106           const ATOMIC z4_real = REAL(in,istride,i+4*m);
 107           const ATOMIC z4_imag = IMAG(in,istride,i+4*m);
 108           const ATOMIC z5_real = REAL(in,istride,i+5*m);
 109           const ATOMIC z5_imag = IMAG(in,istride,i+5*m);
 110
 111           /* compute x = W(6) z */
 112
 113           /* W(6) is a combination of sums and differences of W(3) acting
 114              on the even and odd elements of z */
 115
 116           /* ta1 = z2 + z4 */
 117           const ATOMIC ta1_real = z2_real + z4_real;
 118           const ATOMIC ta1_imag = z2_imag + z4_imag;
 119
 120           /* ta2 = z0 - ta1/2 */
 121           const ATOMIC ta2_real = z0_real - ta1_real / 2;
 122           const ATOMIC ta2_imag = z0_imag - ta1_imag / 2;
 123
 124           /* ta3 = (+/-) sin(pi/3)*(z2 - z4) */
 125           const ATOMIC ta3_real = ((int) sign) * tau * (z2_real - z4_real);
 126           const ATOMIC ta3_imag = ((int) sign) * tau * (z2_imag - z4_imag);
 127
 128           /* a0 = z0 + ta1 */
 129           const ATOMIC a0_real = z0_real + ta1_real;
 130           const ATOMIC a0_imag = z0_imag + ta1_imag;
 131
 132           /* a1 = ta2 + i ta3 */
 133           const ATOMIC a1_real = ta2_real - ta3_imag;
 134           const ATOMIC a1_imag = ta2_imag + ta3_real;
 135
 136           /* a2 = ta2 - i ta3 */
 137           const ATOMIC a2_real = ta2_real + ta3_imag;
 138           const ATOMIC a2_imag = ta2_imag - ta3_real;
 139
 140           /* tb1 = z5 + z1 */
 141           const ATOMIC tb1_real = z5_real + z1_real;
 142           const ATOMIC tb1_imag = z5_imag + z1_imag;
 143
 144           /* tb2 = z3 - tb1/2 */
 145           const ATOMIC tb2_real = z3_real - tb1_real / 2;
 146           const ATOMIC tb2_imag = z3_imag - tb1_imag / 2;
 147
 148           /* tb3 = (+/-) sin(pi/3)*(z5 - z1) */
 149           const ATOMIC tb3_real = ((int) sign) * tau * (z5_real - z1_real);
 150           const ATOMIC tb3_imag = ((int) sign) * tau * (z5_imag - z1_imag);
 151
 152           /* b0 = z3 + tb1 */
 153           const ATOMIC b0_real = z3_real + tb1_real;
 154           const ATOMIC b0_imag = z3_imag + tb1_imag;
 155
 156           /* b1 = tb2 + i tb3 */
 157           const ATOMIC b1_real = tb2_real - tb3_imag;
 158           const ATOMIC b1_imag = tb2_imag + tb3_real;
 159
 160           /* b2 = tb2 - i tb3 */
 161           const ATOMIC b2_real = tb2_real + tb3_imag;
 162           const ATOMIC b2_imag = tb2_imag - tb3_real;
 163
 164           /* x0 = a0 + b0 */
 165           const ATOMIC x0_real = a0_real + b0_real;
 166           const ATOMIC x0_imag = a0_imag + b0_imag;
 167
 168           /* x4 = a1 + b1 */
 169           const ATOMIC x4_real = a1_real + b1_real;
 170           const ATOMIC x4_imag = a1_imag + b1_imag;
 171
 172           /* x2 = a2 + b2 */
 173           const ATOMIC x2_real = a2_real + b2_real;
 174           const ATOMIC x2_imag = a2_imag + b2_imag;
 175
 176           /* x3 = a0 - b0 */
 177           const ATOMIC x3_real = a0_real - b0_real;
 178           const ATOMIC x3_imag = a0_imag - b0_imag;
 179
 180           /* x1 = a1 - b1 */
 181           const ATOMIC x1_real = a1_real - b1_real;
 182           const ATOMIC x1_imag = a1_imag - b1_imag;
 183
 184           /* x5 = a2 - b2 */
 185           const ATOMIC x5_real = a2_real - b2_real;
 186           const ATOMIC x5_imag = a2_imag - b2_imag;
 187
 188           /* apply twiddle factors */
 189
 190           /* to0 = 1 * x0 */
 191           REAL(out,ostride,j) = x0_real;
 192           IMAG(out,ostride,j) = x0_imag;
 193
 194           /* to1 = w1 * x1 */
 195           REAL(out,ostride,j+p_1) = w1_real * x1_real - w1_imag * x1_imag;
 196           IMAG(out,ostride,j+p_1) = w1_real * x1_imag + w1_imag * x1_real;
 197
 198           /* to2 = w2 * x2 */
 199           REAL(out,ostride,j+2*p_1) = w2_real * x2_real - w2_imag * x2_imag;
 200           IMAG(out,ostride,j+2*p_1) = w2_real * x2_imag + w2_imag * x2_real;
 201
 202           /* to3 = w3 * x3 */
 203           REAL(out,ostride,j+3*p_1) = w3_real * x3_real - w3_imag * x3_imag;
 204           IMAG(out,ostride,j+3*p_1) = w3_real * x3_imag + w3_imag * x3_real;
 205
 206           /* to4 = w4 * x4 */
 207           REAL(out,ostride,j+4*p_1) = w4_real * x4_real - w4_imag * x4_imag;
 208           IMAG(out,ostride,j+4*p_1) = w4_real * x4_imag + w4_imag * x4_real;
 209
 210           /* to5 = w5 * x5 */
 211           REAL(out,ostride,j+5*p_1) = w5_real * x5_real - w5_imag * x5_imag;
 212           IMAG(out,ostride,j+5*p_1) = w5_real * x5_imag + w5_imag * x5_real;
 213
 214           i++;
 215           j++;
 216         }
 217       j += jump;
 218     }
 219   return 0;
 220 }