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

   1 /* fft/real_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 void
  21 FUNCTION(fft_real,pass_n) (const BASE in[],
  22                            const size_t istride,
  23                            BASE out[],
  24                            const size_t ostride,
  25                            const size_t factor,
  26                            const size_t product,
  27                            const size_t n,
  28                            const TYPE(gsl_complex) twiddle[])
  29 {
  30   size_t k, k1;
  31
  32   const size_t m = n / factor;
  33   const size_t q = n / product;
  34   const size_t product_1 = product / factor;
  35
  36   size_t e1, e2;
  37
  38   const double d_theta = 2.0 * M_PI / ((double) factor);
  39   const ATOMIC cos_d_theta = cos (d_theta);
  40   const ATOMIC sin_d_theta = sin (d_theta);
  41
  42   for (k1 = 0; k1 < q; k1++)
  43     {
  44       /* compute x = W(factor) z, for z real */
  45
  46       ATOMIC dw_real = 1.0, dw_imag = 0.0;
  47
  48       for (e1 = 0; e1 <= factor - e1; e1++)
  49         {
  50           ATOMIC sum_real = 0.0;
  51           ATOMIC sum_imag = 0.0;
  52
  53           ATOMIC w_real = 1.0, w_imag = 0.0;
  54
  55           if (e1 > 0)
  56             {
  57               ATOMIC tmp_real = dw_real * cos_d_theta + dw_imag * sin_d_theta;
  58               ATOMIC tmp_imag = -dw_real * sin_d_theta + dw_imag * cos_d_theta;
  59               dw_real = tmp_real;
  60               dw_imag = tmp_imag;
  61             }
  62
  63           for (e2 = 0; e2 < factor; e2++)
  64             {
  65               ATOMIC z_real = VECTOR(in,istride,k1 * product_1 + e2 * m);
  66
  67               if (e2 > 0)
  68                 {
  69                   ATOMIC tmp_real = dw_real * w_real - dw_imag * w_imag;
  70                   ATOMIC tmp_imag = dw_real * w_imag + dw_imag * w_real;
  71                   w_real = tmp_real;
  72                   w_imag = tmp_imag;
  73                 }
  74
  75               sum_real += w_real * z_real;
  76               sum_imag += w_imag * z_real;
  77
  78             }
  79           if (e1 == 0)
  80             {
  81               const size_t to0 = product * k1;
  82               VECTOR(out,ostride,to0) = sum_real;
  83             }
  84           else if (e1 < factor - e1)
  85             {
  86               const size_t to0 = k1 * product + 2 * e1 * product_1 - 1;
  87               VECTOR(out,ostride,to0) = sum_real;
  88               VECTOR(out,ostride,to0 + 1) = sum_imag;
  89             }
  90           else if (e1 == factor - e1)
  91             {
  92               const size_t to0 = k1 * product + 2 * e1 * product_1 - 1;
  93               VECTOR(out,ostride,to0) = sum_real;
  94             }
  95
  96         }
  97     }
  98
  99   if (product_1 == 1)
 100     return;
 101
 102   for (k = 1; k < (product_1 + 1) / 2; k++)
 103     {
 104       for (k1 = 0; k1 < q; k1++)
 105         {
 106
 107           ATOMIC dw_real = 1.0, dw_imag = 0.0;
 108
 109           for (e1 = 0; e1 < factor; e1++)
 110             {
 111               ATOMIC sum_real = 0.0, sum_imag = 0.0;
 112
 113               ATOMIC w_real = 1.0, w_imag = 0.0;
 114
 115               if (e1 > 0)
 116                 {
 117                   const ATOMIC tmp_real = dw_real * cos_d_theta + dw_imag * sin_d_theta;
 118                   const ATOMIC tmp_imag = -dw_real * sin_d_theta + dw_imag * cos_d_theta;
 119                   dw_real = tmp_real;
 120                   dw_imag = tmp_imag;
 121                 }
 122
 123               for (e2 = 0; e2 < factor; e2++)
 124                 {
 125
 126                   int tskip = (product_1 + 1) / 2 - 1;
 127                   const size_t from0 = k1 * product_1 + 2 * k + e2 * m - 1;
 128                   ATOMIC tw_real, tw_imag;
 129                   ATOMIC z_real, z_imag;
 130
 131                   if (e2 == 0)
 132                     {
 133                       tw_real = 1.0;
 134                       tw_imag = 0.0;
 135                     }
 136                   else
 137                     {
 138                       const size_t t_index = (k - 1) + (e2 - 1) * tskip;
 139                       tw_real = GSL_REAL(twiddle[t_index]);
 140                       tw_imag = -GSL_IMAG(twiddle[t_index]);
 141                     }
 142
 143                   {
 144                     const ATOMIC f0_real = VECTOR(in,istride,from0);
 145                     const ATOMIC f0_imag = VECTOR(in,istride,from0 + 1);
 146
 147                     z_real = tw_real * f0_real - tw_imag * f0_imag;
 148                     z_imag = tw_real * f0_imag + tw_imag * f0_real;
 149                   }
 150
 151                   if (e2 > 0)
 152                     {
 153                       const ATOMIC tmp_real = dw_real * w_real - dw_imag * w_imag;
 154                       const ATOMIC tmp_imag = dw_real * w_imag + dw_imag * w_real;
 155                       w_real = tmp_real;
 156                       w_imag = tmp_imag;
 157                     }
 158
 159                   sum_real += w_real * z_real - w_imag * z_imag;
 160                   sum_imag += w_real * z_imag + w_imag * z_real;
 161                 }
 162
 163               if (e1 < factor - e1)
 164                 {
 165                   const size_t to0 = k1 * product - 1 + 2 * e1 * product_1 + 2 * k;
 166                   VECTOR(out,ostride,to0) = sum_real;
 167                   VECTOR(out,ostride,to0 + 1) = sum_imag;
 168                 }
 169               else
 170                 {
 171                   const size_t to0 = k1 * product - 1 + 2 * (factor - e1) * product_1 - 2 * k;
 172                   VECTOR(out,ostride,to0) = sum_real;
 173                   VECTOR(out,ostride,to0 + 1) = -sum_imag;
 174                 }
 175
 176             }
 177         }
 178     }
 179
 180
 181   if (product_1 % 2 == 1)
 182     return;
 183
 184   {
 185     double tw_arg = M_PI / ((double) factor);
 186     ATOMIC cos_tw_arg = cos (tw_arg);
 187     ATOMIC sin_tw_arg = -sin (tw_arg);
 188
 189     for (k1 = 0; k1 < q; k1++)
 190       {
 191         ATOMIC dw_real = 1.0, dw_imag = 0.0;
 192
 193         for (e1 = 0; e1 < factor; e1++)
 194           {
 195             ATOMIC z_real, z_imag;
 196
 197             ATOMIC sum_real = 0.0;
 198             ATOMIC sum_imag = 0.0;
 199
 200             ATOMIC w_real = 1.0, w_imag = 0.0;
 201             ATOMIC tw_real = 1.0, tw_imag = 0.0;
 202
 203             if (e1 > 0)
 204               {
 205                 ATOMIC t_real = dw_real * cos_d_theta + dw_imag * sin_d_theta;
 206                 ATOMIC t_imag = -dw_real * sin_d_theta + dw_imag * cos_d_theta;
 207                 dw_real = t_real;
 208                 dw_imag = t_imag;
 209               }
 210
 211             for (e2 = 0; e2 < factor; e2++)
 212               {
 213
 214                 if (e2 > 0)
 215                   {
 216                     ATOMIC tmp_real = tw_real * cos_tw_arg - tw_imag * sin_tw_arg;
 217                     ATOMIC tmp_imag = tw_real * sin_tw_arg + tw_imag * cos_tw_arg;
 218                     tw_real = tmp_real;
 219                     tw_imag = tmp_imag;
 220                   }
 221
 222                 if (e2 > 0)
 223                   {
 224                     ATOMIC tmp_real = dw_real * w_real - dw_imag * w_imag;
 225                     ATOMIC tmp_imag = dw_real * w_imag + dw_imag * w_real;
 226                     w_real = tmp_real;
 227                     w_imag = tmp_imag;
 228                   }
 229
 230
 231                 {
 232                   const size_t from0 = k1 * product_1 + 2 * k + e2 * m - 1;
 233                   const ATOMIC f0_real = VECTOR(in,istride,from0);
 234                   z_real = tw_real * f0_real;
 235                   z_imag = tw_imag * f0_real;
 236                 }
 237
 238                 sum_real += w_real * z_real - w_imag * z_imag;
 239                 sum_imag += w_real * z_imag + w_imag * z_real;
 240               }
 241
 242             if (e1 + 1 < factor - e1)
 243               {
 244                 const size_t to0 = k1 * product - 1 + 2 * e1 * product_1 + 2 * k;
 245                 VECTOR(out,ostride,to0) = sum_real;
 246                 VECTOR(out,ostride,to0 + 1) = sum_imag;
 247               }
 248             else if (e1 + 1 == factor - e1)
 249               {
 250                 const size_t to0 = k1 * product - 1 + 2 * e1 * product_1 + 2 * k;
 251                 VECTOR(out,ostride,to0) = sum_real;
 252               }
 253             else
 254               {
 255                 const size_t to0 = k1 * product - 1 + 2 * (factor - e1) * product_1 - 2 * k;
 256                 VECTOR(out,ostride,to0) = sum_real;
 257                 VECTOR(out,ostride,to0 + 1) = -sum_imag;
 258               }
 259
 260           }
 261       }
 262   }
 263   return;
 264 }