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

   1 /* fft/c_pass_n.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_n) (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 factor,
  27                               const size_t product,
  28                               const size_t n,
  29                               const TYPE(gsl_complex) twiddle[])
  30 {
  31   size_t i = 0, j = 0;
  32   size_t k, k1;
  33
  34   const size_t m = n / factor;
  35   const size_t q = n / product;
  36   const size_t p_1 = product / factor;
  37   const size_t jump = (factor - 1) * p_1;
  38
  39   size_t e, e1;
  40
  41   for (i = 0; i < m; i++)
  42     {
  43       REAL(out,ostride,i) = REAL(in,istride,i);
  44       IMAG(out,ostride,i) = IMAG(in,istride,i);
  45     }
  46
  47   for (e = 1; e < (factor - 1) / 2 + 1; e++)
  48     {
  49       for (i = 0; i < m; i++)
  50         {
  51           const size_t idx = i + e * m;
  52           const size_t idxc = i + (factor - e) * m;
  53           REAL(out,ostride,idx) = REAL(in,istride,idx) + REAL(in,istride,idxc);
  54           IMAG(out,ostride,idx) = IMAG(in,istride,idx) + IMAG(in,istride,idxc);
  55           REAL(out,ostride,idxc) = REAL(in,istride,idx) - REAL(in,istride,idxc);
  56           IMAG(out,ostride,idxc) = IMAG(in,istride,idx) - IMAG(in,istride,idxc);
  57         }
  58     }
  59
  60   /* e = 0 */
  61
  62   for (i=0 ; i<m; i++)
  63     {
  64       REAL(in,istride,i) = REAL(out,ostride,i);
  65       IMAG(in,istride,i) = IMAG(out,ostride,i);
  66     }
  67
  68   for (e1 = 1; e1 < (factor - 1) / 2 + 1; e1++)
  69     {
  70       for (i = 0; i < m; i++)
  71         {
  72           REAL(in,istride,i) += REAL(out,ostride,i + e1*m) ;
  73           IMAG(in,istride,i) += IMAG(out,ostride,i + e1*m) ;
  74         }
  75     }
  76
  77   for (e = 1; e < (factor-1)/2 + 1; e++)
  78     {
  79       size_t idx = e*q ;
  80       const size_t idx_step = e * q ;
  81       ATOMIC w_real, w_imag ;
  82
  83       const size_t em = e * m ;
  84       const size_t ecm = (factor - e) * m ;
  85
  86       for (i = 0; i < m; i++)
  87         {
  88           REAL(in,istride,i+em) = REAL(out,ostride,i) ;
  89           IMAG(in,istride,i+em) = IMAG(out,ostride,i) ;
  90           REAL(in,istride,i+ecm) = REAL(out,ostride,i) ;
  91           IMAG(in,istride,i+ecm) = IMAG(out,ostride,i) ;
  92         }
  93
  94       for (e1 = 1; e1 < (factor - 1) / 2 + 1; e1++)
  95         {
  96           if (idx == 0) {
  97             w_real = 1 ;
  98             w_imag = 0 ;
  99           } else {
 100             if (sign == gsl_fft_forward) {
 101               w_real = GSL_REAL(twiddle[idx - 1]) ;
 102               w_imag = GSL_IMAG(twiddle[idx - 1]) ;
 103             } else {
 104               w_real = GSL_REAL(twiddle[idx - 1]) ;
 105               w_imag = -GSL_IMAG(twiddle[idx - 1]) ;
 106             }
 107           }
 108
 109           for (i = 0; i < m; i++)
 110             {
 111               const ATOMIC xp_real = REAL(out,ostride,i + e1 * m);
 112               const ATOMIC xp_imag = IMAG(out,ostride,i + e1 * m);
 113               const ATOMIC xm_real = REAL(out,ostride,i + (factor - e1) *m);
 114               const ATOMIC xm_imag = IMAG(out,ostride,i + (factor - e1) *m);
 115
 116               const ATOMIC ap = w_real * xp_real ;
 117               const ATOMIC am = w_imag * xm_imag ;
 118
 119               ATOMIC sum_real = ap - am;
 120               ATOMIC sumc_real = ap + am;
 121
 122               const ATOMIC bp = w_real * xp_imag ;
 123               const ATOMIC bm = w_imag * xm_real ;
 124
 125               ATOMIC sum_imag = bp + bm;
 126               ATOMIC sumc_imag = bp - bm;
 127
 128               REAL(in,istride,i + em) += sum_real;
 129               IMAG(in,istride,i + em) += sum_imag;
 130               REAL(in,istride,i + ecm) += sumc_real;
 131               IMAG(in,istride,i + ecm) += sumc_imag;
 132             }
 133           idx += idx_step ;
 134           idx %= factor * q ;
 135         }
 136     }
 137
 138   i = 0;
 139   j = 0;
 140
 141   /* k = 0 */
 142   for (k1 = 0; k1 < p_1; k1++)
 143     {
 144       REAL(out,ostride,k1) = REAL(in,istride,k1);
 145       IMAG(out,ostride,k1) = IMAG(in,istride,k1);
 146     }
 147
 148   for (e1 = 1; e1 < factor; e1++)
 149     {
 150       for (k1 = 0; k1 < p_1; k1++)
 151         {
 152           REAL(out,ostride,k1 + e1 * p_1) = REAL(in,istride,k1 + e1 * m) ;
 153           IMAG(out,ostride,k1 + e1 * p_1) = IMAG(in,istride,k1 + e1 * m) ;
 154         }
 155     }
 156
 157   i = p_1 ;
 158   j = product ;
 159
 160   for (k = 1; k < q; k++)
 161     {
 162       for (k1 = 0; k1 < p_1; k1++)
 163         {
 164           REAL(out,ostride,j) = REAL(in,istride,i);
 165           IMAG(out,ostride,j) = IMAG(in,istride,i);
 166           i++;
 167           j++;
 168         }
 169       j += jump;
 170     }
 171
 172   i = p_1 ;
 173   j = product ;
 174
 175   for (k = 1; k < q; k++)
 176     {
 177       for (k1 = 0; k1 < p_1; k1++)
 178         {
 179           for (e1 = 1; e1 < factor; e1++)
 180             {
 181               ATOMIC x_real = REAL(in, istride,i + e1 * m);
 182               ATOMIC x_imag = IMAG(in, istride,i + e1 * m);
 183
 184               ATOMIC w_real, w_imag ;
 185               if (sign == gsl_fft_forward) {
 186                 w_real = GSL_REAL(twiddle[(e1-1)*q + k-1]) ;
 187                 w_imag = GSL_IMAG(twiddle[(e1-1)*q + k-1]) ;
 188               } else {
 189                 w_real = GSL_REAL(twiddle[(e1-1)*q + k-1]) ;
 190                 w_imag = -GSL_IMAG(twiddle[(e1-1)*q + k-1]) ;
 191               }
 192
 193               REAL(out,ostride,j + e1 * p_1) = w_real * x_real - w_imag * x_imag;
 194               IMAG(out,ostride,j + e1 * p_1) = w_real * x_imag + w_imag * x_real;
 195             }
 196           i++;
 197           j++;
 198         }
 199       j += jump;
 200     }
 201
 202   return 0;
 203 }
 204