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

   1 /* fft/dft_source.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 int
  21 FUNCTION(gsl_dft_complex,forward) (const BASE data[],
  22                                    const size_t stride, const size_t n,
  23                                    BASE result[])
  24 {
  25   gsl_fft_direction sign = gsl_fft_forward;
  26   int status = FUNCTION(gsl_dft_complex,transform) (data, stride, n, result, sign);
  27   return status;
  28 }
  29
  30 int
  31 FUNCTION(gsl_dft_complex,backward) (const BASE data[],
  32                                     const size_t stride, const size_t n,
  33                                     BASE result[])
  34 {
  35   gsl_fft_direction sign = gsl_fft_backward;
  36   int status = FUNCTION(gsl_dft_complex,transform) (data, stride, n, result, sign);
  37   return status;
  38 }
  39
  40
  41 int
  42 FUNCTION(gsl_dft_complex,inverse) (const BASE data[],
  43                                    const size_t stride, const size_t n,
  44                                    BASE result[])
  45 {
  46   gsl_fft_direction sign = gsl_fft_backward;
  47   int status = FUNCTION(gsl_dft_complex,transform) (data, stride, n, result, sign);
  48
  49   /* normalize inverse fft with 1/n */
  50
  51   {
  52     const ATOMIC norm = ONE / (ATOMIC)n;
  53     size_t i;
  54     for (i = 0; i < n; i++)
  55       {
  56         REAL(result,stride,i) *= norm;
  57         IMAG(result,stride,i) *= norm;
  58       }
  59   }
  60   return status;
  61 }
  62
  63 int
  64 FUNCTION(gsl_dft_complex,transform) (const BASE data[],
  65                                      const size_t stride, const size_t n,
  66                                      BASE result[],
  67                                      const gsl_fft_direction sign)
  68 {
  69
  70   size_t i, j, exponent;
  71
  72   const double d_theta = 2.0 * ((int) sign) * M_PI / (double) n;
  73
  74   /* FIXME: check that input length == output length and give error */
  75
  76   for (i = 0; i < n; i++)
  77     {
  78       ATOMIC sum_real = 0;
  79       ATOMIC sum_imag = 0;
  80
  81       exponent = 0;
  82
  83       for (j = 0; j < n; j++)
  84         {
  85           double theta = d_theta * (double) exponent;
  86           /* sum = exp(i theta) * data[j] */
  87
  88           ATOMIC w_real = (ATOMIC) cos (theta);
  89           ATOMIC w_imag = (ATOMIC) sin (theta);
  90
  91           ATOMIC data_real = REAL(data,stride,j);
  92           ATOMIC data_imag = IMAG(data,stride,j);
  93
  94           sum_real += w_real * data_real - w_imag * data_imag;
  95           sum_imag += w_real * data_imag + w_imag * data_real;
  96
  97           exponent = (exponent + i) % n;
  98         }
  99       REAL(result,stride,i) = sum_real;
 100       IMAG(result,stride,i) = sum_imag;
 101     }
 102
 103   return 0;
 104 }