Olga Sergijenko / gammalib

Commit 768912b8cc23b7776465d72b9d04ddb3de92515a

Authored by Jürgen Knödlseder
1 parent 32eaa5e9

Add GFftWavetable and GFft classes 

The GFft class implements a one-dimensional Fast Fourier transformation. The GFftWavetable class is a support class for GFft and holds the trigonometric coefficients for a factorisation.
Inline Side-by-side
Showing 21 changed files with 3390 additions and 27 deletions
ChangeLog
  View file @768912b
1   -2016-10-08
  1 +2016-10-16
2 2  
3 3 * Version 1.2.0 released
4 4 =========================
5 5  
  6 + Add GFft and GFftWavetable classes (#1731)
6 7 Add GNdarray class (#1768)
7 8 Increase precision of gammalib::str() floating point functions
8 9 Add GFits::publish() methods (#1231)
... ...
1 1 New Features and Important Changes in GammaLib 1.2.0
2 2  
3   -8 October 2016
  3 +16 October 2016
4 4  
5 5  
6 6 1. Introduction
... ... @@ -17,6 +17,8 @@ The following classes have been added:
17 17 - GModelSpectralMultiplicative
18 18 - GVOTable
19 19 - GNdarray
  20 +- GFft
  21 +- GFftWavetable
20 22  
21 23 The following classes have been removed:
22 24 - ...
... ... @@ -159,6 +161,11 @@ product of spectral model components as a new spectral model (#1861).
159 161 The GNdarray class has been added that implements a n-dimensional double
160 162 precision floating point array (#1768).
161 163  
  164 +The GFft class has been added that performs a Fast Fourier Transformation
  165 +of a n-dimensional array. The GFftWavetable class is a support class for
  166 +GFft that provides the trigonometric coefficients for a factorisation. For
  167 +the moment, GFft does only support 1-dimensional arrays (#1731).
  168 +
162 169  
163 170 10. Observation module
164 171 ----------------------
... ...
include/GFft.hpp 0 → 100644
  View file @768912b
  1 +/***************************************************************************
  2 + * GFft.hpp - Fast Fourier transformation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GFft.hpp
  23 + * @brief Fast Fourier transformation class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GFFT_HPP
  28 +#define GFFT_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <vector>
  33 +#include <complex>
  34 +#include "GBase.hpp"
  35 +#include "GFftWavetable.hpp"
  36 +
  37 +/* __ Forward declarations _______________________________________________ */
  38 +class GNdarray;
  39 +//class GFftWavetable;
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GFft
  44 + *
  45 + * @brief Fast Fourier Transformation class
  46 + *
  47 + * This class implements a Fast Fourier Transformation of an n-dimensional
  48 + * double precision floating point array.
  49 + *
  50 + * The class implementation is based on the Fast Fourier Transformation
  51 + * functions that are provided by the GNU Scientific Library (GSL)
  52 + * (version 2.2.1). The mixed-radix routines for complex numbers have been
  53 + * adopted that work for any array lengths. The mixed-radix routines are a
  54 + * reimplementation of the FFTPACK library of Paul Swarztrauber.
  55 + *
  56 + * The GSL is documented at https://www.gnu.org/software/gsl/manual
  57 + ***************************************************************************/
  58 +class GFft : public GBase {
  59 +
  60 +public:
  61 + // Constructors and destructors
  62 + GFft(void);
  63 + explicit GFft(const GNdarray& array);
  64 + GFft(const GFft& fft);
  65 + virtual ~GFft(void);
  66 +
  67 + // FFT access operators
  68 + std::complex<double>& operator()(const int& ix);
  69 + const std::complex<double>& operator()(const int& ix) const;
  70 +
  71 + // Operators
  72 + GFft& operator=(const GFft& fft);
  73 + GFft& operator+=(const GFft& fft);
  74 + GFft& operator-=(const GFft& fft);
  75 + GFft& operator*=(const GFft& fft);
  76 + GFft& operator/=(const GFft& fft);
  77 + GFft operator-(void) const;
  78 +
  79 + // Methods
  80 + void clear(void);
  81 + GFft* clone(void) const;
  82 + std::string classname(void) const;
  83 + int dim() const;
  84 + int size(void) const;
  85 + const std::vector<int>& shape(void) const;
  86 + const std::vector<int>& strides(void) const;
  87 + void forward(const GNdarray& array);
  88 + GNdarray backward(void) const;
  89 + std::string print(const GChatter& chatter = NORMAL) const;
  90 +
  91 +protected:
  92 + // Protected methods
  93 + void init_members(void);
  94 + void copy_members(const GFft& fft);
  95 + void free_members(void);
  96 + void set_data(const GNdarray& array);
  97 +
  98 + // Low-level FFT methods
  99 + void transform(std::vector<std::complex<double> >& data,
  100 + const int& stride,
  101 + const int& n,
  102 + const GFftWavetable& wavetable,
  103 + const bool& forward = true);
  104 + void factor2(const std::vector<std::complex<double> >& in,
  105 + const int& istride,
  106 + std::vector<std::complex<double> >& out,
  107 + const int& ostride,
  108 + const GFftWavetable& wavetable,
  109 + const int& sign,
  110 + const int& product,
  111 + const int& n,
  112 + const int& index);
  113 + void factor3(const std::vector<std::complex<double> >& in,
  114 + const int& istride,
  115 + std::vector<std::complex<double> >& out,
  116 + const int& ostride,
  117 + const GFftWavetable& wavetable,
  118 + const int& sign,
  119 + const int& product,
  120 + const int& n,
  121 + const int& index);
  122 + void factor4(const std::vector<std::complex<double> >& in,
  123 + const int& istride,
  124 + std::vector<std::complex<double> >& out,
  125 + const int& ostride,
  126 + const GFftWavetable& wavetable,
  127 + const int& sign,
  128 + const int& product,
  129 + const int& n,
  130 + const int& index);
  131 + void factor5(const std::vector<std::complex<double> >& in,
  132 + const int& istride,
  133 + std::vector<std::complex<double> >& out,
  134 + const int& ostride,
  135 + const GFftWavetable& wavetable,
  136 + const int& sign,
  137 + const int& product,
  138 + const int& n,
  139 + const int& index);
  140 + void factor6(const std::vector<std::complex<double> >& in,
  141 + const int& istride,
  142 + std::vector<std::complex<double> >& out,
  143 + const int& ostride,
  144 + const GFftWavetable& wavetable,
  145 + const int& sign,
  146 + const int& product,
  147 + const int& n,
  148 + const int& index);
  149 + void factor7(const std::vector<std::complex<double> >& in,
  150 + const int& istride,
  151 + std::vector<std::complex<double> >& out,
  152 + const int& ostride,
  153 + const GFftWavetable& wavetable,
  154 + const int& sign,
  155 + const int& product,
  156 + const int& n,
  157 + const int& index);
  158 + void factorn(std::vector<std::complex<double> >& in,
  159 + const int& istride,
  160 + std::vector<std::complex<double> >& out,
  161 + const int& ostride,
  162 + const GFftWavetable& wavetable,
  163 + const int& sign,
  164 + const int& factor,
  165 + const int& product,
  166 + const int& n,
  167 + const int& index);
  168 +
  169 + // Protected members
  170 + std::vector<int> m_shape; //!< Array dimensions
  171 + std::vector<int> m_strides; //!< Steps in each dimension
  172 + std::vector<std::complex<double> > m_data; //!< Array data
  173 + std::vector<GFftWavetable> m_wavetable; //!< Trigonometric coefficients
  174 +};
  175 +
  176 +
  177 +/***********************************************************************//**
  178 + * @brief Return class name
  179 + *
  180 + * @return String containing the class name ("GFft").
  181 + ***************************************************************************/
  182 +inline
  183 +std::string GFft::classname(void) const
  184 +{
  185 + return ("GFft");
  186 +}
  187 +
  188 +
  189 +/***********************************************************************//**
  190 + * @brief 1-dimensional FFT element access operator
  191 + *
  192 + * @param[in] ix Element index [0,...,shape(0)-1].
  193 + * @return Reference to FFT element.
  194 + ***************************************************************************/
  195 +inline
  196 +std::complex<double>& GFft::operator()(const int& ix)
  197 +{
  198 + // Return array element
  199 + return m_data[ix];
  200 +}
  201 +
  202 +
  203 +/***********************************************************************//**
  204 + * @brief 1-dimensional FFT element access operator (const variant)
  205 + *
  206 + * @param[in] ix Element index [0,...,shape(0)-1].
  207 + * @return Const reference to FFT element.
  208 + ***************************************************************************/
  209 +inline
  210 +const std::complex<double>& GFft::operator()(const int& ix) const
  211 +{
  212 + // Return array element
  213 + return m_data[ix];
  214 +}
  215 +
  216 +
  217 +/***********************************************************************//**
  218 + * @brief Return dimension of Fast Fourier Transformation
  219 + *
  220 + * @return Dimension of Fast Fourier Transformation.
  221 + *
  222 + * Returns the dimension of the Fast Fourier Transformation.
  223 + ***************************************************************************/
  224 +inline
  225 +int GFft::dim() const
  226 +{
  227 + return m_shape.size();
  228 +}
  229 +
  230 +
  231 +/***********************************************************************//**
  232 + * @brief Return number of elements in array
  233 + *
  234 + * @return Number of elements in array
  235 + *
  236 + * Returns the number of elements in the array.
  237 + ***************************************************************************/
  238 +inline
  239 +int GFft::size() const
  240 +{
  241 + return m_data.size();
  242 +}
  243 +
  244 +
  245 +/***********************************************************************//**
  246 + * @brief Return shape of array
  247 + *
  248 + * @return Shape of array
  249 + *
  250 + * Returns the shape of the array.
  251 + ***************************************************************************/
  252 +inline
  253 +const std::vector<int>& GFft::shape(void) const
  254 +{
  255 + return m_shape;
  256 +}
  257 +
  258 +
  259 +/***********************************************************************//**
  260 + * @brief Return strides of array
  261 + *
  262 + * @return Strides of array
  263 + *
  264 + * Returns the strides of the array.
  265 + ***************************************************************************/
  266 +inline
  267 +const std::vector<int>& GFft::strides(void) const
  268 +{
  269 + return m_strides;
  270 +}
  271 +
  272 +
  273 +/***********************************************************************//**
  274 + * @brief Return sum of two Fast Fourier Transformations
  275 + *
  276 + * @param[in] a First Fast Fourier Transformation.
  277 + * @param[in] b Second Fast Fourier Transformation.
  278 + * @return Sum of the Fast Fourier Transformations @p a and @p b.
  279 + *
  280 + * Returns the sum of the Fast Fourier Transformations @p a and @p b.
  281 + ***************************************************************************/
  282 +inline
  283 +GFft operator+(const GFft& a, const GFft& b)
  284 +{
  285 + GFft result = a;
  286 + result += b;
  287 + return result;
  288 +}
  289 +
  290 +
  291 +/***********************************************************************//**
  292 + * @brief Return difference of two Fast Fourier Transformations
  293 + *
  294 + * @param[in] a First Fast Fourier Transformation.
  295 + * @param[in] b Second Fast Fourier Transformation.
  296 + * @return Difference of the Fast Fourier Transformations @p a and @p b.
  297 + *
  298 + * Returns the difference of the Fast Fourier Transformations @p a and @p b.
  299 + ***************************************************************************/
  300 +inline
  301 +GFft operator-(const GFft& a, const GFft& b)
  302 +{
  303 + GFft result = a;
  304 + result -= b;
  305 + return result;
  306 +}
  307 +
  308 +
  309 +/***********************************************************************//**
  310 + * @brief Return product of two Fast Fourier Transformations
  311 + *
  312 + * @param[in] a First Fast Fourier Transformation.
  313 + * @param[in] b Second Fast Fourier Transformation.
  314 + * @return Product of the Fast Fourier Transformations @p a and @p b.
  315 + *
  316 + * Returns the product of the Fast Fourier Transformations @p a and @p b.
  317 + ***************************************************************************/
  318 +inline
  319 +GFft operator*(const GFft& a, const GFft& b)
  320 +{
  321 + GFft result = a;
  322 + result *= b;
  323 + return result;
  324 +}
  325 +
  326 +
  327 +/***********************************************************************//**
  328 + * @brief Return quotient of two Fast Fourier Transformations
  329 + *
  330 + * @param[in] a First Fast Fourier Transformation.
  331 + * @param[in] b Second Fast Fourier Transformation.
  332 + * @return Quotient of the Fast Fourier Transformations @p a and @p b.
  333 + *
  334 + * Returns the quotient of the Fast Fourier Transformations @p a and @p b.
  335 + ***************************************************************************/
  336 +inline
  337 +GFft operator/(const GFft& a, const GFft& b)
  338 +{
  339 + GFft result = a;
  340 + result /= b;
  341 + return result;
  342 +}
  343 +
  344 +#endif /* GFFT_HPP */
... ...
include/GFftWavetable.hpp 0 → 100644
  View file @768912b
  1 +/***************************************************************************
  2 + * GFftWavetable.hpp - Lookup table class for Fast Fourier transformation *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GFftWavetable.hpp
  23 + * @brief Lookup table class interface definition for Fast Fourier transformation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GFFTWAVETABLE_HPP
  28 +#define GFFTWAVETABLE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <vector>
  33 +#include <complex>
  34 +#include "GBase.hpp"
  35 +
  36 +/* __ Forward declarations _______________________________________________ */
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class GFftWavetable
  41 + *
  42 + * @brief Lookup table class for Fast Fourier Transformation
  43 + ***************************************************************************/
  44 +class GFftWavetable : public GBase {
  45 +
  46 +public:
  47 + // Constructors and destructors
  48 + GFftWavetable(void);
  49 + explicit GFftWavetable(const int& size);
  50 + GFftWavetable(const GFftWavetable& wavetable);
  51 + virtual ~GFftWavetable(void);
  52 +
  53 + // Operators
  54 + GFftWavetable& operator=(const GFftWavetable& wavetable);
  55 + std::complex<double>& operator[](const int& index);
  56 + const std::complex<double>& operator[](const int& index) const;
  57 +
  58 + // Methods
  59 + void clear(void);
  60 + GFftWavetable* clone(void) const;
  61 + std::string classname(void) const;
  62 + int size(void) const;
  63 + const int& index(const int& factor) const;
  64 + int factors(void) const;
  65 + int factor(const int& index) const;
  66 + std::string print(const GChatter& chatter = NORMAL) const;
  67 +
  68 +protected:
  69 + // Protected methods
  70 + void init_members(void);
  71 + void copy_members(const GFftWavetable& wavetable);
  72 + void free_members(void);
  73 + void set_members(const int& n);
  74 + void set_factors(const int& n);
  75 +
  76 + // Protected members
  77 + std::vector<int> m_factors; //!< Wavetable factors
  78 + std::vector<int> m_twiddle; //!< Start index of factors
  79 + std::vector<std::complex<double> > m_trig; //!< Trigonometric coefficients
  80 +};
  81 +
  82 +
  83 +/***********************************************************************//**
  84 + * @brief Return class name
  85 + *
  86 + * @return String containing the class name ("GFftWavetable").
  87 + ***************************************************************************/
  88 +inline
  89 +std::string GFftWavetable::classname(void) const
  90 +{
  91 + return ("GFftWavetable");
  92 +}
  93 +
  94 +
  95 +/***********************************************************************//**
  96 + * @brief Return reference to trigonometric coefficient
  97 + *
  98 + * @param[in] index Trigonometric coefficient index [0,...,size()-1].
  99 + *
  100 + * Returns a reference to the trigonometric coefficient with the specified
  101 + * @p index.
  102 + ***************************************************************************/
  103 +inline
  104 +std::complex<double>& GFftWavetable::operator[](const int& index)
  105 +{
  106 + return (m_trig[index]);
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Return reference to trigonometric coefficient (const version)
  112 + *
  113 + * @param[in] index Trigonometric coefficient index [0,...,size()-1].
  114 + *
  115 + * Returns a const reference to the trigonometric coefficient with the
  116 + * specified @p index.
  117 + ***************************************************************************/
  118 +inline
  119 +const std::complex<double>& GFftWavetable::operator[](const int& index) const
  120 +{
  121 + return (m_trig[index]);
  122 +}
  123 +
  124 +
  125 +/***********************************************************************//**
  126 + * @brief Return start index for a given factor
  127 + *
  128 + * @param[in] factor Factor index [0,...,factors()-1].
  129 + *
  130 + * Returns the index of the first trigonometric coefficient for a given
  131 + * @p factor.
  132 + ***************************************************************************/
  133 +inline
  134 +const int& GFftWavetable::index(const int& factor) const
  135 +{
  136 + return (m_twiddle[factor]);
  137 +}
  138 +
  139 +
  140 +/***********************************************************************//**
  141 + * @brief Return number of trigonometric coefficients
  142 + *
  143 + * @return Number of trigonometric coefficients.
  144 + ***************************************************************************/
  145 +inline
  146 +int GFftWavetable::size(void) const
  147 +{
  148 + return (m_trig.size());
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
  153 + * @brief Return number of factorisation factors
  154 + *
  155 + * @return Number of factorisation factors.
  156 + ***************************************************************************/
  157 +inline
  158 +int GFftWavetable::factors(void) const
  159 +{
  160 + return (m_factors.size());
  161 +}
  162 +
  163 +#endif /* GFFTWAVETABLE_HPP */
... ...
include/GNdarray.hpp
  View file @768912b
... ... @@ -38,7 +38,7 @@
38 38 *
39 39 * @brief N-dimensional array class
40 40 *
41   - * This class implement a n-dimensional double precision floating point
  41 + * This class implements a n-dimensional double precision floating point
42 42 * array.
43 43 ***************************************************************************/
44 44 class GNdarray : public GBase {
... ... @@ -105,6 +105,7 @@ public:
105 105 int dim() const;
106 106 int size(void) const;
107 107 const std::vector<int>& shape(void) const;
  108 + const std::vector<int>& strides(void) const;
108 109 double& at(const int& ix);
109 110 double& at(const int& ix, const int& iy);
110 111 double& at(const int& ix, const int& iy, const int& iz);
... ... @@ -115,9 +116,7 @@ public:
115 116 const double& at(const std::vector<int>& i) const;
116 117 std::string print(const GChatter& chatter = NORMAL) const;
117 118  
118   -#ifndef SWIG
119 119 protected:
120   -#endif
121 120 // Protected methods
122 121 void init_members(void);
123 122 void copy_members(const GNdarray& array);
... ... @@ -306,6 +305,20 @@ const std::vector&lt;int&gt;&amp; GNdarray::shape(void) const
306 305  
307 306  
308 307 /***********************************************************************//**
  308 + * @brief Return strides of array
  309 + *
  310 + * @return Strides of array
  311 + *
  312 + * Returns the strides of the array.
  313 + ***************************************************************************/
  314 +inline
  315 +const std::vector<int>& GNdarray::strides(void) const
  316 +{
  317 + return m_strides;
  318 +}
  319 +
  320 +
  321 +/***********************************************************************//**
309 322 * @brief 1-dimensional array element access with range checking
310 323 *
311 324 * @param[in] ix Element index [0,...,shape(0)-1].
... ...
include/GTestSuite.hpp
  View file @768912b
... ... @@ -34,6 +34,7 @@ class GTestSuite;
34 34 typedef void (GTestSuite::*pfunction)(void);
35 35  
36 36 /* __ Includes ___________________________________________________________ */
  37 +#include <complex>
37 38 #include <vector>
38 39 #include <string>
39 40 #include "GLog.hpp"
... ... @@ -92,6 +93,15 @@ public:
92 93 const double& eps,
93 94 const std::string& name = "",
94 95 const std::string& message = "");
  96 + void test_value(const std::complex<double>& value,
  97 + const std::complex<double>& expected,
  98 + const std::string& name = "",
  99 + const std::string& message = "");
  100 + void test_value(const std::complex<double>& value,
  101 + const std::complex<double>& expected,
  102 + const double& eps,
  103 + const std::string& name = "",
  104 + const std::string& message = "");
95 105 void test_value(const std::string& value,
96 106 const std::string& expected,
97 107 const std::string& name = "",
... ...
include/GTools.hpp
  View file @768912b
... ... @@ -30,6 +30,7 @@
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <vector>
32 32 #include <string>
  33 +#include <complex>
33 34 #include <cmath>
34 35 #include <cfloat>
35 36  
... ... @@ -69,6 +70,8 @@ namespace gammalib {
69 70 std::string str(const long long int& value);
70 71 std::string str(const float& value, const int& precision = 0);
71 72 std::string str(const double& value, const int& precision = 0);
  73 + std::string str(const std::complex<double>& value,
  74 + const int& precision = 0);
72 75 char* tochar(const std::string& arg);
73 76 short toshort(const std::string& arg);
74 77 unsigned short toushort(const std::string& arg);
... ...
include/GammaLib.hpp
  View file @768912b
... ... @@ -64,6 +64,8 @@
64 64 #include "GFunction.hpp"
65 65 #include "GMath.hpp"
66 66 #include "GNdarray.hpp"
  67 +#include "GFft.hpp"
  68 +#include "GFftWavetable.hpp"
67 69  
68 70 /* __ FITS module ________________________________________________________ */
69 71 #include "GFits.hpp"
... ...
include/Makefile.am
  View file @768912b
... ... @@ -45,6 +45,8 @@ pkginclude_HEADERS = GammaLib.hpp \
45 45 GFunction.hpp \
46 46 GMath.hpp \
47 47 GNdarray.hpp \
  48 + GFft.hpp \
  49 + GFftWavetable.hpp \
48 50 GFits.hpp \
49 51 GFitsHDU.hpp \
50 52 GFitsHeader.hpp \
... ...
pyext/GFft.i 0 → 100755
  View file @768912b
  1 +/***************************************************************************
  2 + * GFft.i - Fast Fourier transformation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GFft.i
  23 + * @brief Fast Fourier transformation class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "GFft.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GFft
  35 + *
  36 + * @brief Fast Fourier Transformation class
  37 + ***************************************************************************/
  38 +class GFft : public GBase {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + GFft(void);
  43 + explicit GFft(const GNdarray& array);
  44 + GFft(const GFft& fft);
  45 + virtual ~GFft(void);
  46 +
  47 + // Operators
  48 + GFft& operator+=(const GFft& fft);
  49 + GFft& operator-=(const GFft& fft);
  50 + GFft& operator*=(const GFft& fft);
  51 + GFft& operator/=(const GFft& fft);
  52 + GFft operator-(void) const;
  53 +
  54 + // Methods
  55 + void clear(void);
  56 + GFft* clone(void) const;
  57 + std::string classname(void) const;
  58 + int dim() const;
  59 + int size(void) const;
  60 + const std::vector<int>& shape(void) const;
  61 + const std::vector<int>& strides(void) const;
  62 + void forward(const GNdarray& array);
  63 + GNdarray backward(void) const;
  64 +};
  65 +
  66 +
  67 +/***********************************************************************//**
  68 + * @brief GFft class extension
  69 + ***************************************************************************/
  70 +%extend GFft {
  71 + std::complex<double> __getitem__(int GTuple[]) {
  72 + // Handle first the special case of a single index. In Python we want
  73 + // to support accessing a n-dimensional array with a single index so
  74 + // that we can iterate over the array
  75 + if (GTuple[0] == 1) {
  76 + if (GTuple[1] >= 0 && GTuple[1] < self->size()) {
  77 + return ((*self)(GTuple[1]));
  78 + }
  79 + else {
  80 + throw GException::out_of_range("__getitem__(int)", "Array "
  81 + "index", GTuple[1], self->size());
  82 + }
  83 + }
  84 + if (GTuple[0] != self->dim()) {
  85 + throw GException::invalid_value("__getitem__(int)", "Invalid "
  86 + "access of "+gammalib::str(self->dim())+"-dimensional "
  87 + "array with "+gammalib::str(GTuple[0])+"-dimensional "
  88 + "access operator.");
  89 + }
  90 + int index = 0;
  91 + for (int i = 0; i < self->shape().size(); ++i) {
  92 + if (GTuple[i+1] < 0 || GTuple[i+1] >= self->shape()[i]) {
  93 + throw GException::out_of_range("__getitem__(int)", "Dimension "+
  94 + gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
  95 + }
  96 + index += self->strides()[i] * GTuple[i+1];
  97 + }
  98 + return ((*self)(index));
  99 + }
  100 + void __setitem__(int GTuple[], const std::complex<double>& value) {
  101 + if (GTuple[0] != self->dim()) {
  102 + throw GException::invalid_value("__setitem__(int)", "Invalid "
  103 + "access of "+gammalib::str(self->dim())+"-dimensional "
  104 + "array with "+gammalib::str(GTuple[0])+"-dimensional "
  105 + "access operator.");
  106 + }
  107 + int index = 0;
  108 + for (int i = 0; i < self->shape().size(); ++i) {
  109 + if (GTuple[i+1] < 0 || GTuple[i+1] >= self->shape()[i]) {
  110 + throw GException::out_of_range("__setitem__(int)", "Dimension "+
  111 + gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
  112 + }
  113 + index += self->strides()[i] * GTuple[i+1];
  114 + }
  115 + (*self)(index) = value;
  116 + }
  117 + GFft copy() {
  118 + return (*self);
  119 + }
  120 +};
... ...
pyext/GFftWavetable.i 0 → 100755
  View file @768912b
  1 +/***************************************************************************
  2 + * GFftWavetable.i - Lookup table class for Fast Fourier transformation *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GFftWavetable.i
  23 + * @brief Lookup table class interface definition for Fast Fourier transformation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "GFftWavetable.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GFftWavetable
  35 + *
  36 + * @brief Lookup table class for Fast Fourier Transformation
  37 + ***************************************************************************/
  38 +class GFftWavetable : public GBase {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + GFftWavetable(void);
  43 + explicit GFftWavetable(const int& size);
  44 + GFftWavetable(const GFftWavetable& wavetable);
  45 + virtual ~GFftWavetable(void);
  46 +
  47 + // Methods
  48 + void clear(void);
  49 + GFftWavetable* clone(void) const;
  50 + std::string classname(void) const;
  51 + int size(void) const;
  52 + int factors(void) const;
  53 + int factor(const int& index) const;
  54 +};
  55 +
  56 +
  57 +/***********************************************************************//**
  58 + * @brief GFftWavetable class extension
  59 + ***************************************************************************/
  60 +%extend GFftWavetable {
  61 + GFftWavetable copy() {
  62 + return (*self);
  63 + }
  64 +};
... ...
pyext/GNdarray.i
  View file @768912b
... ... @@ -26,7 +26,6 @@
26 26  
27 27 %{
28 28 /* Put headers and other declarations here that are needed for compilation */
29   -#define SWIG
30 29 #include "GNdarray.hpp"
31 30 %}
32 31  
... ... @@ -68,6 +67,7 @@ public:
68 67 int dim() const;
69 68 int size(void) const;
70 69 const std::vector<int>& shape(void) const;
  70 + const std::vector<int>& strides(void) const;
71 71 };
72 72  
73 73  
... ... @@ -81,7 +81,7 @@ public:
81 81 // that we can iterate over the array
82 82 if (GTuple[0] == 1) {
83 83 if (GTuple[1] >= 0 && GTuple[1] < self->size()) {
84   - return (self->m_data[GTuple[1]]);
  84 + return ((*self)(GTuple[1]));
85 85 }
86 86 else {
87 87 throw GException::out_of_range("__getitem__(int)", "Array "
... ... @@ -100,9 +100,9 @@ public:
100 100 throw GException::out_of_range("__getitem__(int)", "Dimension "+
101 101 gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
102 102 }
103   - index += self->m_strides[i] * GTuple[i+1];
  103 + index += self->strides()[i] * GTuple[i+1];
104 104 }
105   - return (self->m_data[index]);
  105 + return ((*self)(index));
106 106 }
107 107 void __setitem__(int GTuple[], const double& value) {
108 108 if (GTuple[0] != self->dim()) {
... ... @@ -117,9 +117,9 @@ public:
117 117 throw GException::out_of_range("__setitem__(int)", "Dimension "+
118 118 gammalib::str(i)+" index", GTuple[i+1], self->shape()[i]);
119 119 }
120   - index += self->m_strides[i] * GTuple[i+1];
  120 + index += self->strides()[i] * GTuple[i+1];
121 121 }
122   - self->m_data[index] = value;
  122 + (*self)(index) = value;
123 123 }
124 124 GNdarray __add__(const GNdarray &a) {
125 125 return (*self) + a;
... ...
pyext/gammalib/numerics.i
  View file @768912b
... ... @@ -39,7 +39,9 @@
39 39 /* __ Include standard typemaps for vectors and strings __________________ */
40 40 %include stl.i
41 41 %include std_vector.i
  42 +%include std_complex.i
42 43 %template(vectori) std::vector<int>;
  44 +%template(vectorc) std::vector<std::complex<double> >;
43 45  
44 46 /* __ Include GammaLib typemaps __________________________________________ */
45 47 %include GTypemaps.i
... ... @@ -57,3 +59,5 @@
57 59 %include "GIntegral.i"
58 60 %include "GMath.i"
59 61 %include "GNdarray.i"
  62 +%include "GFft.i"
  63 +%include "GFftWavetable.i"
... ...
src/numerics/GFft.cpp 0 → 100644
  View file @768912b
  1 +/***************************************************************************
  2 + * GFft.hpp - Fast Fourier transformation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GFft.cpp
  23 + * @brief Fast Fourier transformation class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GMath.hpp"
  33 +#include "GException.hpp"
  34 +#include "GFft.hpp"
  35 +#include "GNdarray.hpp"
  36 +
  37 +/* __ Method name definitions ____________________________________________ */
  38 +
  39 +
  40 +/*==========================================================================
  41 + = =
  42 + = Constructors/destructors =
  43 + = =
  44 + ==========================================================================*/
  45 +
  46 +/***********************************************************************//**
  47 + * @brief Void constructor
  48 + ***************************************************************************/
  49 +GFft::GFft(void)
  50 +{
  51 + // Initialise class members
  52 + init_members();
  53 +
  54 + // Return
  55 + return;
  56 +}
  57 +
  58 +
  59 +/***********************************************************************//**
  60 + * @brief N-dimensional array constructor
  61 + *
  62 + * @param[in] array N-dimensional array.
  63 + *
  64 + * Constructs a Fast Fourier Transformation from a n-dimensional array.
  65 + ***************************************************************************/
  66 +GFft::GFft(const GNdarray& array)
  67 +{
  68 + // Initialise class members
  69 + init_members();
  70 +
  71 + // Perform foward transformation
  72 + forward(array);
  73 +
  74 + // Return
  75 + return;
  76 +}
  77 +
  78 +
  79 +/***********************************************************************//**
  80 + * @brief Copy constructor
  81 + *
  82 + * @param[in] fft Fast Fourier Transform.
  83 + ***************************************************************************/
  84 +GFft::GFft(const GFft& fft)
  85 +{
  86 + // Initialise class members
  87 + init_members();
  88 +
  89 + // Copy members
  90 + copy_members(fft);
  91 +
  92 + // Return
  93 + return;
  94 +}
  95 +
  96 +
  97 +/***********************************************************************//**
  98 + * @brief Destructor
  99 + ***************************************************************************/
  100 +GFft::~GFft(void)
  101 +{
  102 + // Free members
  103 + free_members();
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/*==========================================================================
  111 + = =
  112 + = Operators =
  113 + = =
  114 + ==========================================================================*/
  115 +
  116 +/***********************************************************************//**
  117 + * @brief Assignment operator
  118 + *
  119 + * @param[in] fft Fast Fourier Transform.
  120 + * @return Fast Fourier Transform.
  121 + ***************************************************************************/
  122 +GFft& GFft::operator=(const GFft& fft)
  123 +{
  124 + // Execute only if object is not identical
  125 + if (this != &fft) {
  126 +
  127 + // Free members
  128 + free_members();
  129 +
  130 + // Initialise private members
  131 + init_members();
  132 +
  133 + // Copy members
  134 + copy_members(fft);
  135 +
  136 + } // endif: object was not identical
  137 +
  138 + // Return this object
  139 + return *this;
  140 +}
  141 +
  142 +
  143 +/***********************************************************************//**
  144 + * @brief Unary addition operator
  145 + *
  146 + * @param[in] fft Fast Fourier Transform.
  147 + * @return Fast Fourier Transform.
  148 + *
  149 + * Adds a Fast Fourier Transform to another Fast Fourier Transform.
  150 + *
  151 + * @todo Implement operator
  152 + ***************************************************************************/
  153 +GFft& GFft::operator+=(const GFft& fft)
  154 +{
  155 + // TODO: Implement operator
  156 +
  157 + // Return Fast Fourier Transform
  158 + return *this;
  159 +}
  160 +
  161 +
  162 +/***********************************************************************//**
  163 + * @brief Unary subtraction operator
  164 + *
  165 + * @param[in] fft Fast Fourier transform.
  166 + * @return Fast Fourier transform.
  167 + *
  168 + * Subtracts a Fast Fourier Transform from another Fast Fourier Transform.
  169 + *
  170 + * @todo Implement operator
  171 + ***************************************************************************/
  172 +GFft& GFft::operator-=(const GFft& fft)
  173 +{
  174 + // TODO: Implement operator
  175 +
  176 + // Return Fast Fourier Transform
  177 + return *this;
  178 +}
  179 +
  180 +
  181 +/***********************************************************************//**
  182 + * @brief Unary multiplication operator
  183 + *
  184 + * @param[in] fft Fast Fourier transform.
  185 + * @return Fast Fourier transform.
  186 + *
  187 + * Multiplies a Fast Fourier Transform to another Fast Fourier Transform.
  188 + *
  189 + * @todo Implement operator
  190 + ***************************************************************************/
  191 +GFft& GFft::operator*=(const GFft& fft)
  192 +{
  193 + // TODO: Implement operator
  194 +
  195 + // Return Fast Fourier Transform
  196 + return *this;
  197 +}
  198 +
  199 +
  200 +/***********************************************************************//**
  201 + * @brief Unary division operator
  202 + *
  203 + * @param[in] fft Fast Fourier transform.
  204 + * @return Fast Fourier transform.
  205 + *
  206 + * Divides a Fast Fourier Transform by another Fast Fourier Transform.
  207 + *
  208 + * @todo Implement operator
  209 + ***************************************************************************/
  210 +GFft& GFft::operator/=(const GFft& fft)
  211 +{
  212 + // TODO: Implement operator
  213 +
  214 + // Return Fast Fourier Transform
  215 + return *this;
  216 +}
  217 +
  218 +
  219 +/***********************************************************************//**
  220 + * @brief Unary minus operator
  221 + *
  222 + * @return Fast Fourier transformation.
  223 + *
  224 + * Negates all elements of the Fast Fourier Transformation.
  225 + ***************************************************************************/
  226 +GFft GFft::operator-(void) const
  227 +{
  228 + // Copy FFT
  229 + GFft fft = *this;
  230 +
  231 + // Negate all elements
  232 + for (int i = 0; i < m_data.size(); ++i) {
  233 + fft.m_data[i] = -fft.m_data[i];
  234 + }
  235 +
  236 + // Return FFT
  237 + return fft;
  238 +}
  239 +
  240 +
  241 +/*==========================================================================
  242 + = =
  243 + = Public methods =
  244 + = =
  245 + ==========================================================================*/
  246 +
  247 +/***********************************************************************//**
  248 + * @brief Clear Fast Fourier Transform
  249 + ***************************************************************************/
  250 +void GFft::clear(void)
  251 +{
  252 + // Free members
  253 + free_members();
  254 +
  255 + // Initialise private members
  256 + init_members();
  257 +
  258 + // Return
  259 + return;
  260 +}
  261 +
  262 +
  263 +/***********************************************************************//**
  264 + * @brief Clone Fast Fourier Transform
  265 + *
  266 + * @return Pointer to deep copy of Fast Fourier Transform.
  267 + ***************************************************************************/
  268 +GFft* GFft::clone(void) const
  269 +{
  270 + // Clone FFT
  271 + return new GFft(*this);
  272 +}
  273 +
  274 +
  275 +/***********************************************************************//**
  276 + * @brief Forward Fast Fourier Transform
  277 + *
  278 + * @param[in] array N-dimensional array.
  279 + *
  280 + * @todo Method should support N-dim arrays, only supports 1-dim so far
  281 + ***************************************************************************/
  282 +void GFft::forward(const GNdarray& array)
  283 +{
  284 + // Set data
  285 + set_data(array);
  286 +
  287 + // Continue only if there are elements in the data member
  288 + if (m_data.size() > 0) {
  289 +
  290 + // Perform FFT (is only correct for 1D for now)
  291 + transform(m_data, 1, m_shape[0], m_wavetable[0], true);
  292 +
  293 + } // endif: there were elements in the FFT data member
  294 +
  295 + // Return
  296 + return;
  297 +}
  298 +
  299 +
  300 +/***********************************************************************//**
  301 + * @brief Backward Fast Fourier Transform
  302 + *
  303 + * @return N-dimensional array.
  304 + *
  305 + * @todo Method should support N-dim arrays, only supports 1-dim so far
  306 + ***************************************************************************/
  307 +GNdarray GFft::backward(void) const
  308 +{
  309 + // Initialise n-dimensional array
  310 + GNdarray array(shape());
  311 +
  312 + // Continue only if there are elements in the array
  313 + if (array.size() > 0) {
  314 +
  315 + // Create copy of complex array
  316 + std::vector<std::complex<double> > data = m_data;
  317 +
  318 + // Perform FFT (circumvent const correctness)
  319 + const_cast<GFft*>(this)->transform(data, 1, m_shape[0], m_wavetable[0],
  320 + false);
  321 +
  322 + // Get normalisation factor
  323 + const double norm = 1.0 / double(array.size());
  324 +
  325 + // Extract N-dimensional array and normalise with 1/n
  326 + for (int i = 0; i < array.size(); ++i) {
  327 + array(i) = data[i].real() * norm;
  328 + }
  329 +
  330 + } // endif: there were elements in the array
  331 +
  332 + // Return N-dimensional array
  333 + return array;
  334 +}
  335 +
  336 +
  337 +/***********************************************************************//**
  338 + * @brief Print Fast Fourier Transform information
  339 + *
  340 + * @param[in] chatter Chattiness.
  341 + * @return String containing Fast Fourier Transform information.
  342 + ***************************************************************************/
  343 +std::string GFft::print(const GChatter& chatter) const
  344 +{
  345 + // Initialise result string
  346 + std::string result;
  347 +
  348 + // Continue only if chatter is not silent
  349 + if (chatter != SILENT) {
  350 +
  351 + // Append header
  352 + result.append("=== GFft ===");
  353 +
  354 + // Append array information
  355 + result.append("\n"+gammalib::parformat("Dimension"));
  356 + result.append(gammalib::str(dim()));
  357 + result.append("\n"+gammalib::parformat("Shape"));
  358 + result.append("(");
  359 + for (int i = 0; i < dim(); ++i) {
  360 + if (i > 0) {
  361 + result.append(",");
  362 + }
  363 + result.append(gammalib::str(m_shape[i]));
  364 + }
  365 + result.append(")");
  366 + result.append("\n"+gammalib::parformat("Size"));
  367 + result.append(gammalib::str(size()));
  368 +
  369 + // VERBOSE: Put all FFT elements in string
  370 + if (chatter >= VERBOSE) {
  371 + result.append("\n"+gammalib::parformat("Elements"));
  372 + for (int i = 0; i < m_data.size(); ++i) {
  373 + if (i > 0) {
  374 + result.append(",");
  375 + }
  376 + result.append(gammalib::str(m_data[i]));
  377 + }
  378 + result.append(")");
  379 + }
  380 +
  381 + } // endif: chatter was not silent
  382 +
  383 + // Return result
  384 + return result;
  385 +}
  386 +
  387 +
  388 +/*==========================================================================
  389 + = =
  390 + = Private methods =
  391 + = =
  392 + ==========================================================================*/
  393 +
  394 +/***********************************************************************//**
  395 + * @brief Initialise class members
  396 + ***************************************************************************/
  397 +void GFft::init_members(void)
  398 +{
  399 + // Initialise members
  400 + m_shape.clear();
  401 + m_strides.clear();
  402 + m_data.clear();
  403 + m_wavetable.clear();
  404 +
  405 + // Return
  406 + return;
  407 +}
  408 +
  409 +
  410 +/***********************************************************************//**
  411 + * @brief Copy class members
  412 + *
  413 + * @param[in] fft Fast Fourier Transform.
  414 + ***************************************************************************/
  415 +void GFft::copy_members(const GFft& fft)
  416 +{
  417 + // Copy members
  418 + m_shape = fft.m_shape;
  419 + m_strides = fft.m_strides;
  420 + m_data = fft.m_data;
  421 + m_wavetable = fft.m_wavetable;
  422 +
  423 + // Return
  424 + return;
  425 +}
  426 +
  427 +
  428 +/***********************************************************************//**
  429 + * @brief Delete class members
  430 + ***************************************************************************/
  431 +void GFft::free_members(void)
  432 +{
  433 + // Return
  434 + return;
  435 +}
  436 +
  437 +
  438 +/***********************************************************************//**
  439 + * @brief Set data from n-dimensional array
  440 + *
  441 + * @param[in] array N-dimensional array
  442 + ***************************************************************************/
  443 +void GFft::set_data(const GNdarray& array)
  444 +{
  445 + // Clear information
  446 + m_shape.clear();
  447 + m_strides.clear();
  448 + m_data.clear();
  449 + m_wavetable.clear();
  450 +
  451 + // Continue only if there are elements in the array
  452 + if (array.size() > 0) {
  453 +
  454 + // Allocate data and initialise data to zero
  455 + std::complex<double> zero(0.0, 0.0);
  456 + m_data.assign(array.size(), zero);
  457 +
  458 + // Set real part of data from n-dimensional array
  459 + for (int i = 0; i < array.size(); ++i) {
  460 + m_data[i].real(array(i));
  461 + }
  462 +
  463 + // Save shape
  464 + m_shape = array.shape();
  465 +
  466 + // Set strides
  467 + int stride = 1;
  468 + for (int i = 0; i < m_shape.size(); ++i) {
  469 + m_strides.push_back(stride);
  470 + stride *= m_shape[i];
  471 + }
  472 +
  473 + // Set trigonometric coefficients
  474 + for (int i = 0; i < m_shape.size(); ++i) {
  475 + GFftWavetable wavetable(m_shape[i]);
  476 + m_wavetable.push_back(wavetable);
  477 + }
  478 +
  479 + } // endif: there were elements in the array
  480 +
  481 + // Return
  482 + return;
  483 +}
  484 +
  485 +
  486 +/***********************************************************************//**
  487 + * @brief Perform Fast Fourier Transform
  488 + *
  489 + * @param[in] data FFT array.
  490 + * @param[in] stride Step size when traversing FFT array in one dimension.
  491 + * @param[in] n Number of elements in the dimension.
  492 + * @param[in] wavetable Trigonometric coefficients.
  493 + * @param[in] forward Forward transform if true, otherwise backward transform.
  494 + ***************************************************************************/
  495 +void GFft::transform(std::vector<std::complex<double> >& data,
  496 + const int& stride,
  497 + const int& n,
  498 + const GFftWavetable& wavetable,
  499 + const bool& forward)
  500 +{
  501 + // Allocate scratch space
  502 + std::complex<double> zero(0.0, 0.0);
  503 + std::vector<std::complex<double> > scratch(n, zero);
  504 +
  505 + // Initialise transformation state
  506 + int state = 0;
  507 +
  508 + // Initialise factorisation product
  509 + int product = 1;
  510 +
  511 + // Set sign
  512 + int sign = (forward) ? -1 : +1;
  513 +
  514 + // Loop over all factors
  515 + for (int i = 0; i < wavetable.factors(); ++i) {
  516 +
  517 + // Get factorisation factor and start index
  518 + int factor = wavetable.factor(i);
  519 + int index = wavetable.index(i);
  520 +
  521 + // Update factorisation product
  522 + product *= factor;
  523 +
  524 + // Compute number of coefficients per factor
  525 + int q = n / product;
  526 +
  527 + // Set state dependent stuff
  528 + std::vector<std::complex<double> >* in = &data;
  529 + std::vector<std::complex<double> >* out = &scratch;
  530 + int istride;
  531 + int ostride;
  532 + if (state == 0) {
  533 + in = &data;
  534 + istride = stride;
  535 + out = &scratch;
  536 + ostride = 1;
  537 + state = 1;
  538 + }
  539 + else {
  540 + in = &scratch;
  541 + istride = 1;
  542 + out = &data;
  543 + ostride = stride;
  544 + state = 0;
  545 + }
  546 +
  547 + // Call factor dependent method
  548 + if (factor == 2) {
  549 + factor2(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  550 + }
  551 + else if (factor == 3) {
  552 + factor3(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  553 + }
  554 + else if (factor == 4) {
  555 + factor4(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  556 + }
  557 + else if (factor == 5) {
  558 + factor5(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  559 + }
  560 + else if (factor == 6) {
  561 + factor6(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  562 + }
  563 + else if (factor == 7) {
  564 + factor7(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  565 + }
  566 + else {
  567 + factorn(*in, istride, *out, ostride, wavetable, sign, factor, product, n, index);
  568 + }
  569 +
  570 + } // endfor: looped over all factors
  571 +
  572 + // In case the loop was exited in state 1 the results are in the scratch
  573 + // array and we need to copy the results back from the scratch array to
  574 + // the data array
  575 + if (state == 1) {
  576 + for (int i = 0; i < n; ++i) {
  577 + data[stride*i] = scratch[i];
  578 + }
  579 + }
  580 +
  581 + // Return
  582 + return;
  583 +}
  584 +
  585 +
  586 +/***********************************************************************//**
  587 + * @brief Compute FFT for factor 2
  588 + *
  589 + * @param[in] in Input array.
  590 + * @param[in] istride Step size when traversing input array.
  591 + * @param[in] out Output array.
  592 + * @param[in] ostride Step size when traversing output array.
  593 + * @param[in] wavetable Trigonometric coefficients.
  594 + * @param[in] sign Forward (-1) or backward (+1).
  595 + * @param[in] product ...
  596 + * @param[in] n Logical array length.
  597 + * @param[in] index Start index of trigonometric coefficients.
  598 + ***************************************************************************/
  599 +void GFft::factor2(const std::vector<std::complex<double> >& in,
  600 + const int& istride,
  601 + std::vector<std::complex<double> >& out,
  602 + const int& ostride,
  603 + const GFftWavetable& wavetable,
  604 + const int& sign,
  605 + const int& product,
  606 + const int& n,
  607 + const int& index)
  608 +{
  609 + // Compute ...
  610 + const int factor = 2;
  611 + const int m = n / factor;
  612 + const int q = n / product;
  613 + const int p_1 = product / factor;
  614 + const int jump = (factor - 1) * p_1;
  615 +
  616 + // Loop over ...
  617 + for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
  618 +
  619 + // Allocate variables
  620 + double w1_real, w1_imag;
  621 +
  622 + // Set trigonometric coefficients for k=0 since they are not stored
  623 + // in the wavetable object
  624 + if (k == 0) {
  625 + w1_real = 1.0;
  626 + w1_imag = 0.0;
  627 + }
  628 + else {
  629 +
  630 + // Compute indices
  631 + int twiddle = index + k - 1;
  632 +
  633 + // Set trigonometric coefficients for forward transform
  634 + if (sign == -1) {
  635 + w1_real = wavetable[twiddle].real();
  636 + w1_imag = wavetable[twiddle].imag();
  637 + }
  638 +
  639 + // ... otherwise set trigonometric coefficients for backward
  640 + // tranform: w -> conjugate(w)
  641 + else {
  642 + w1_real = wavetable[twiddle].real();
  643 + w1_imag = -wavetable[twiddle].imag();
  644 + }
  645 +
  646 + } // endelse: k != 0
  647 +
  648 + // Loop over ...
  649 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  650 +
  651 + const double z0_real = in[istride*i].real();
  652 + const double z0_imag = in[istride*i].imag();
  653 + const double z1_real = in[istride*(i+m)].real();
  654 + const double z1_imag = in[istride*(i+m)].imag();
  655 +
  656 + // compute x = W(2) z
  657 +
  658 + // x0 = z0 + z1
  659 + const double x0_real = z0_real + z1_real;
  660 + const double x0_imag = z0_imag + z1_imag;
  661 +
  662 + // x1 = z0 - z1
  663 + const double x1_real = z0_real - z1_real;
  664 + const double x1_imag = z0_imag - z1_imag;
  665 +
  666 + // out0 = 1 * x0
  667 + out[ostride*j].real(x0_real);
  668 + out[ostride*j].imag(x0_imag);
  669 +
  670 + // out1 = w1 * x1
  671 + out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
  672 + out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
  673 +
  674 + } // endfor: k1
  675 +
  676 + } // endfor: k
  677 +
  678 + // Return
  679 + return;
  680 +}
  681 +
  682 +
  683 +/***********************************************************************//**
  684 + * @brief Compute FFT for factor 3
  685 + *
  686 + * @param[in] in Input array.
  687 + * @param[in] istride Step size when traversing input array.
  688 + * @param[in] out Output array.
  689 + * @param[in] ostride Step size when traversing output array.
  690 + * @param[in] wavetable Trigonometric coefficients.
  691 + * @param[in] sign Forward (-1) or backward (+1).
  692 + * @param[in] product ...
  693 + * @param[in] n Logical array length.
  694 + * @param[in] index Start index of trigonometric coefficients.
  695 + ***************************************************************************/
  696 +void GFft::factor3(const std::vector<std::complex<double> >& in,
  697 + const int& istride,
  698 + std::vector<std::complex<double> >& out,
  699 + const int& ostride,
  700 + const GFftWavetable& wavetable,
  701 + const int& sign,
  702 + const int& product,
  703 + const int& n,
  704 + const int& index)
  705 +{
  706 + // Compute ...
  707 + const int factor = 3;
  708 + const int m = n / factor;
  709 + const int q = n / product;
  710 + const int p_1 = product / factor;
  711 + const int jump = (factor - 1) * p_1;
  712 +
  713 +
  714 + // Precompute some factors
  715 + const double tau = std::sqrt(3.0) / 2.0;
  716 +
  717 + // Loop over ...
  718 + for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
  719 +
  720 + // Allocate variables
  721 + double w1_real, w1_imag, w2_real, w2_imag;
  722 +
  723 + // Set trigonometric coefficients for k=0 since they are not stored
  724 + // in the wavetable object
  725 + if (k == 0) {
  726 + w1_real = 1.0;
  727 + w1_imag = 0.0;
  728 + w2_real = 1.0;
  729 + w2_imag = 0.0;
  730 + }
  731 + else {
  732 +
  733 + // Compute indices
  734 + int twiddle1 = index + k - 1;
  735 + int twiddle2 = twiddle1 + q;
  736 +
  737 + // Set trigonometric coefficients for forward transform
  738 + if (sign == -1) {
  739 + w1_real = wavetable[twiddle1].real();
  740 + w1_imag = wavetable[twiddle1].imag();
  741 + w2_real = wavetable[twiddle2].real();
  742 + w2_imag = wavetable[twiddle2].imag();
  743 + }
  744 +
  745 + // ... otherwise set trigonometric coefficients for backward
  746 + // tranform: w -> conjugate(w)
  747 + else {
  748 + w1_real = wavetable[twiddle1].real();
  749 + w1_imag = -wavetable[twiddle1].imag();
  750 + w2_real = wavetable[twiddle2].real();
  751 + w2_imag = -wavetable[twiddle2].imag();
  752 + }
  753 + } // endelse: k != 0
  754 +
  755 + // ...
  756 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  757 +
  758 + // Get z
  759 + const double z0_real = in[istride*i].real();
  760 + const double z0_imag = in[istride*i].imag();
  761 + const double z1_real = in[istride*(i+m)].real();
  762 + const double z1_imag = in[istride*(i+m)].imag();
  763 + const double z2_real = in[istride*(i+2*m)].real();
  764 + const double z2_imag = in[istride*(i+2*m)].imag();
  765 +
  766 + // compute x = W(3) z
  767 +
  768 + // t1 = z1 + z2
  769 + const double t1_real = z1_real + z2_real;
  770 + const double t1_imag = z1_imag + z2_imag;
  771 +
  772 + // t2 = z0 - t1/2
  773 + const double t2_real = z0_real - t1_real / 2.0;
  774 + const double t2_imag = z0_imag - t1_imag / 2.0;
  775 +
  776 + // t3 = (+/-) sin(pi/3)*(z1 - z2)
  777 + const double t3_real = ((int) sign) * tau * (z1_real - z2_real);
  778 + const double t3_imag = ((int) sign) * tau * (z1_imag - z2_imag);
  779 +
  780 + // x0 = z0 + t1
  781 + const double x0_real = z0_real + t1_real;
  782 + const double x0_imag = z0_imag + t1_imag;
  783 +
  784 + // x1 = t2 + i t3
  785 + const double x1_real = t2_real - t3_imag;
  786 + const double x1_imag = t2_imag + t3_real;
  787 +
  788 + // x2 = t2 - i t3
  789 + const double x2_real = t2_real + t3_imag;
  790 + const double x2_imag = t2_imag - t3_real;
  791 +
  792 + // out0 = 1 * x0
  793 + out[ostride*j].real(x0_real);
  794 + out[ostride*j].imag(x0_imag);
  795 +
  796 + // out1 = w1 * x1
  797 + out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
  798 + out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
  799 +
  800 + // out2 = w2 * x2
  801 + out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
  802 + out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
  803 +
  804 + } // endfor: k1
  805 +
  806 + } // endfor: k
  807 +
  808 + // Return
  809 + return;
  810 +}
  811 +
  812 +
  813 +/***********************************************************************//**
  814 + * @brief Compute FFT for factor 4
  815 + *
  816 + * @param[in] in Input array.
  817 + * @param[in] istride Step size when traversing input array.
  818 + * @param[in] out Output array.
  819 + * @param[in] ostride Step size when traversing output array.
  820 + * @param[in] wavetable Trigonometric coefficients.
  821 + * @param[in] sign Forward (-1) or backward (+1).
  822 + * @param[in] product ...
  823 + * @param[in] n Logical array length.
  824 + * @param[in] index Start index of trigonometric coefficients.
  825 + ***************************************************************************/
  826 +void GFft::factor4(const std::vector<std::complex<double> >& in,
  827 + const int& istride,
  828 + std::vector<std::complex<double> >& out,
  829 + const int& ostride,
  830 + const GFftWavetable& wavetable,
  831 + const int& sign,
  832 + const int& product,
  833 + const int& n,
  834 + const int& index)
  835 +{
  836 + // Compute ...
  837 + const int factor = 4;
  838 + const int m = n / factor;
  839 + const int q = n / product;
  840 + const int p_1 = product / factor;
  841 + const int jump = (factor - 1) * p_1;
  842 +
  843 + // Loop over ...
  844 + for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
  845 +
  846 + // Allocate w
  847 + double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag;
  848 +
  849 + // Set trigonometric coefficients for k=0 since they are not stored
  850 + // in the wavetable object
  851 + if (k == 0) {
  852 + w1_real = 1.0;
  853 + w1_imag = 0.0;
  854 + w2_real = 1.0;
  855 + w2_imag = 0.0;
  856 + w3_real = 1.0;
  857 + w3_imag = 0.0;
  858 + }
  859 + else
  860 + {
  861 +
  862 + // Compute indices
  863 + int twiddle1 = index + k - 1;
  864 + int twiddle2 = twiddle1 + q;
  865 + int twiddle3 = twiddle2 + q;
  866 +
  867 + // Set trigonometric coefficients for forward transform
  868 + if (sign == -1) {
  869 + w1_real = wavetable[twiddle1].real();
  870 + w1_imag = wavetable[twiddle1].imag();
  871 + w2_real = wavetable[twiddle2].real();
  872 + w2_imag = wavetable[twiddle2].imag();
  873 + w3_real = wavetable[twiddle3].real();
  874 + w3_imag = wavetable[twiddle3].imag();
  875 + }
  876 +
  877 + // ... otherwise set trigonometric coefficients for backward
  878 + // tranform: w -> conjugate(w)
  879 + else {
  880 + w1_real = wavetable[twiddle1].real();
  881 + w1_imag = -wavetable[twiddle1].imag();
  882 + w2_real = wavetable[twiddle2].real();
  883 + w2_imag = -wavetable[twiddle2].imag();
  884 + w3_real = wavetable[twiddle3].real();
  885 + w3_imag = -wavetable[twiddle3].imag();
  886 + }
  887 +
  888 + } // endelse: k != 0
  889 +
  890 + // Loop over ...
  891 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  892 +
  893 +
  894 + // Get z
  895 + const double z0_real = in[istride*i].real();
  896 + const double z0_imag = in[istride*i].imag();
  897 + const double z1_real = in[istride*(i+m)].real();
  898 + const double z1_imag = in[istride*(i+m)].imag();
  899 + const double z2_real = in[istride*(i+2*m)].real();
  900 + const double z2_imag = in[istride*(i+2*m)].imag();
  901 + const double z3_real = in[istride*(i+3*m)].real();
  902 + const double z3_imag = in[istride*(i+3*m)].imag();
  903 +
  904 + // compute x = W(4) z
  905 +
  906 + // t1 = z0 + z2
  907 + const double t1_real = z0_real + z2_real;
  908 + const double t1_imag = z0_imag + z2_imag;
  909 +
  910 + // t2 = z1 + z3
  911 + const double t2_real = z1_real + z3_real;
  912 + const double t2_imag = z1_imag + z3_imag;
  913 +
  914 + // t3 = z0 - z2
  915 + const double t3_real = z0_real - z2_real;
  916 + const double t3_imag = z0_imag - z2_imag;
  917 +
  918 + // t4 = (+/-) (z1 - z3)
  919 + const double t4_real = ((int) sign) * (z1_real - z3_real);
  920 + const double t4_imag = ((int) sign) * (z1_imag - z3_imag);
  921 +
  922 + // x0 = t1 + t2
  923 + const double x0_real = t1_real + t2_real;
  924 + const double x0_imag = t1_imag + t2_imag;
  925 +
  926 + // x1 = t3 + i t4
  927 + const double x1_real = t3_real - t4_imag;
  928 + const double x1_imag = t3_imag + t4_real;
  929 +
  930 + // x2 = t1 - t2
  931 + const double x2_real = t1_real - t2_real;
  932 + const double x2_imag = t1_imag - t2_imag;
  933 +
  934 + // x3 = t3 - i t4
  935 + const double x3_real = t3_real + t4_imag;
  936 + const double x3_imag = t3_imag - t4_real;
  937 +
  938 + // out0 = 1 * x0
  939 + out[ostride*j].real(x0_real);
  940 + out[ostride*j].imag(x0_imag);
  941 +
  942 + // out1 = w1 * x1
  943 + out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
  944 + out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
  945 +
  946 + // out2 = w2 * x2
  947 + out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
  948 + out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
  949 +
  950 + // out3 = w3 * x3
  951 + out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
  952 + out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
  953 +
  954 + } // endfor: k1
  955 +
  956 + } // endfor: k
  957 +
  958 + // Return
  959 + return;
  960 +}
  961 +
  962 +
  963 +/***********************************************************************//**
  964 + * @brief Compute FFT for factor 5
  965 + *
  966 + * @param[in] in Input array.
  967 + * @param[in] istride Step size when traversing input array.
  968 + * @param[in] out Output array.
  969 + * @param[in] ostride Step size when traversing output array.
  970 + * @param[in] wavetable Trigonometric coefficients.
  971 + * @param[in] sign Forward (-1) or backward (+1).
  972 + * @param[in] product ...
  973 + * @param[in] n Logical array length.
  974 + * @param[in] index Start index of trigonometric coefficients.
  975 + ***************************************************************************/
  976 +void GFft::factor5(const std::vector<std::complex<double> >& in,
  977 + const int& istride,
  978 + std::vector<std::complex<double> >& out,
  979 + const int& ostride,
  980 + const GFftWavetable& wavetable,
  981 + const int& sign,
  982 + const int& product,
  983 + const int& n,
  984 + const int& index)
  985 +{
  986 + // Compute ...
  987 + const int factor = 5;
  988 + const int m = n / factor;
  989 + const int q = n / product;
  990 + const int p_1 = product / factor;
  991 + const int jump = (factor - 1) * p_1;
  992 +
  993 + // Precompute some factors
  994 + const double sin_2pi_by_5 = std::sin(gammalib::twopi / 5.0);
  995 + const double sin_2pi_by_10 = std::sin(gammalib::twopi / 10.0);
  996 + const double tau = std::sqrt(5.0) / 4.0;
  997 +
  998 + // Loop over ...
  999 + for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
  1000 +
  1001 + // Allocate w
  1002 + double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag, w4_real,
  1003 + w4_imag;
  1004 +
  1005 + // Set trigonometric coefficients for k=0 since they are not stored
  1006 + // in the wavetable object
  1007 + if (k == 0) {
  1008 + w1_real = 1.0;
  1009 + w1_imag = 0.0;
  1010 + w2_real = 1.0;
  1011 + w2_imag = 0.0;
  1012 + w3_real = 1.0;
  1013 + w3_imag = 0.0;
  1014 + w4_real = 1.0;
  1015 + w4_imag = 0.0;
  1016 + }
  1017 + else {
  1018 +
  1019 + // Compute indices
  1020 + int twiddle1 = index + k - 1;
  1021 + int twiddle2 = twiddle1 + q;
  1022 + int twiddle3 = twiddle2 + q;
  1023 + int twiddle4 = twiddle3 + q;
  1024 +
  1025 + // Set trigonometric coefficients for forward transform
  1026 + if (sign == -1) {
  1027 + w1_real = wavetable[twiddle1].real();
  1028 + w1_imag = wavetable[twiddle1].imag();
  1029 + w2_real = wavetable[twiddle2].real();
  1030 + w2_imag = wavetable[twiddle2].imag();
  1031 + w3_real = wavetable[twiddle3].real();
  1032 + w3_imag = wavetable[twiddle3].imag();
  1033 + w4_real = wavetable[twiddle4].real();
  1034 + w4_imag = wavetable[twiddle4].imag();
  1035 + }
  1036 +
  1037 + // ... otherwise set trigonometric coefficients for backward
  1038 + // tranform: w -> conjugate(w)
  1039 + else {
  1040 + w1_real = wavetable[twiddle1].real();
  1041 + w1_imag = -wavetable[twiddle1].imag();
  1042 + w2_real = wavetable[twiddle2].real();
  1043 + w2_imag = -wavetable[twiddle2].imag();
  1044 + w3_real = wavetable[twiddle3].real();
  1045 + w3_imag = -wavetable[twiddle3].imag();
  1046 + w4_real = wavetable[twiddle4].real();
  1047 + w4_imag = -wavetable[twiddle4].imag();
  1048 + }
  1049 +
  1050 + } // endelse: k != 0
  1051 +
  1052 + // Loop over ...
  1053 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  1054 +
  1055 + // Get z
  1056 + const double z0_real = in[istride*i].real();
  1057 + const double z0_imag = in[istride*i].imag();
  1058 + const double z1_real = in[istride*(i+m)].real();
  1059 + const double z1_imag = in[istride*(i+m)].imag();
  1060 + const double z2_real = in[istride*(i+2*m)].real();
  1061 + const double z2_imag = in[istride*(i+2*m)].imag();
  1062 + const double z3_real = in[istride*(i+3*m)].real();
  1063 + const double z3_imag = in[istride*(i+3*m)].imag();
  1064 + const double z4_real = in[istride*(i+4*m)].real();
  1065 + const double z4_imag = in[istride*(i+4*m)].imag();
  1066 +
  1067 + // compute x = W(5) z
  1068 +
  1069 + // t1 = z1 + z4
  1070 + const double t1_real = z1_real + z4_real;
  1071 + const double t1_imag = z1_imag + z4_imag;
  1072 +
  1073 + // t2 = z2 + z3
  1074 + const double t2_real = z2_real + z3_real;
  1075 + const double t2_imag = z2_imag + z3_imag;
  1076 +
  1077 + // t3 = z1 - z4
  1078 + const double t3_real = z1_real - z4_real;
  1079 + const double t3_imag = z1_imag - z4_imag;
  1080 +
  1081 + // t4 = z2 - z3
  1082 + const double t4_real = z2_real - z3_real;
  1083 + const double t4_imag = z2_imag - z3_imag;
  1084 +
  1085 + // t5 = t1 + t2
  1086 + const double t5_real = t1_real + t2_real;
  1087 + const double t5_imag = t1_imag + t2_imag;
  1088 +
  1089 + // t6 = (sqrt(5)/4)(t1 - t2)
  1090 + const double t6_real = (tau) * (t1_real - t2_real);
  1091 + const double t6_imag = (tau) * (t1_imag - t2_imag);
  1092 +
  1093 + // t7 = z0 - ((t5)/4)
  1094 + const double t7_real = z0_real - t5_real / 4.0;
  1095 + const double t7_imag = z0_imag - t5_imag / 4.0;
  1096 +
  1097 + // t8 = t7 + t6
  1098 + const double t8_real = t7_real + t6_real;
  1099 + const double t8_imag = t7_imag + t6_imag;
  1100 +
  1101 + // t9 = t7 - t6
  1102 + const double t9_real = t7_real - t6_real;
  1103 + const double t9_imag = t7_imag - t6_imag;
  1104 +
  1105 + // t10 = sin(2 pi/5) t3 + sin(2 pi/10) t4 */
  1106 + const double t10_real = ((int) sign) * (sin_2pi_by_5 * t3_real +
  1107 + sin_2pi_by_10 * t4_real);
  1108 + const double t10_imag = ((int) sign) * (sin_2pi_by_5 * t3_imag +
  1109 + sin_2pi_by_10 * t4_imag);
  1110 +
  1111 + // t11 = sin(2 pi/10) t3 - sin(2 pi/5) t4
  1112 + const double t11_real = ((int) sign) * (sin_2pi_by_10 * t3_real -
  1113 + sin_2pi_by_5 * t4_real);
  1114 + const double t11_imag = ((int) sign) * (sin_2pi_by_10 * t3_imag -
  1115 + sin_2pi_by_5 * t4_imag);
  1116 +
  1117 + // x0 = z0 + t5
  1118 + const double x0_real = z0_real + t5_real;
  1119 + const double x0_imag = z0_imag + t5_imag;
  1120 +
  1121 + // x1 = t8 + i t10
  1122 + const double x1_real = t8_real - t10_imag;
  1123 + const double x1_imag = t8_imag + t10_real;
  1124 +
  1125 + // x2 = t9 + i t11
  1126 + const double x2_real = t9_real - t11_imag;
  1127 + const double x2_imag = t9_imag + t11_real;
  1128 +
  1129 + // x3 = t9 - i t11
  1130 + const double x3_real = t9_real + t11_imag;
  1131 + const double x3_imag = t9_imag - t11_real;
  1132 +
  1133 + // x4 = t8 - i t10
  1134 + const double x4_real = t8_real + t10_imag;
  1135 + const double x4_imag = t8_imag - t10_real;
  1136 +
  1137 + // out0 = 1 * x0
  1138 + out[ostride*j].real(x0_real);
  1139 + out[ostride*j].imag(x0_imag);
  1140 +
  1141 + // out1 = w1 * x1
  1142 + out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
  1143 + out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
  1144 +
  1145 + // out2 = w2 * x2
  1146 + out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
  1147 + out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
  1148 +
  1149 + // out3 = w3 * x3
  1150 + out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
  1151 + out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
  1152 +
  1153 + // out4 = w4 * x4
  1154 + out[ostride*(j+4*p_1)].real(w4_real * x4_real - w4_imag * x4_imag);
  1155 + out[ostride*(j+4*p_1)].imag(w4_real * x4_imag + w4_imag * x4_real);
  1156 +
  1157 + } // endfor: k1
  1158 +
  1159 + } // endfor: k
  1160 +
  1161 + // Return
  1162 + return;
  1163 +}
  1164 +
  1165 +
  1166 +/***********************************************************************//**
  1167 + * @brief Compute FFT for factor 6
  1168 + *
  1169 + * @param[in] in Input array.
  1170 + * @param[in] istride Step size when traversing input array.
  1171 + * @param[in] out Output array.
  1172 + * @param[in] ostride Step size when traversing output array.
  1173 + * @param[in] wavetable Trigonometric coefficients.
  1174 + * @param[in] sign Forward (-1) or backward (+1).
  1175 + * @param[in] product ...
  1176 + * @param[in] n Logical array length.
  1177 + * @param[in] index Start index of trigonometric coefficients.
  1178 + ***************************************************************************/
  1179 +void GFft::factor6(const std::vector<std::complex<double> >& in,
  1180 + const int& istride,
  1181 + std::vector<std::complex<double> >& out,
  1182 + const int& ostride,
  1183 + const GFftWavetable& wavetable,
  1184 + const int& sign,
  1185 + const int& product,
  1186 + const int& n,
  1187 + const int& index)
  1188 +{
  1189 + // Compute ...
  1190 + const int factor = 6;
  1191 + const int m = n / factor;
  1192 + const int q = n / product;
  1193 + const int p_1 = product / factor;
  1194 + const int jump = (factor - 1) * p_1;
  1195 +
  1196 + // Precompute some factors
  1197 + const double tau = std::sqrt(3.0) / 2.0;
  1198 +
  1199 + // Loop over ...
  1200 + for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
  1201 +
  1202 + // Allocate w
  1203 + double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag, w4_real,
  1204 + w4_imag, w5_real, w5_imag;
  1205 +
  1206 + // Set trigonometric coefficients for k=0 since they are not stored
  1207 + // in the wavetable object
  1208 + if (k == 0) {
  1209 + w1_real = 1.0;
  1210 + w1_imag = 0.0;
  1211 + w2_real = 1.0;
  1212 + w2_imag = 0.0;
  1213 + w3_real = 1.0;
  1214 + w3_imag = 0.0;
  1215 + w4_real = 1.0;
  1216 + w4_imag = 0.0;
  1217 + w5_real = 1.0;
  1218 + w5_imag = 0.0;
  1219 + }
  1220 + else {
  1221 +
  1222 + // Compute indices
  1223 + int twiddle1 = index + k - 1;
  1224 + int twiddle2 = twiddle1 + q;
  1225 + int twiddle3 = twiddle2 + q;
  1226 + int twiddle4 = twiddle3 + q;
  1227 + int twiddle5 = twiddle4 + q;
  1228 +
  1229 + // Set trigonometric coefficients for forward transform
  1230 + if (sign == -1) {
  1231 + w1_real = wavetable[twiddle1].real();
  1232 + w1_imag = wavetable[twiddle1].imag();
  1233 + w2_real = wavetable[twiddle2].real();
  1234 + w2_imag = wavetable[twiddle2].imag();
  1235 + w3_real = wavetable[twiddle3].real();
  1236 + w3_imag = wavetable[twiddle3].imag();
  1237 + w4_real = wavetable[twiddle4].real();
  1238 + w4_imag = wavetable[twiddle4].imag();
  1239 + w5_real = wavetable[twiddle5].real();
  1240 + w5_imag = wavetable[twiddle5].imag();
  1241 + }
  1242 +
  1243 + // ... otherwise set trigonometric coefficients for backward
  1244 + // tranform: w -> conjugate(w)
  1245 + else {
  1246 + w1_real = wavetable[twiddle1].real();
  1247 + w1_imag = -wavetable[twiddle1].imag();
  1248 + w2_real = wavetable[twiddle2].real();
  1249 + w2_imag = -wavetable[twiddle2].imag();
  1250 + w3_real = wavetable[twiddle3].real();
  1251 + w3_imag = -wavetable[twiddle3].imag();
  1252 + w4_real = wavetable[twiddle4].real();
  1253 + w4_imag = -wavetable[twiddle4].imag();
  1254 + w5_real = wavetable[twiddle5].real();
  1255 + w5_imag = -wavetable[twiddle5].imag();
  1256 + }
  1257 +
  1258 + } // endelse: k != 0
  1259 +
  1260 + // ...
  1261 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  1262 +
  1263 + // Get z
  1264 + const double z0_real = in[istride*i].real();
  1265 + const double z0_imag = in[istride*i].imag();
  1266 + const double z1_real = in[istride*(i+m)].real();
  1267 + const double z1_imag = in[istride*(i+m)].imag();
  1268 + const double z2_real = in[istride*(i+2*m)].real();
  1269 + const double z2_imag = in[istride*(i+2*m)].imag();
  1270 + const double z3_real = in[istride*(i+3*m)].real();
  1271 + const double z3_imag = in[istride*(i+3*m)].imag();
  1272 + const double z4_real = in[istride*(i+4*m)].real();
  1273 + const double z4_imag = in[istride*(i+4*m)].imag();
  1274 + const double z5_real = in[istride*(i+5*m)].real();
  1275 + const double z5_imag = in[istride*(i+5*m)].imag();
  1276 +
  1277 + // compute x = W(6) z
  1278 +
  1279 + // W(6) is a combination of sums and differences of W(3) acting
  1280 + // on the even and odd elements of z
  1281 +
  1282 + // ta1 = z2 + z4
  1283 + const double ta1_real = z2_real + z4_real;
  1284 + const double ta1_imag = z2_imag + z4_imag;
  1285 +
  1286 + // ta2 = z0 - ta1/2
  1287 + const double ta2_real = z0_real - ta1_real / 2.0;
  1288 + const double ta2_imag = z0_imag - ta1_imag / 2.0;
  1289 +
  1290 + // ta3 = (+/-) sin(pi/3)*(z2 - z4)
  1291 + const double ta3_real = ((int) sign) * tau * (z2_real - z4_real);
  1292 + const double ta3_imag = ((int) sign) * tau * (z2_imag - z4_imag);
  1293 +
  1294 + // a0 = z0 + ta1
  1295 + const double a0_real = z0_real + ta1_real;
  1296 + const double a0_imag = z0_imag + ta1_imag;
  1297 +
  1298 + // a1 = ta2 + i ta3
  1299 + const double a1_real = ta2_real - ta3_imag;
  1300 + const double a1_imag = ta2_imag + ta3_real;
  1301 +
  1302 + // a2 = ta2 - i ta3
  1303 + const double a2_real = ta2_real + ta3_imag;
  1304 + const double a2_imag = ta2_imag - ta3_real;
  1305 +
  1306 + // tb1 = z5 + z1
  1307 + const double tb1_real = z5_real + z1_real;
  1308 + const double tb1_imag = z5_imag + z1_imag;
  1309 +
  1310 + // tb2 = z3 - tb1/2
  1311 + const double tb2_real = z3_real - tb1_real / 2.0;
  1312 + const double tb2_imag = z3_imag - tb1_imag / 2.0;
  1313 +
  1314 + // tb3 = (+/-) sin(pi/3)*(z5 - z1)
  1315 + const double tb3_real = ((int) sign) * tau * (z5_real - z1_real);
  1316 + const double tb3_imag = ((int) sign) * tau * (z5_imag - z1_imag);
  1317 +
  1318 + // b0 = z3 + tb1
  1319 + const double b0_real = z3_real + tb1_real;
  1320 + const double b0_imag = z3_imag + tb1_imag;
  1321 +
  1322 + // b1 = tb2 + i tb3
  1323 + const double b1_real = tb2_real - tb3_imag;
  1324 + const double b1_imag = tb2_imag + tb3_real;
  1325 +
  1326 + // b2 = tb2 - i tb3
  1327 + const double b2_real = tb2_real + tb3_imag;
  1328 + const double b2_imag = tb2_imag - tb3_real;
  1329 +
  1330 + // x0 = a0 + b0
  1331 + const double x0_real = a0_real + b0_real;
  1332 + const double x0_imag = a0_imag + b0_imag;
  1333 +
  1334 + // x4 = a1 + b1
  1335 + const double x4_real = a1_real + b1_real;
  1336 + const double x4_imag = a1_imag + b1_imag;
  1337 +
  1338 + // x2 = a2 + b2
  1339 + const double x2_real = a2_real + b2_real;
  1340 + const double x2_imag = a2_imag + b2_imag;
  1341 +
  1342 + // x3 = a0 - b0
  1343 + const double x3_real = a0_real - b0_real;
  1344 + const double x3_imag = a0_imag - b0_imag;
  1345 +
  1346 + // x1 = a1 - b1
  1347 + const double x1_real = a1_real - b1_real;
  1348 + const double x1_imag = a1_imag - b1_imag;
  1349 +
  1350 + // x5 = a2 - b2
  1351 + const double x5_real = a2_real - b2_real;
  1352 + const double x5_imag = a2_imag - b2_imag;
  1353 +
  1354 + // out0 = 1 * x0
  1355 + out[ostride*j].real(x0_real);
  1356 + out[ostride*j].imag(x0_imag);
  1357 +
  1358 + // out1 = w1 * x1
  1359 + out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
  1360 + out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
  1361 +
  1362 + // out2 = w2 * x2
  1363 + out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
  1364 + out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
  1365 +
  1366 + // out3 = w3 * x3
  1367 + out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
  1368 + out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
  1369 +
  1370 + // out4 = w4 * x4
  1371 + out[ostride*(j+4*p_1)].real(w4_real * x4_real - w4_imag * x4_imag);
  1372 + out[ostride*(j+4*p_1)].imag(w4_real * x4_imag + w4_imag * x4_real);
  1373 +
  1374 + // out5 = w5 * x5
  1375 + out[ostride*(j+5*p_1)].real(w5_real * x5_real - w5_imag * x5_imag);
  1376 + out[ostride*(j+5*p_1)].imag(w5_real * x5_imag + w5_imag * x5_real);
  1377 +
  1378 + } // endfor: k1
  1379 +
  1380 + } // endfor: k
  1381 +
  1382 + // Return
  1383 + return;
  1384 +}
  1385 +
  1386 +
  1387 +/***********************************************************************//**
  1388 + * @brief Compute FFT for factor 7
  1389 + *
  1390 + * @param[in] in Input array.
  1391 + * @param[in] istride Step size when traversing input array.
  1392 + * @param[in] out Output array.
  1393 + * @param[in] ostride Step size when traversing output array.
  1394 + * @param[in] wavetable Trigonometric coefficients.
  1395 + * @param[in] sign Forward (-1) or backward (+1).
  1396 + * @param[in] product ...
  1397 + * @param[in] n Logical array length.
  1398 + * @param[in] index Start index of trigonometric coefficients.
  1399 + ***************************************************************************/
  1400 +void GFft::factor7(const std::vector<std::complex<double> >& in,
  1401 + const int& istride,
  1402 + std::vector<std::complex<double> >& out,
  1403 + const int& ostride,
  1404 + const GFftWavetable& wavetable,
  1405 + const int& sign,
  1406 + const int& product,
  1407 + const int& n,
  1408 + const int& index)
  1409 +{
  1410 + // Compute ...
  1411 + const int factor = 7;
  1412 + const int m = n / factor;
  1413 + const int q = n / product;
  1414 + const int p_1 = product / factor;
  1415 + const int jump = (factor - 1) * p_1;
  1416 +
  1417 + // Precompute some factors
  1418 + static const double twopi7 = gammalib::twopi / 7.0;
  1419 + static const double c1 = std::cos(1.0 * twopi7);
  1420 + static const double c2 = std::cos(2.0 * twopi7);
  1421 + static const double c3 = std::cos(3.0 * twopi7);
  1422 + static const double s1 = std::sin(1.0 * twopi7);
  1423 + static const double s2 = std::sin(2.0 * twopi7);
  1424 + static const double s3 = std::sin(3.0 * twopi7);
  1425 + static const double tau1 = (c1 + c2 + c3) / 3.0 - 1.0;
  1426 + static const double tau2 = (2.0 * c1 - c2 - c3) / 3.0;
  1427 + static const double tau3 = (c1 - 2.0*c2 + c3) / 3.0;
  1428 + static const double tau4 = (c1 + c2 - 2.0 * c3) / 3.0;
  1429 + static const double tau5 = (s1 + s2 - s3) / 3.0;
  1430 + static const double tau6 = (2.0 * s1 - s2 + s3) / 3.0;
  1431 + static const double tau7 = (s1 - 2.0 * s2 - s3) / 3.0;
  1432 + static const double tau8 = (s1 + s2 + 2.0 * s3) / 3.0;
  1433 +
  1434 + // Loop over ...
  1435 + for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
  1436 +
  1437 + //
  1438 + double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag, w4_real,
  1439 + w4_imag, w5_real, w5_imag, w6_real, w6_imag;
  1440 +
  1441 + // Set trigonometric coefficients for k=0 since they are not stored
  1442 + // in the wavetable object
  1443 + if (k == 0) {
  1444 + w1_real = 1.0;
  1445 + w1_imag = 0.0;
  1446 + w2_real = 1.0;
  1447 + w2_imag = 0.0;
  1448 + w3_real = 1.0;
  1449 + w3_imag = 0.0;
  1450 + w4_real = 1.0;
  1451 + w4_imag = 0.0;
  1452 + w5_real = 1.0;
  1453 + w5_imag = 0.0;
  1454 + w6_real = 1.0;
  1455 + w6_imag = 0.0;
  1456 + }
  1457 + else {
  1458 +
  1459 + // Compute indices
  1460 + int twiddle1 = index + k - 1;
  1461 + int twiddle2 = twiddle1 + q;
  1462 + int twiddle3 = twiddle2 + q;
  1463 + int twiddle4 = twiddle3 + q;
  1464 + int twiddle5 = twiddle4 + q;
  1465 + int twiddle6 = twiddle5 + q;
  1466 +
  1467 + // Set trigonometric coefficients for forward transform
  1468 + if (sign == -1) {
  1469 + w1_real = wavetable[twiddle1].real();
  1470 + w1_imag = wavetable[twiddle1].imag();
  1471 + w2_real = wavetable[twiddle2].real();
  1472 + w2_imag = wavetable[twiddle2].imag();
  1473 + w3_real = wavetable[twiddle3].real();
  1474 + w3_imag = wavetable[twiddle3].imag();
  1475 + w4_real = wavetable[twiddle4].real();
  1476 + w4_imag = wavetable[twiddle4].imag();
  1477 + w5_real = wavetable[twiddle5].real();
  1478 + w5_imag = wavetable[twiddle5].imag();
  1479 + w6_real = wavetable[twiddle6].real();
  1480 + w6_imag = wavetable[twiddle6].imag();
  1481 + }
  1482 +
  1483 + // ... otherwise set trigonometric coefficients for backward
  1484 + // tranform: w -> conjugate(w)
  1485 + else {
  1486 + w1_real = wavetable[twiddle1].real();
  1487 + w1_imag = -wavetable[twiddle1].imag();
  1488 + w2_real = wavetable[twiddle2].real();
  1489 + w2_imag = -wavetable[twiddle2].imag();
  1490 + w3_real = wavetable[twiddle3].real();
  1491 + w3_imag = -wavetable[twiddle3].imag();
  1492 + w4_real = wavetable[twiddle4].real();
  1493 + w4_imag = -wavetable[twiddle4].imag();
  1494 + w5_real = wavetable[twiddle5].real();
  1495 + w5_imag = -wavetable[twiddle5].imag();
  1496 + w6_real = wavetable[twiddle6].real();
  1497 + w6_imag = -wavetable[twiddle6].imag();
  1498 + }
  1499 +
  1500 + } // endelse: k != 0
  1501 +
  1502 + // ...
  1503 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  1504 +
  1505 + // Get z
  1506 + const double z0_real = in[istride*i].real();
  1507 + const double z0_imag = in[istride*i].imag();
  1508 + const double z1_real = in[istride*(i+m)].real();
  1509 + const double z1_imag = in[istride*(i+m)].imag();
  1510 + const double z2_real = in[istride*(i+2*m)].real();
  1511 + const double z2_imag = in[istride*(i+2*m)].imag();
  1512 + const double z3_real = in[istride*(i+3*m)].real();
  1513 + const double z3_imag = in[istride*(i+3*m)].imag();
  1514 + const double z4_real = in[istride*(i+4*m)].real();
  1515 + const double z4_imag = in[istride*(i+4*m)].imag();
  1516 + const double z5_real = in[istride*(i+5*m)].real();
  1517 + const double z5_imag = in[istride*(i+5*m)].imag();
  1518 + const double z6_real = in[istride*(i+6*m)].real();
  1519 + const double z6_imag = in[istride*(i+6*m)].imag();
  1520 +
  1521 + // Compute x = W(7) z
  1522 +
  1523 + // t0 = z1 + z6
  1524 + const double t0_real = z1_real + z6_real;
  1525 + const double t0_imag = z1_imag + z6_imag;
  1526 +
  1527 + // t1 = z1 - z6
  1528 + const double t1_real = z1_real - z6_real;
  1529 + const double t1_imag = z1_imag - z6_imag;
  1530 +
  1531 + // t2 = z2 + z5
  1532 + const double t2_real = z2_real + z5_real;
  1533 + const double t2_imag = z2_imag + z5_imag;
  1534 +
  1535 + // t3 = z2 - z5
  1536 + const double t3_real = z2_real - z5_real;
  1537 + const double t3_imag = z2_imag - z5_imag;
  1538 +
  1539 + // t4 = z4 + z3
  1540 + const double t4_real = z4_real + z3_real;
  1541 + const double t4_imag = z4_imag + z3_imag;
  1542 +
  1543 + // t5 = z4 - z3
  1544 + const double t5_real = z4_real - z3_real;
  1545 + const double t5_imag = z4_imag - z3_imag;
  1546 +
  1547 + // t6 = t2 + t0
  1548 + const double t6_real = t2_real + t0_real;
  1549 + const double t6_imag = t2_imag + t0_imag;
  1550 +
  1551 + // t7 = t5 + t3
  1552 + const double t7_real = t5_real + t3_real;
  1553 + const double t7_imag = t5_imag + t3_imag;
  1554 +
  1555 + // b0 = z0 + t6 + t4
  1556 + const double b0_real = z0_real + t6_real + t4_real;
  1557 + const double b0_imag = z0_imag + t6_imag + t4_imag;
  1558 +
  1559 + // b1 = ((cos(2pi/7) + cos(4pi/7) + cos(6pi/7))/3-1) (t6 + t4)
  1560 + const double b1_real = tau1 * (t6_real + t4_real);
  1561 + const double b1_imag = tau1 * (t6_imag + t4_imag);
  1562 +
  1563 + // b2 = ((2*cos(2pi/7) - cos(4pi/7) - cos(6pi/7))/3) (t0 - t4)
  1564 + const double b2_real = tau2 * (t0_real - t4_real);
  1565 + const double b2_imag = tau2 * (t0_imag - t4_imag);
  1566 +
  1567 + // b3 = ((cos(2pi/7) - 2*cos(4pi/7) + cos(6pi/7))/3) (t4 - t2)
  1568 + const double b3_real = tau3 * (t4_real - t2_real);
  1569 + const double b3_imag = tau3 * (t4_imag - t2_imag);
  1570 +
  1571 + // b4 = ((cos(2pi/7) + cos(4pi/7) - 2*cos(6pi/7))/3) (t2 - t0)
  1572 + const double b4_real = tau4 * (t2_real - t0_real);
  1573 + const double b4_imag = tau4 * (t2_imag - t0_imag);
  1574 +
  1575 + // b5 = sign * ((sin(2pi/7) + sin(4pi/7) - sin(6pi/7))/3) (t7 + t1)
  1576 + const double b5_real = (-(int)sign) * tau5 * (t7_real + t1_real);
  1577 + const double b5_imag = (-(int)sign) * tau5 * (t7_imag + t1_imag);
  1578 +
  1579 + // b6 = sign * ((2sin(2pi/7) - sin(4pi/7) + sin(6pi/7))/3) (t1 - t5)
  1580 + const double b6_real = (-(int)sign) * tau6 * (t1_real - t5_real);
  1581 + const double b6_imag = (-(int)sign) * tau6 * (t1_imag - t5_imag);
  1582 +
  1583 + // b7 = sign * ((sin(2pi/7) - 2sin(4pi/7) - sin(6pi/7))/3) (t5 - t3)
  1584 + const double b7_real = (-(int)sign) * tau7 * (t5_real - t3_real);
  1585 + const double b7_imag = (-(int)sign) * tau7 * (t5_imag - t3_imag);
  1586 +
  1587 + // b8 = sign * ((sin(2pi/7) + sin(4pi/7) + 2sin(6pi/7))/3) (t3 - t1) */
  1588 + const double b8_real = (-(int)sign) * tau8 * (t3_real - t1_real);
  1589 + const double b8_imag = (-(int)sign) * tau8 * (t3_imag - t1_imag);
  1590 +
  1591 + // T0 = b0 + b1
  1592 + const double T0_real = b0_real + b1_real;
  1593 + const double T0_imag = b0_imag + b1_imag;
  1594 +
  1595 + // T1 = b2 + b3
  1596 + const double T1_real = b2_real + b3_real;
  1597 + const double T1_imag = b2_imag + b3_imag;
  1598 +
  1599 + // T2 = b4 - b3
  1600 + const double T2_real = b4_real - b3_real;
  1601 + const double T2_imag = b4_imag - b3_imag;
  1602 +
  1603 + // T3 = -b2 - b4
  1604 + const double T3_real = -b2_real - b4_real;
  1605 + const double T3_imag = -b2_imag - b4_imag;
  1606 +
  1607 + // T4 = b6 + b7
  1608 + const double T4_real = b6_real + b7_real;
  1609 + const double T4_imag = b6_imag + b7_imag;
  1610 +
  1611 + // T5 = b8 - b7
  1612 + const double T5_real = b8_real - b7_real;
  1613 + const double T5_imag = b8_imag - b7_imag;
  1614 +
  1615 + // T6 = -b8 - b6
  1616 + const double T6_real = -b8_real - b6_real;
  1617 + const double T6_imag = -b8_imag - b6_imag;
  1618 +
  1619 + // T7 = T0 + T1
  1620 + const double T7_real = T0_real + T1_real;
  1621 + const double T7_imag = T0_imag + T1_imag;
  1622 +
  1623 + // T8 = T0 + T2
  1624 + const double T8_real = T0_real + T2_real;
  1625 + const double T8_imag = T0_imag + T2_imag;
  1626 +
  1627 + // T9 = T0 + T3
  1628 + const double T9_real = T0_real + T3_real;
  1629 + const double T9_imag = T0_imag + T3_imag;
  1630 +
  1631 + // T10 = T4 + b5
  1632 + const double T10_real = T4_real + b5_real;
  1633 + const double T10_imag = T4_imag + b5_imag;
  1634 +
  1635 + // T11 = T5 + b5
  1636 + const double T11_real = T5_real + b5_real;
  1637 + const double T11_imag = T5_imag + b5_imag;
  1638 +
  1639 + // T12 = T6 + b5
  1640 + const double T12_real = T6_real + b5_real;
  1641 + const double T12_imag = T6_imag + b5_imag;
  1642 +
  1643 + // x0 = b0
  1644 + const double x0_real = b0_real;
  1645 + const double x0_imag = b0_imag;
  1646 +
  1647 + // x1 = T7 - i T10
  1648 + const double x1_real = T7_real + T10_imag;
  1649 + const double x1_imag = T7_imag - T10_real;
  1650 +
  1651 + // x2 = T9 - i T12
  1652 + const double x2_real = T9_real + T12_imag;
  1653 + const double x2_imag = T9_imag - T12_real;
  1654 +
  1655 + // x3 = T8 + i T11
  1656 + const double x3_real = T8_real - T11_imag;
  1657 + const double x3_imag = T8_imag + T11_real;
  1658 +
  1659 + // x4 = T8 - i T11
  1660 + const double x4_real = T8_real + T11_imag;
  1661 + const double x4_imag = T8_imag - T11_real;
  1662 +
  1663 + // x5 = T9 + i T12
  1664 + const double x5_real = T9_real - T12_imag;
  1665 + const double x5_imag = T9_imag + T12_real;
  1666 +
  1667 + // x6 = T7 + i T10
  1668 + const double x6_real = T7_real - T10_imag;
  1669 + const double x6_imag = T7_imag + T10_real;
  1670 +
  1671 + // out0 = 1 * x0
  1672 + out[ostride*j].real(x0_real);
  1673 + out[ostride*j].imag(x0_imag);
  1674 +
  1675 + // out1 = w1 * x1
  1676 + out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
  1677 + out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
  1678 +
  1679 + // out2 = w2 * x2
  1680 + out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
  1681 + out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
  1682 +
  1683 + // out3 = w3 * x3
  1684 + out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
  1685 + out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
  1686 +
  1687 + // out4 = w4 * x4
  1688 + out[ostride*(j+4*p_1)].real(w4_real * x4_real - w4_imag * x4_imag);
  1689 + out[ostride*(j+4*p_1)].imag(w4_real * x4_imag + w4_imag * x4_real);
  1690 +
  1691 + // out5 = w5 * x5
  1692 + out[ostride*(j+5*p_1)].real(w5_real * x5_real - w5_imag * x5_imag);
  1693 + out[ostride*(j+5*p_1)].imag(w5_real * x5_imag + w5_imag * x5_real);
  1694 +
  1695 + // out6 = w6 * x6
  1696 + out[ostride*(j+6*p_1)].real(w6_real * x6_real - w6_imag * x6_imag);
  1697 + out[ostride*(j+6*p_1)].imag(w6_real * x6_imag + w6_imag * x6_real);
  1698 +
  1699 + } // endfor: k1
  1700 +
  1701 + } // endfor: k
  1702 +
  1703 + // Return
  1704 + return;
  1705 +}
  1706 +
  1707 +
  1708 +/***********************************************************************//**
  1709 + * @brief Compute FFT for arbitrary factor
  1710 + *
  1711 + * @param[in] in Input array.
  1712 + * @param[in] istride Step size when traversing input array.
  1713 + * @param[in] out Output array.
  1714 + * @param[in] ostride Step size when traversing output array.
  1715 + * @param[in] wavetable Trigonometric coefficients.
  1716 + * @param[in] sign Forward (-1) or backward (+1).
  1717 + * @param[in] factor Factor.
  1718 + * @param[in] product ...
  1719 + * @param[in] n Logical array length.
  1720 + * @param[in] index Start index of trigonometric coefficients.
  1721 + ***************************************************************************/
  1722 +void GFft::factorn(std::vector<std::complex<double> >& in,
  1723 + const int& istride,
  1724 + std::vector<std::complex<double> >& out,
  1725 + const int& ostride,
  1726 + const GFftWavetable& wavetable,
  1727 + const int& sign,
  1728 + const int& factor,
  1729 + const int& product,
  1730 + const int& n,
  1731 + const int& index)
  1732 +{
  1733 + // Compute ...
  1734 + const int m = n / factor;
  1735 + const int q = n / product;
  1736 + const int p_1 = product / factor;
  1737 + const int jump = (factor - 1) * p_1;
  1738 +
  1739 + // ...
  1740 + for (int i = 0; i < m; ++i) {
  1741 + out[ostride*i] = in[istride*i];
  1742 + }
  1743 +
  1744 + // ...
  1745 + for (int e = 1; e < (factor - 1) / 2 + 1; ++e) {
  1746 + for (int i = 0; i < m; ++i) {
  1747 + const int idx = i + e * m;
  1748 + const int idxc = i + (factor - e) * m;
  1749 + out[ostride*idx] = in[istride*idx] + in[istride*idxc];
  1750 + out[ostride*idxc] = in[istride*idx] - in[istride*idxc];
  1751 + }
  1752 + }
  1753 +
  1754 + // e = 0
  1755 + for (int i = 0; i < m; ++i) {
  1756 + in[istride*i] = out[ostride*i];
  1757 + }
  1758 + for (int e = 1; e < (factor - 1) / 2 + 1; ++e) {
  1759 + for (int i = 0; i < m; ++i) {
  1760 + const int idx = i + e * m;
  1761 + in[istride*i] += out[ostride*idx];
  1762 + }
  1763 + }
  1764 +
  1765 + // ...
  1766 + for (int e = 1; e < (factor-1)/2 + 1; ++e) {
  1767 +
  1768 + //
  1769 + int idx = e * q;
  1770 + const int idx_step = e * q;
  1771 +
  1772 + //
  1773 + double w_real, w_imag;
  1774 +
  1775 + //
  1776 + const int em = e * m;
  1777 + const int ecm = (factor - e) * m;
  1778 +
  1779 + // ...
  1780 + for (int i = 0; i < m; ++i) {
  1781 + in[istride*(i+em)] = out[ostride*i];
  1782 + in[istride*(i+ecm)] = out[ostride*i];
  1783 + }
  1784 +
  1785 + // ...
  1786 + for (int e1 = 1; e1 < (factor - 1) / 2 + 1; ++e1) {
  1787 +
  1788 + //
  1789 + if (idx == 0) {
  1790 + w_real = 1.0;
  1791 + w_imag = 0.0;
  1792 + }
  1793 + else {
  1794 +
  1795 + // Compute indices
  1796 + int twiddle = index + idx - 1;
  1797 +
  1798 + // Set trigonometric coefficients for forward transform
  1799 + if (sign == -1) {
  1800 + w_real = wavetable[twiddle].real();
  1801 + w_imag = wavetable[twiddle].imag();
  1802 + }
  1803 + // ... otherwise set trigonometric coefficients for backward
  1804 + // tranform: w -> conjugate(w)
  1805 + else {
  1806 + w_real = wavetable[twiddle].real();
  1807 + w_imag = -wavetable[twiddle].imag();
  1808 + }
  1809 + }
  1810 +
  1811 + // Loop over ...
  1812 + for (int i = 0; i < m; ++i) {
  1813 +
  1814 + //
  1815 + const double xp_real = out[ostride*(i + e1 * m)].real();
  1816 + const double xp_imag = out[ostride*(i + e1 * m)].imag();
  1817 + const double xm_real = out[ostride*(i + (factor-e1)*m)].real();
  1818 + const double xm_imag = out[ostride*(i + (factor-e1)*m)].imag();
  1819 +
  1820 + //
  1821 + const double ap = w_real * xp_real;
  1822 + const double am = w_imag * xm_imag;
  1823 +
  1824 + //
  1825 + double sum_real = ap - am;
  1826 + double sumc_real = ap + am;
  1827 +
  1828 + //
  1829 + const double bp = w_real * xp_imag;
  1830 + const double bm = w_imag * xm_real;
  1831 +
  1832 + //
  1833 + double sum_imag = bp + bm;
  1834 + double sumc_imag = bp - bm;
  1835 +
  1836 + //
  1837 + in[istride*(i+em)].real(in[istride*(i+em)].real() + sum_real);
  1838 + in[istride*(i+em)].imag(in[istride*(i+em)].imag() + sum_imag);
  1839 + in[istride*(i+ecm)].real(in[istride*(i+ecm)].real() + sumc_real);
  1840 + in[istride*(i+ecm)].imag(in[istride*(i+ecm)].imag() + sumc_imag);
  1841 +
  1842 + } // endfor: i
  1843 +
  1844 + // Increment
  1845 + idx += idx_step ;
  1846 + idx %= factor * q ;
  1847 +
  1848 + } // endfor: e1
  1849 +
  1850 + } // endfor: e
  1851 +
  1852 + // k = 0
  1853 + for (int k1 = 0; k1 < p_1; ++k1) {
  1854 + out[ostride*k1] = in[istride*k1];
  1855 + }
  1856 +
  1857 + // k > 0
  1858 + for (int e1 = 1; e1 < factor; ++e1) {
  1859 + for (int k1 = 0; k1 < p_1; ++k1) {
  1860 + out[ostride*(k1 + e1 * p_1)] = in[istride*(k1 + e1 * m)];
  1861 + }
  1862 + }
  1863 +
  1864 + // e = 0
  1865 + for (int k = 1, i = p_1, j = product; k < q; ++k, j += jump) {
  1866 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  1867 + out[ostride*j] = in[istride*i];
  1868 + }
  1869 + }
  1870 +
  1871 + // e > 0
  1872 + for (int k = 1, i = p_1, j = product; k < q; ++k, j += jump) {
  1873 +
  1874 + for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
  1875 +
  1876 + for (int e1 = 1; e1 < factor; ++e1) {
  1877 +
  1878 + // Get x
  1879 + double x_real = in[istride*(i + e1 * m)].real();
  1880 + double x_imag = in[istride*(i + e1 * m)].imag();
  1881 +
  1882 + // Get w
  1883 + double w_real, w_imag;
  1884 +
  1885 + // Compute indices
  1886 + int twiddle = index + (e1-1)*q + k-1;
  1887 +
  1888 + // Set trigonometric coefficients for forward transform
  1889 + if (sign == -1) {
  1890 + w_real = wavetable[twiddle].real();
  1891 + w_imag = wavetable[twiddle].imag();
  1892 + }
  1893 + // ... otherwise set trigonometric coefficients for backward
  1894 + // tranform: w -> conjugate(w)
  1895 + else {
  1896 + w_real = wavetable[twiddle].real();
  1897 + w_imag = -wavetable[twiddle].imag();
  1898 + }
  1899 +
  1900 + //
  1901 + out[ostride*(j + e1 * p_1)].real(w_real * x_real - w_imag * x_imag);
  1902 + out[ostride*(j + e1 * p_1)].imag(w_real * x_imag + w_imag * x_real);
  1903 +
  1904 + } // endfor: e1
  1905 +
  1906 + } // endfor: k1
  1907 +
  1908 + } // endfor: k
  1909 +
  1910 + // Return
  1911 + return;
  1912 +}
... ...
src/numerics/GFftWavetable.cpp 0 → 100644
  View file @768912b
  1 +/***************************************************************************
  2 + * GFftWavetable.cpp - Lookup table class for Fast Fourier transformation *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GFftWavetable.cpp
  23 + * @brief Lookup table class implementation for Fast Fourier transformation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GMath.hpp"
  33 +#include "GException.hpp"
  34 +#include "GFftWavetable.hpp"
  35 +
  36 +
  37 +/* __ Method name definitions ____________________________________________ */
  38 +#define G_FACTOR "GFftWavetable::factor(int&)"
  39 +#define G_SET_MEMBERS "GFftWavetable::set_members(int&)"
  40 +#define G_SET_FACTORS "GFftWavetable::set_factors(int&)"
  41 +
  42 +
  43 +/*==========================================================================
  44 + = =
  45 + = Constructors/destructors =
  46 + = =
  47 + ==========================================================================*/
  48 +
  49 +/***********************************************************************//**
  50 + * @brief Void constructor
  51 + ***************************************************************************/
  52 +GFftWavetable::GFftWavetable(void)
  53 +{
  54 + // Initialise class members
  55 + init_members();
  56 +
  57 + // Return
  58 + return;
  59 +}
  60 +
  61 +
  62 +/***********************************************************************//**
  63 + * @brief Wave table constructor
  64 + *
  65 + * @param[in] size Array size.
  66 + *
  67 + * Constructs the lookup table for Fast Fourier Transformation for a given
  68 + * array @p size.
  69 + ***************************************************************************/
  70 +GFftWavetable::GFftWavetable(const int& size)
  71 +{
  72 + // Initialise class members
  73 + init_members();
  74 +
  75 + // Set members
  76 + set_members(size);
  77 +
  78 + // Return
  79 + return;
  80 +}
  81 +
  82 +
  83 +/***********************************************************************//**
  84 + * @brief Copy constructor
  85 + *
  86 + * @param[in] wavetable Lookup table for Fast Fourier Transform.
  87 + ***************************************************************************/
  88 +GFftWavetable::GFftWavetable(const GFftWavetable& wavetable)
  89 +{
  90 + // Initialise class members
  91 + init_members();
  92 +
  93 + // Copy members
  94 + copy_members(wavetable);
  95 +
  96 + // Return
  97 + return;
  98 +}
  99 +
  100 +
  101 +/***********************************************************************//**
  102 + * @brief Destructor
  103 + ***************************************************************************/
  104 +GFftWavetable::~GFftWavetable(void)
  105 +{
  106 + // Free members
  107 + free_members();
  108 +
  109 + // Return
  110 + return;
  111 +}
  112 +
  113 +
  114 +/*==========================================================================
  115 + = =
  116 + = Operators =
  117 + = =
  118 + ==========================================================================*/
  119 +
  120 +/***********************************************************************//**
  121 + * @brief Assignment operator
  122 + *
  123 + * @param[in] wavetable Lookup table for Fast Fourier Transform.
  124 + * @return Lookup table for Fast Fourier Transform.
  125 + ***************************************************************************/
  126 +GFftWavetable& GFftWavetable::operator=(const GFftWavetable& wavetable)
  127 +{
  128 + // Execute only if object is not identical
  129 + if (this != &wavetable) {
  130 +
  131 + // Free members
  132 + free_members();
  133 +
  134 + // Initialise private members
  135 + init_members();
  136 +
  137 + // Copy members
  138 + copy_members(wavetable);
  139 +
  140 + } // endif: object was not identical
  141 +
  142 + // Return this object
  143 + return *this;
  144 +}
  145 +
  146 +
  147 +/*==========================================================================
  148 + = =
  149 + = Public methods =
  150 + = =
  151 + ==========================================================================*/
  152 +
  153 +/***********************************************************************//**
  154 + * @brief Clear lookup table for Fast Fourier Transform
  155 + ***************************************************************************/
  156 +void GFftWavetable::clear(void)
  157 +{
  158 + // Free members
  159 + free_members();
  160 +
  161 + // Initialise private members
  162 + init_members();
  163 +
  164 + // Return
  165 + return;
  166 +}
  167 +
  168 +
  169 +/***********************************************************************//**
  170 + * @brief Clone lookup table for Fast Fourier Transform
  171 + *
  172 + * @return Pointer to deep copy of lookup table for Fast Fourier Transform.
  173 + ***************************************************************************/
  174 +GFftWavetable* GFftWavetable::clone(void) const
  175 +{
  176 + // Clone FFT
  177 + return new GFftWavetable(*this);
  178 +}
  179 +
  180 +
  181 +/***********************************************************************//**
  182 + * @brief Return factorisation factor
  183 + *
  184 + * @param[in] index Factorisation index [0,...,factors()-1].
  185 + * @return Factorisation factor.
  186 + *
  187 + * @exception GException::out_of_range
  188 + * Index is outside valid range.
  189 + *
  190 + * Returns factorisation factor.
  191 + ***************************************************************************/
  192 +int GFftWavetable::factor(const int& index) const
  193 +{
  194 + // Throw an exception if index is out of range
  195 + if (index < 0 || index >= factors()) {
  196 + throw GException::out_of_range(G_FACTOR, "Factorisation factor index",
  197 + index, factors());
  198 + }
  199 +
  200 + // Return factorisation factor
  201 + return (m_factors[index]);
  202 +}
  203 +
  204 +
  205 +/***********************************************************************//**
  206 + * @brief Print lookup table for Fast Fourier Transform information
  207 + *
  208 + * @param[in] chatter Chattiness.
  209 + * @return String containing lookup table for Fast Fourier Transform
  210 + * information.
  211 + ***************************************************************************/
  212 +std::string GFftWavetable::print(const GChatter& chatter) const
  213 +{
  214 + // Initialise result string
  215 + std::string result;
  216 +
  217 + // Continue only if chatter is not silent
  218 + if (chatter != SILENT) {
  219 +
  220 + // Append header
  221 + result.append("=== GFftWavetable ===");
  222 +
  223 + // Append information
  224 + result.append("\n"+gammalib::parformat("Coefficients"));
  225 + result.append(gammalib::str(size()));
  226 + result.append("\n"+gammalib::parformat("Factorisation"));
  227 + if (factors() > 0) {
  228 + for (int i = 0; i < factors(); ++i) {
  229 + if (i > 0) {
  230 + result.append("*");
  231 + }
  232 + result.append(gammalib::str(factor(i)));
  233 + }
  234 + }
  235 + else {
  236 + result.append("none");
  237 + }
  238 +
  239 + } // endif: chatter was not silent
  240 +
  241 + // Return result
  242 + return result;
  243 +}
  244 +
  245 +
  246 +/*==========================================================================
  247 + = =
  248 + = Private methods =
  249 + = =
  250 + ==========================================================================*/
  251 +
  252 +/***********************************************************************//**
  253 + * @brief Initialise class members
  254 + ***************************************************************************/
  255 +void GFftWavetable::init_members(void)
  256 +{
  257 + // Initialise members
  258 + m_factors.clear();
  259 + m_twiddle.clear();
  260 + m_trig.clear();
  261 +
  262 + // Return
  263 + return;
  264 +}
  265 +
  266 +
  267 +/***********************************************************************//**
  268 + * @brief Copy class members
  269 + *
  270 + * @param[in] wavetable Lookup table for Fast Fourier Transform.
  271 + ***************************************************************************/
  272 +void GFftWavetable::copy_members(const GFftWavetable& wavetable)
  273 +{
  274 + // Copy members
  275 + m_factors = wavetable.m_factors;
  276 + m_twiddle = wavetable.m_twiddle;
  277 + m_trig = wavetable.m_trig;
  278 +
  279 + // Return
  280 + return;
  281 +}
  282 +
  283 +
  284 +/***********************************************************************//**
  285 + * @brief Delete class members
  286 + ***************************************************************************/
  287 +void GFftWavetable::free_members(void)
  288 +{
  289 + // Return
  290 + return;
  291 +}
  292 +
  293 +
  294 +/***********************************************************************//**
  295 + * @brief Set wavetable
  296 + *
  297 + * @param[in] n Length.
  298 + *
  299 + * Computes the coefficients
  300 + *
  301 + * \f[
  302 + * \frac{-i 2 \pi j k}{n}
  303 + * \f]
  304 + *
  305 + * of the discrete Fourier transformation
  306 + *
  307 + * \f[
  308 + * x_j = \sum_{k=0}^{n-1} z_k \exp \left( \frac{-i 2 \pi j k}{n} \right)
  309 + * \f]
  310 + *
  311 + * The method is a C++ implementation of the gsl_fft_complex_wavetable_alloc()
  312 + * function of the GSL library (version 2.2.1), defined in the file
  313 + * fft/c_init.c.
  314 + ***************************************************************************/
  315 +void GFftWavetable::set_members(const int& n)
  316 +{
  317 + // Initialise trigonometric lookup table
  318 + init_members();
  319 +
  320 + // Make sure that length is positive
  321 + if (n < 1) {
  322 + std::string msg = "Invalid array length "+gammalib::str(n)+". Array "
  323 + "needs to be of positive length.";
  324 + throw GException::invalid_argument(G_SET_MEMBERS, msg);
  325 + }
  326 +
  327 + // Compute factorisation
  328 + set_factors(n);
  329 +
  330 + // Compute delta theta angle (2 pi / N)
  331 + double d_theta = -gammalib::twopi / double(n);
  332 +
  333 + // Initialise factorisation product
  334 + int product = 1;
  335 +
  336 + // Loop over all factors in the factorisation
  337 + for (int i = 0; i < m_factors.size(); ++i) {
  338 +
  339 + // Set start index of the current factor
  340 + m_twiddle.push_back(m_trig.size());
  341 +
  342 + // Store current factorisation product
  343 + int product_1 = product;
  344 +
  345 + // Compute next factorisation product
  346 + product *= m_factors[i];
  347 +
  348 + // Compute number of coefficients for current factor
  349 + int q = n / product;
  350 +/*
  351 +std::cout << "i=" << i << " q=" << q << " factor=" << m_factors[i];
  352 +std::cout << " product=" << product << " n=" << n << std::endl;
  353 +*/
  354 + // Loop over all output coefficients
  355 + for (int j = 1; j < m_factors[i]; ++j) {
  356 +
  357 + // Loop over all coefficients
  358 + for (int k = 0, m = 0; k < q; ++k) {
  359 +
  360 + // Compute theta value
  361 + m = m + j * product_1;
  362 + m = m % n;
  363 + double theta = d_theta * m; /* d_theta*j*k*p_(i-1) */
  364 +/*
  365 +std::cout << " j=" << j << " k=" << k << " theta=" << theta << std::endl;
  366 +*/
  367 + // Compute trigonometric lookup table value
  368 + std::complex<double> value(std::cos(theta), std::sin(theta));
  369 +
  370 + // Push back value on lookup table
  371 + m_trig.push_back(value);
  372 +
  373 + } // endfor: loop over all elements for a given factor
  374 +
  375 + } // endfor
  376 +
  377 + } // endfor: computed ...
  378 +
  379 + // Return
  380 + return;
  381 +}
  382 +
  383 +
  384 +/***********************************************************************//**
  385 + * @brief Compute FFT factorisation
  386 + *
  387 + * @param[in] n Length of array to be transformed.
  388 + *
  389 + * @exception GException::invalid_argument
  390 + * Non positive array length specified
  391 + * @expection GException::invalid_value
  392 + * Product of factorisation factors are not equal to array length
  393 + *
  394 + * Computes the factorisation
  395 + *
  396 + * \f[
  397 + * n = \prod_i p_i
  398 + * \f]
  399 + *
  400 + * of an array of length @p n.
  401 + *
  402 + * The method is a C++ implementation of the fft_complex_factorize() and the
  403 + * fft_factorize() functions of the GSL library (version 2.2.1), defined in
  404 + * the file fft/factorize.c.
  405 + ***************************************************************************/
  406 +void GFftWavetable::set_factors(const int& n)
  407 +{
  408 + // Make sure that length is positive
  409 + if (n < 1) {
  410 + std::string msg = "Invalid array length "+gammalib::str(n)+". Array "
  411 + "needs to be of positive length.";
  412 + throw GException::invalid_argument(G_SET_FACTORS, msg);
  413 + }
  414 +
  415 + // Set FFT factorisations. Other factors can be added here if their
  416 + // transform modules are implemented.
  417 + std::vector<int> subtransforms;
  418 + subtransforms.push_back(7);
  419 + subtransforms.push_back(6);
  420 + subtransforms.push_back(5);
  421 + subtransforms.push_back(4);
  422 + subtransforms.push_back(3);
  423 + subtransforms.push_back(2);
  424 +
  425 + // If array has unit length then set the first and only factor to 1 ...
  426 + if (n == 1) {
  427 + m_factors.push_back(1);
  428 + }
  429 +
  430 + // ... otherwise
  431 + else {
  432 +
  433 + // Save array length
  434 + int ntest = n;
  435 +
  436 + // Deal with the implemented factors first
  437 + for (int i = 0; i < subtransforms.size() && ntest > 1; ++i) {
  438 + int factor = subtransforms[i];
  439 + while ((ntest % factor) == 0) {
  440 + ntest /= factor;
  441 + m_factors.push_back(factor);
  442 + }
  443 + }
  444 + /*
  445 + int i = 0;
  446 + while (subtransforms[i] && ntest != 1) {
  447 + int factor = subtransforms[i];
  448 + while ((ntest % factor) == 0) {
  449 + ntest /= factor;
  450 + m_factors.push_back(factor);
  451 + }
  452 + i++;
  453 + }
  454 + */
  455 +
  456 + // Deal with any other even prime factors (there is only one)
  457 + int factor = 2;
  458 + while ((ntest % factor) == 0 && (ntest != 1)) {
  459 + ntest /= factor;
  460 + m_factors.push_back(factor);
  461 + }
  462 +
  463 + // Deal with any other odd prime factors
  464 + factor = 3;
  465 + while (ntest != 1) {
  466 + while ((ntest % factor) != 0) {
  467 + factor += 2;
  468 + }
  469 + ntest /= factor;
  470 + m_factors.push_back(factor);
  471 + }
  472 +
  473 + // Check that the factorization is correct
  474 + int product = 1;
  475 + for (int i = 0; i < m_factors.size(); ++i) {
  476 + product *= m_factors[i];
  477 + }
  478 + if (product != n) {
  479 + std::string msg = "Product "+gammalib::str(product)+" of "
  480 + "factorisation factors differs from array "
  481 + "length "+gammalib::str(n)+".";
  482 + throw GException::invalid_value(G_SET_FACTORS, msg);
  483 + }
  484 +
  485 + } // endelse: array length is larger than 1
  486 +
  487 + // Return
  488 + return;
  489 +}
... ...
src/numerics/GNdarray.cpp
  View file @768912b
... ... @@ -159,31 +159,38 @@ GNdarray::GNdarray(const int&amp; nx, const int&amp; ny, const int&amp; nz)
159 159 *
160 160 * @param[in] n Array dimensions.
161 161 *
162   - * Constructs a n-dimensional array.
  162 + * Constructs a n-dimensional array. If the array dimension vector @p n is
  163 + * empty an empty n-dimensional array is constructed (similar to the void
  164 + * constructor).
163 165 ***************************************************************************/
164 166 GNdarray::GNdarray(const std::vector<int>& n)
165 167 {
166 168 // Initialise class members
167 169 init_members();
168 170  
169   - // Compute size of array
170   - int size = n[0];
171   - for (int i = 1; i < n.size(); ++i) {
172   - size *= n[i];
173   - }
  171 + // Continue only if there are dimensions in vector
  172 + if (n.size() > 0) {
174 173  
175   - // Set shape
176   - m_shape = n;
  174 + // Compute size of array
  175 + int size = n[0];
  176 + for (int i = 1; i < n.size(); ++i) {
  177 + size *= n[i];
  178 + }
177 179  
178   - // Set strides
179   - int stride = 1;
180   - for (int i = 0; i < n.size(); ++i) {
181   - m_strides.push_back(stride);
182   - stride *= n[i];
183   - }
  180 + // Set shape
  181 + m_shape = n;
184 182  
185   - // Initialise data array with zeros
186   - m_data.assign(size, 0.0);
  183 + // Set strides
  184 + int stride = 1;
  185 + for (int i = 0; i < n.size(); ++i) {
  186 + m_strides.push_back(stride);
  187 + stride *= n[i];
  188 + }
  189 +
  190 + // Initialise data array with zeros
  191 + m_data.assign(size, 0.0);
  192 +
  193 + } // endif: there were dimensions in vector
187 194  
188 195 // Return
189 196 return;
... ...
src/numerics/Makefile.am
  View file @768912b
... ... @@ -30,6 +30,8 @@ sources = GIntegral.cpp \
30 30 GFunction.cpp \
31 31 GMath.cpp \
32 32 GNdarray.cpp \
  33 + GFft.cpp \
  34 + GFftWavetable.cpp \
33 35 GException_numerics.cpp
34 36  
35 37 # Build libtool library
... ...
src/support/GTools.cpp
  View file @768912b
... ... @@ -549,6 +549,51 @@ std::string gammalib::str(const double&amp; value, const int&amp; precision)
549 549  
550 550  
551 551 /***********************************************************************//**
  552 + * @brief Convert complex value into string
  553 + *
  554 + * @param[in] value Complex value to be converted into string.
  555 + * @param[in] precision Floating point precision.
  556 + * @return String with complex value.
  557 + *
  558 + * Converts a complex value into a string. Any positive @p precision argument
  559 + * specifies the exact number of digits after the comma.
  560 + ***************************************************************************/
  561 +std::string gammalib::str(const std::complex<double>& value,
  562 + const int& precision)
  563 +{
  564 + // Allocate output stream
  565 + std::ostringstream s_value;
  566 +
  567 + // If specified then set the requested fixed point precision. Otherwise
  568 + // use a precision that should be sufficient for floating point values.
  569 + if (precision > 0) {
  570 + s_value.precision(precision);
  571 + s_value.setf(std::ios::fixed, std::ios::floatfield);
  572 + }
  573 + else {
  574 + s_value.precision(15);
  575 + }
  576 +
  577 + // Put double precision floating point value in stream
  578 + s_value << value.real();
  579 + if (value.imag() < 0.0) {
  580 + s_value << "-";
  581 + }
  582 + else {
  583 + s_value << "+";
  584 + }
  585 + s_value << std::abs(value.imag()) << "j";
  586 +
  587 + // Convert to a string
  588 + std::string result = s_value.str();
  589 +
  590 + // Return result
  591 + return result;
  592 +}
  593 +
  594 +
  595 +
  596 +/***********************************************************************//**
552 597 * @brief Convert string to C string
553 598 *
554 599 * @param[in] arg String to be converted.
... ...
src/test/GTestSuite.cpp
  View file @768912b
... ... @@ -548,6 +548,103 @@ void GTestSuite::test_value(const double&amp; value,
548 548  
549 549  
550 550 /***********************************************************************//**
  551 + * @brief Test a complex value
  552 + *
  553 + * @param[in] value Complex value to test.
  554 + * @param[in] expected Expected complex value.
  555 + * @param[in] name Test case name.
  556 + * @param[in] message Test case message.
  557 + *
  558 + * Test if the @p value is equal to the @p expected value within a relative
  559 + * precision of 1.0e-7.
  560 + ***************************************************************************/
  561 +void GTestSuite::test_value(const std::complex<double>& value,
  562 + const std::complex<double>& expected,
  563 + const std::string& name,
  564 + const std::string& message)
  565 +{
  566 + // Compute precision
  567 + double eps = (expected != 0.0) ? 1.0e-7 * std::abs(expected) : 1.0e-7;
  568 +
  569 + // Test double precision value
  570 + test_value(value, expected, eps, name, message);
  571 +
  572 + // Return
  573 + return;
  574 +}
  575 +
  576 +
  577 +/***********************************************************************//**
  578 + * @brief Test a complex value
  579 + *
  580 + * @param[in] value Complex value to test.
  581 + * @param[in] expected Expected complex value.
  582 + * @param[in] eps Precision of the test.
  583 + * @param[in] name Test case name.
  584 + * @param[in] message Test case message.
  585 + *
  586 + * Test if the value is comprised in the interval
  587 + * [expected-eps, expected+eps].
  588 + ***************************************************************************/
  589 +void GTestSuite::test_value(const std::complex<double>& value,
  590 + const std::complex<double>& expected,
  591 + const double& eps,
  592 + const std::string& name,
  593 + const std::string& message)
  594 +{
  595 + // Set test case name. If no name is specify then build the name from
  596 + // the actual test parameters.
  597 + std::string formated_name;
  598 + if (name != "") {
  599 + formated_name = format_name(name);
  600 + }
  601 + else {
  602 + formated_name = format_name("Test if "+gammalib::str(value)+
  603 + " is comprised within "+
  604 + gammalib::str(expected)+" +/- "+
  605 + gammalib::str(eps));
  606 + }
  607 +
  608 + // Create a test case of failure type
  609 + GTestCase* testcase = new GTestCase(GTestCase::FAIL_TEST, formated_name);
  610 +
  611 + // If value is not between in interval [expected-eps, expected+eps]
  612 + // then signal test as failed and increment the number of failures
  613 + if ((value.real() > expected.real() + eps) ||
  614 + (value.real() < expected.real() - eps) ||
  615 + (value.imag() > expected.imag() + eps) ||
  616 + (value.imag() < expected.imag() - eps)) {
  617 + testcase->has_passed(false);
  618 + m_failures++;
  619 + }
  620 +
  621 + // If no message is specified then build message from test result
  622 + std::string formated_message;
  623 + if (message != "") {
  624 + formated_message = message;
  625 + }
  626 + else {
  627 + formated_message = "Value "+gammalib::str(value)+" not within "+
  628 + gammalib::str(expected)+" +/- "+gammalib::str(eps)+
  629 + " (value-expected = "+gammalib::str(value-expected)+
  630 + ").";
  631 + }
  632 +
  633 + // Set message
  634 + testcase->message(formated_message);
  635 +
  636 + // Log the result (".","F" or, "E")
  637 + std::cout << testcase->print();
  638 +
  639 + // Add test case to test suite
  640 + m_tests.push_back(testcase);
  641 +
  642 + // Return
  643 + return;
  644 +}
  645 +
  646 +
  647 +/***********************************************************************//**
551 648 * @brief Test a string value
552 649 *
553 650 * @param[in] value String value to test.
... ...
test/test_GNumerics.cpp
  View file @768912b
... ... @@ -52,6 +52,8 @@ void TestGNumerics::set(void)
52 52 // Append tests
53 53 append(static_cast<pfunction>(&TestGNumerics::test_ndarray),
54 54 "Test GNdarray");
  55 + append(static_cast<pfunction>(&TestGNumerics::test_fft),
  56 + "Test GFft");
55 57 append(static_cast<pfunction>(&TestGNumerics::test_integral),
56 58 "Test GIntegral");
57 59 append(static_cast<pfunction>(&TestGNumerics::test_romberg_integration),
... ... @@ -276,6 +278,78 @@ void TestGNumerics::test_ndarray(void)
276 278  
277 279  
278 280 /***********************************************************************//**
  281 + * @brief Test Fast Fourier Transform
  282 + ***************************************************************************/
  283 +void TestGNumerics::test_fft(void)
  284 +{
  285 + // Loop over array lengths
  286 + for (int n = 2; n < 12; ++n) {
  287 +
  288 + // Skip n=8,9,10 since they can be factorised
  289 + if (n > 7 && n < 11) {
  290 + continue;
  291 + }
  292 +
  293 + // Allocate 1-dimensional array
  294 + GNdarray array(n);
  295 +
  296 + // Set one elements
  297 + array(n/2) = 1.0;
  298 +
  299 + // Allocate FFT
  300 + GFft fft(array);
  301 +
  302 + // Backward transformation
  303 + GNdarray back = fft.backward();
  304 +
  305 + // Check FFT results
  306 + check_fft(array, fft, back);
  307 +
  308 + } // endfor: looped over array lengths
  309 +
  310 + // Exit test
  311 + return;
  312 +}
  313 +
  314 +
  315 +/***********************************************************************//**
  316 + * @brief Check Fast Fourier Transform
  317 + *
  318 + * @param[in] array Input array.
  319 + * @param[in] fft Fast Fourier Transform of input array.
  320 + * @param[in] back Backward transform of FFT (should be equal to input array)
  321 + ***************************************************************************/
  322 +void TestGNumerics::check_fft(const GNdarray& array,
  323 + const GFft& fft,
  324 + const GNdarray& back)
  325 +{
  326 + // Get normalisation factor
  327 + double norm = 1.0 / double(array.size());
  328 +
  329 + // Compute reference value and compare with FFT result
  330 + for (int k = 0; k < array.size(); ++k) {
  331 + double ref_real = 0.0;
  332 + double ref_imag = 0.0;
  333 + for (int n = 0; n < array.size(); ++n) {
  334 + double x = -gammalib::twopi * k * n * norm;
  335 + ref_real += array(n) * std::cos(x);
  336 + ref_imag += array(n) * std::sin(x);
  337 + }
  338 + std::complex<double> ref(ref_real, ref_imag);
  339 + test_value(fft(k), ref);
  340 + }
  341 +
  342 + // Test backward transformation
  343 + for (int k = 0; k < array.size(); ++k) {
  344 + test_value(back(k), array(k));
  345 + }
  346 +
  347 + // Exit test
  348 + return;
  349 +}
  350 +
  351 +
  352 +/***********************************************************************//**
279 353 * @brief Test model parameter handling
280 354 ***************************************************************************/
281 355 void TestGNumerics::test_integral(void)
... ...
test/test_GNumerics.hpp
  View file @768912b
... ... @@ -69,6 +69,10 @@ public:
69 69 virtual TestGNumerics* clone(void) const;
70 70 virtual std::string classname(void) const { return "TestGNumerics"; }
71 71 void test_ndarray(void);
  72 + void test_fft(void);
  73 + void check_fft(const GNdarray& array,
  74 + const GFft& fft,
  75 + const GNdarray& back);
72 76 void test_integral(void);
73 77 void test_romberg_integration(void);
74 78 void test_adaptive_simpson_integration(void);
... ...