GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GFft.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GFft.hpp - Fast Fourier transformation class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2016-2018 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 
40 
41 /***********************************************************************//**
42  * @class GFft
43  *
44  * @brief Fast Fourier Transformation class
45  *
46  * This class implements a Fast Fourier Transformation of an n-dimensional
47  * double precision floating point array.
48  *
49  * The class implementation is based on the Fast Fourier Transformation
50  * functions that are provided by the GNU Scientific Library (GSL)
51  * (version 2.2.1). The mixed-radix routines for complex numbers have been
52  * adopted that work for any array lengths. The mixed-radix routines are a
53  * reimplementation of the FFTPACK library of Paul Swarztrauber.
54  *
55  * The GSL is documented at https://www.gnu.org/software/gsl/manual
56  ***************************************************************************/
57 class GFft : public GBase {
58 
59 public:
60  // Constructors and destructors
61  GFft(void);
62  explicit GFft(const GNdarray& array);
63  GFft(const GFft& fft);
64  virtual ~GFft(void);
65 
66  // FFT access operators
67  std::complex<double>& operator()(const int& ix);
68  std::complex<double>& operator()(const int& ix, const int& iy);
69  const std::complex<double>& operator()(const int& ix) const;
70  const std::complex<double>& operator()(const int& ix, const int& iy) const;
71 
72  // Operators
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/=(const GFft& fft);
78  GFft operator-(void) const;
79 
80  // Methods
81  void clear(void);
82  GFft* clone(void) const;
83  std::string classname(void) const;
84  int dim(void) const;
85  int size(void) const;
86  const std::vector<int>& shape(void) const;
87  const std::vector<int>& strides(void) const;
88  void forward(const GNdarray& array);
89  GNdarray backward(void) const;
90  std::string print(const GChatter& chatter = NORMAL) const;
91 
92 protected:
93  // Protected methods
94  void init_members(void);
95  void copy_members(const GFft& fft);
96  void free_members(void);
97  void set_data(const GNdarray& array);
98  bool has_same_shape(const GFft& fft) const;
99  void require_same_shape(const std::string& method, const GFft& fft) const;
100 
101  // Low-level FFT methods
102  void transform(std::complex<double>* data,
103  const int& stride,
104  const int& n,
105  const GFftWavetable& wavetable,
106  const bool& forward = true);
107  void factor2(const std::complex<double>* in,
108  const int& istride,
109  std::complex<double>* out,
110  const int& ostride,
111  const GFftWavetable& wavetable,
112  const int& sign,
113  const int& product,
114  const int& n,
115  const int& index);
116  void factor3(const std::complex<double>* in,
117  const int& istride,
118  std::complex<double>* out,
119  const int& ostride,
120  const GFftWavetable& wavetable,
121  const int& sign,
122  const int& product,
123  const int& n,
124  const int& index);
125  void factor4(const std::complex<double>* in,
126  const int& istride,
127  std::complex<double>* out,
128  const int& ostride,
129  const GFftWavetable& wavetable,
130  const int& sign,
131  const int& product,
132  const int& n,
133  const int& index);
134  void factor5(const std::complex<double>* in,
135  const int& istride,
136  std::complex<double>* out,
137  const int& ostride,
138  const GFftWavetable& wavetable,
139  const int& sign,
140  const int& product,
141  const int& n,
142  const int& index);
143  void factor6(const std::complex<double>* in,
144  const int& istride,
145  std::complex<double>* out,
146  const int& ostride,
147  const GFftWavetable& wavetable,
148  const int& sign,
149  const int& product,
150  const int& n,
151  const int& index);
152  void factor7(const std::complex<double>* in,
153  const int& istride,
154  std::complex<double>* out,
155  const int& ostride,
156  const GFftWavetable& wavetable,
157  const int& sign,
158  const int& product,
159  const int& n,
160  const int& index);
161  void factorn(std::complex<double>* in,
162  const int& istride,
163  std::complex<double>* out,
164  const int& ostride,
165  const GFftWavetable& wavetable,
166  const int& sign,
167  const int& factor,
168  const int& product,
169  const int& n,
170  const int& index);
171  std::vector<std::complex<double> > get_w(const GFftWavetable& wavetable,
172  const int& index,
173  const int& k,
174  const int& q,
175  const int& n,
176  const int& sign) const;
177  std::complex<double> timesi(const std::complex<double>& value) const;
178 
179  // Protected members
180  int m_size; //!< Size of data array
181  std::complex<double>* m_data; //!< Pointer on array data
182  std::vector<int> m_shape; //!< Array dimensions
183  std::vector<int> m_strides; //!< Steps in each dimension
184  std::vector<GFftWavetable> m_wavetable; //!< Trigonometric coefficients
185 };
186 
187 
188 /***********************************************************************//**
189  * @brief Return class name
190  *
191  * @return String containing the class name ("GFft").
192  ***************************************************************************/
193 inline
194 std::string GFft::classname(void) const
195 {
196  return ("GFft");
197 }
198 
199 
200 /***********************************************************************//**
201  * @brief 1-dimensional FFT element access operator
202  *
203  * @param[in] ix Element index [0,...,shape(0)-1].
204  * @return Reference to FFT element.
205  ***************************************************************************/
206 inline
207 std::complex<double>& GFft::operator()(const int& ix)
208 {
209  // Return array element
210  return *(m_data + ix);
211 }
212 
213 
214 /***********************************************************************//**
215  * @brief 2-dimensional FFT element access operator
216  *
217  * @param[in] ix Index in first dimension [0,...,shape(0)-1].
218  * @param[in] iy Index in second dimension [0,...,shape(1)-1].
219  * @return Reference to FFT element.
220  ***************************************************************************/
221 inline
222 std::complex<double>& GFft::operator()(const int& ix, const int& iy)
223 {
224  // Return array element
225  return *(m_data + ix + m_strides[1]*iy);
226 }
227 
228 
229 /***********************************************************************//**
230  * @brief 1-dimensional FFT element access operator (const variant)
231  *
232  * @param[in] ix Element index [0,...,shape(0)-1].
233  * @return Const reference to FFT element.
234  ***************************************************************************/
235 inline
236 const std::complex<double>& GFft::operator()(const int& ix) const
237 {
238  // Return array element
239  return *(m_data + ix);
240 }
241 
242 
243 /***********************************************************************//**
244  * @brief 2-dimensional FFT element access operator (const variant)
245  *
246  * @param[in] ix Index in first dimension [0,...,shape(0)-1].
247  * @param[in] iy Index in second dimension [0,...,shape(1)-1].
248  * @return Const reference to FFT element.
249  ***************************************************************************/
250 inline
251 const std::complex<double>& GFft::operator()(const int& ix, const int& iy) const
252 {
253  // Return array element
254  return *(m_data + ix + m_strides[1]*iy);
255 }
256 
257 
258 
259 
260 /***********************************************************************//**
261  * @brief Return dimension of Fast Fourier Transformation
262  *
263  * @return Dimension of Fast Fourier Transformation.
264  *
265  * Returns the dimension of the Fast Fourier Transformation.
266  ***************************************************************************/
267 inline
268 int GFft::dim(void) const
269 {
270  return ((int)m_shape.size());
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Return number of elements in array
276  *
277  * @return Number of elements in array
278  *
279  * Returns the number of elements in the array.
280  ***************************************************************************/
281 inline
282 int GFft::size(void) const
283 {
284  return m_size;
285 }
286 
287 
288 /***********************************************************************//**
289  * @brief Return shape of array
290  *
291  * @return Shape of array
292  *
293  * Returns the shape of the array.
294  ***************************************************************************/
295 inline
296 const std::vector<int>& GFft::shape(void) const
297 {
298  return m_shape;
299 }
300 
301 
302 /***********************************************************************//**
303  * @brief Return strides of array
304  *
305  * @return Strides of array
306  *
307  * Returns the strides of the array.
308  ***************************************************************************/
309 inline
310 const std::vector<int>& GFft::strides(void) const
311 {
312  return m_strides;
313 }
314 
315 
316 /***********************************************************************//**
317  * @brief Return sum of two Fast Fourier Transformations
318  *
319  * @param[in] a First Fast Fourier Transformation.
320  * @param[in] b Second Fast Fourier Transformation.
321  * @return Sum of the Fast Fourier Transformations @p a and @p b.
322  *
323  * Returns the sum of the Fast Fourier Transformations @p a and @p b.
324  ***************************************************************************/
325 inline
326 GFft operator+(const GFft& a, const GFft& b)
327 {
328  GFft result = a;
329  result += b;
330  return result;
331 }
332 
333 
334 /***********************************************************************//**
335  * @brief Return difference of two Fast Fourier Transformations
336  *
337  * @param[in] a First Fast Fourier Transformation.
338  * @param[in] b Second Fast Fourier Transformation.
339  * @return Difference of the Fast Fourier Transformations @p a and @p b.
340  *
341  * Returns the difference of the Fast Fourier Transformations @p a and @p b.
342  ***************************************************************************/
343 inline
344 GFft operator-(const GFft& a, const GFft& b)
345 {
346  GFft result = a;
347  result -= b;
348  return result;
349 }
350 
351 
352 /***********************************************************************//**
353  * @brief Return product of two Fast Fourier Transformations
354  *
355  * @param[in] a First Fast Fourier Transformation.
356  * @param[in] b Second Fast Fourier Transformation.
357  * @return Product of the Fast Fourier Transformations @p a and @p b.
358  *
359  * Returns the product of the Fast Fourier Transformations @p a and @p b.
360  ***************************************************************************/
361 inline
362 GFft operator*(const GFft& a, const GFft& b)
363 {
364  GFft result = a;
365  result *= b;
366  return result;
367 }
368 
369 
370 /***********************************************************************//**
371  * @brief Return quotient of two Fast Fourier Transformations
372  *
373  * @param[in] a First Fast Fourier Transformation.
374  * @param[in] b Second Fast Fourier Transformation.
375  * @return Quotient of the Fast Fourier Transformations @p a and @p b.
376  *
377  * Returns the quotient of the Fast Fourier Transformations @p a and @p b.
378  ***************************************************************************/
379 inline
380 GFft operator/(const GFft& a, const GFft& b)
381 {
382  GFft result = a;
383  result /= b;
384  return result;
385 }
386 
387 
388 /***********************************************************************//**
389  * @brief Return complex value times i
390  *
391  * @param[in] value Complex value.
392  * @return Complex value times i.
393  *
394  * Returns complex value times i.
395  ***************************************************************************/
396 inline
397 std::complex<double> GFft::timesi(const std::complex<double>& value) const
398 {
399  std::complex<double> result(-value.imag(), value.real());
400  return result;
401 }
402 
403 #endif /* GFFT_HPP */
std::vector< int > m_shape
Array dimensions.
Definition: GFft.hpp:182
void factor3(const std::complex< double > *in, const int &istride, std::complex< double > *out, const int &ostride, const GFftWavetable &wavetable, const int &sign, const int &product, const int &n, const int &index)
Compute FFT for factor 3.
Definition: GFft.cpp:885
void factorn(std::complex< double > *in, const int &istride, std::complex< double > *out, const int &ostride, const GFftWavetable &wavetable, const int &sign, const int &factor, const int &product, const int &n, const int &index)
Compute FFT for arbitrary factor.
Definition: GFft.cpp:1343
void set_data(const GNdarray &array)
Set data from n-dimensional array.
Definition: GFft.cpp:603
GArf operator/(const GArf &arf, const double &scale)
Auxiliary Response File vision operator friend.
Definition: GArf.hpp:357
std::complex< double > & operator()(const int &ix)
1-dimensional FFT element access operator
Definition: GFft.hpp:207
void forward(const GNdarray &array)
Forward Fast Fourier Transform.
Definition: GFft.cpp:328
GFft * clone(void) const
Clone Fast Fourier Transform.
Definition: GFft.cpp:294
std::string print(const GChatter &chatter=NORMAL) const
Print Fast Fourier Transform information.
Definition: GFft.cpp:485
std::complex< double > timesi(const std::complex< double > &value) const
Return complex value times i.
Definition: GFft.hpp:397
std::vector< std::complex< double > > get_w(const GFftWavetable &wavetable, const int &index, const int &k, const int &q, const int &n, const int &sign) const
Extract coefficients from wavetable.
Definition: GFft.cpp:1535
void clear(void)
Clear Fast Fourier Transform.
Definition: GFft.cpp:276
std::complex< double > * m_data
Pointer on array data.
Definition: GFft.hpp:181
void factor5(const std::complex< double > *in, const int &istride, std::complex< double > *out, const int &ostride, const GFftWavetable &wavetable, const int &sign, const int &product, const int &n, const int &index)
Compute FFT for factor 5.
Definition: GFft.cpp:1029
Definition of interface for all GammaLib classes.
void init_members(void)
Initialise class members.
Definition: GFft.cpp:539
GFft & operator/=(const GFft &fft)
Unary division operator.
Definition: GFft.cpp:230
virtual ~GFft(void)
Destructor.
Definition: GFft.cpp:106
Lookup table class for Fast Fourier Transformation.
void transform(std::complex< double > *data, const int &stride, const int &n, const GFftWavetable &wavetable, const bool &forward=true)
Perform Fast Fourier Transform.
Definition: GFft.cpp:718
int size(void) const
Return number of elements in array.
Definition: GFft.hpp:282
GNdarray backward(void) const
Backward Fast Fourier Transform.
Definition: GFft.cpp:401
int m_size
Size of data array.
Definition: GFft.hpp:180
GFft & operator-=(const GFft &fft)
Unary subtraction operator.
Definition: GFft.cpp:182
GArf operator+(const GArf &a, const GArf &b)
Auxiliary Response File addition operator friend.
Definition: GArf.hpp:293
bool has_same_shape(const GFft &fft) const
Check if FFT has the same shape.
Definition: GFft.cpp:646
void free_members(void)
Delete class members.
Definition: GFft.cpp:582
GFft & operator=(const GFft &fft)
Assignment operator.
Definition: GFft.cpp:128
GFft & operator+=(const GFft &fft)
Unary addition operator.
Definition: GFft.cpp:158
void copy_members(const GFft &fft)
Copy class members.
Definition: GFft.cpp:558
GFft(void)
Void constructor.
Definition: GFft.cpp:55
std::vector< GFftWavetable > m_wavetable
Trigonometric coefficients.
Definition: GFft.hpp:184
Lookup table class interface definition for Fast Fourier transformation.
GFft & operator*=(const GFft &fft)
Unary multiplication operator.
Definition: GFft.cpp:206
Interface class for all GammaLib classes.
Definition: GBase.hpp:52
Fast Fourier Transformation class.
Definition: GFft.hpp:57
void factor7(const std::complex< double > *in, const int &istride, std::complex< double > *out, const int &ostride, const GFftWavetable &wavetable, const int &sign, const int &product, const int &n, const int &index)
Compute FFT for factor 7.
Definition: GFft.cpp:1214
std::vector< int > m_strides
Steps in each dimension.
Definition: GFft.hpp:183
GNdarray sign(const GNdarray &array)
Computes sign of array elements.
Definition: GNdarray.cpp:1269
GChatter
Definition: GTypemaps.hpp:33
void require_same_shape(const std::string &method, const GFft &fft) const
Throw exception if FFT shapes differ.
Definition: GFft.cpp:676
GArf operator*(const GArf &arf, const double &scale)
Auxiliary Response File scaling operator friend.
Definition: GArf.hpp:325
N-dimensional array class.
Definition: GNdarray.hpp:44
GFft operator-(void) const
Unary minus operator.
Definition: GFft.cpp:252
const std::vector< int > & shape(void) const
Return shape of array.
Definition: GFft.hpp:296
void factor2(const std::complex< double > *in, const int &istride, std::complex< double > *out, const int &ostride, const GFftWavetable &wavetable, const int &sign, const int &product, const int &n, const int &index)
Compute FFT for factor 2.
Definition: GFft.cpp:825
int dim(void) const
Return dimension of Fast Fourier Transformation.
Definition: GFft.hpp:268
void factor6(const std::complex< double > *in, const int &istride, std::complex< double > *out, const int &ostride, const GFftWavetable &wavetable, const int &sign, const int &product, const int &n, const int &index)
Compute FFT for factor 6.
Definition: GFft.cpp:1118
GArf operator-(const GArf &a, const GArf &b)
Auxiliary Response File subtraction operator friend.
Definition: GArf.hpp:309
std::string classname(void) const
Return class name.
Definition: GFft.hpp:194
const std::vector< int > & strides(void) const
Return strides of array.
Definition: GFft.hpp:310
void factor4(const std::complex< double > *in, const int &istride, std::complex< double > *out, const int &ostride, const GFftWavetable &wavetable, const int &sign, const int &product, const int &n, const int &index)
Compute FFT for factor 4.
Definition: GFft.cpp:957