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

   1 /* fft/bitreverse.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 #include <config.h>
  21 #include <gsl/gsl_fft.h>
  22
  23 #include "complex_internal.h"
  24 #include "bitreverse.h"
  25
  26 static int
  27 FUNCTION(fft_complex,bitreverse_order) (BASE data[],
  28                                         const size_t stride,
  29                                         const size_t n,
  30                                         size_t logn)
  31 {
  32   /* This is the Goldrader bit-reversal algorithm */
  33
  34   size_t i;
  35   size_t j = 0;
  36
  37   logn = 0 ; /* not needed for this algorithm */
  38
  39   for (i = 0; i < n - 1; i++)
  40     {
  41       size_t k = n / 2 ;
  42
  43       if (i < j)
  44         {
  45           const BASE tmp_real = REAL(data,stride,i);
  46           const BASE tmp_imag = IMAG(data,stride,i);
  47           REAL(data,stride,i) = REAL(data,stride,j);
  48           IMAG(data,stride,i) = IMAG(data,stride,j);
  49           REAL(data,stride,j) = tmp_real;
  50           IMAG(data,stride,j) = tmp_imag;
  51         }
  52
  53       while (k <= j)
  54         {
  55           j = j - k ;
  56           k = k / 2 ;
  57         }
  58
  59       j += k ;
  60     }
  61
  62   return 0;
  63 }
  64
  65
  66 static int
  67 FUNCTION(fft_real,bitreverse_order) (BASE data[],
  68                                 const size_t stride,
  69                                 const size_t n,
  70                                 size_t logn)
  71 {
  72   /* This is the Goldrader bit-reversal algorithm */
  73
  74   size_t i;
  75   size_t j = 0;
  76
  77   logn = 0 ; /* not needed for this algorithm */
  78
  79   for (i = 0; i < n - 1; i++)
  80     {
  81       size_t k = n / 2 ;
  82
  83       if (i < j)
  84         {
  85           const BASE tmp = VECTOR(data,stride,i);
  86           VECTOR(data,stride,i) = VECTOR(data,stride,j);
  87           VECTOR(data,stride,j) = tmp;
  88         }
  89
  90       while (k <= j)
  91         {
  92           j = j - k ;
  93           k = k / 2 ;
  94         }
  95
  96       j += k ;
  97     }
  98
  99   return 0;
 100 }
 101