GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GFftWavetable.cpp
Go to the documentation of this file.
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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  ***************************************************************************/
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.";
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  // Loop over all output coefficients
352  for (int j = 1; j < m_factors[i]; ++j) {
353 
354  // Loop over all coefficients
355  for (int k = 0, m = 0; k < q; ++k) {
356 
357  // Compute theta value
358  m = m + j * product_1;
359  m = m % n;
360  double theta = d_theta * m; /* d_theta*j*k*p_(i-1) */
361 
362  // Compute trigonometric lookup table value
363  std::complex<double> value(std::cos(theta), std::sin(theta));
364 
365  // Push back value on lookup table
366  m_trig.push_back(value);
367 
368  } // endfor: loop over all elements for a given factor
369 
370  } // endfor
371 
372  } // endfor: computed ...
373 
374  // Return
375  return;
376 }
377 
378 
379 /***********************************************************************//**
380  * @brief Compute FFT factorisation
381  *
382  * @param[in] n Length of array to be transformed.
383  *
384  * @exception GException::invalid_argument
385  * Non positive array length specified
386  * @exception GException::invalid_value
387  * Product of factorisation factors are not equal to array length
388  *
389  * Computes the factorisation
390  *
391  * \f[
392  * n = \prod_i p_i
393  * \f]
394  *
395  * of an array of length @p n.
396  *
397  * The method is a C++ implementation of the fft_complex_factorize() and the
398  * fft_factorize() functions of the GSL library (version 2.2.1), defined in
399  * the file fft/factorize.c.
400  ***************************************************************************/
401 void GFftWavetable::set_factors(const int& n)
402 {
403  // Make sure that length is positive
404  if (n < 1) {
405  std::string msg = "Invalid array length "+gammalib::str(n)+". Array "
406  "needs to be of positive length.";
408  }
409 
410  // Set FFT factorisations. Other factors can be added here if their
411  // transform modules are implemented.
412  std::vector<int> subtransforms;
413  subtransforms.push_back(7);
414  subtransforms.push_back(6);
415  subtransforms.push_back(5);
416  subtransforms.push_back(4);
417  subtransforms.push_back(3);
418  subtransforms.push_back(2);
419 
420  // If array has unit length then set the first and only factor to 1 ...
421  if (n == 1) {
422  m_factors.push_back(1);
423  }
424 
425  // ... otherwise
426  else {
427 
428  // Save array length
429  int ntest = n;
430 
431  // Deal with the implemented factors first
432  for (int i = 0; i < subtransforms.size() && ntest > 1; ++i) {
433  int factor = subtransforms[i];
434  while ((ntest % factor) == 0) {
435  ntest /= factor;
436  m_factors.push_back(factor);
437  }
438  }
439 
440  // Deal with any other even prime factors (there is only one)
441  int factor = 2;
442  while ((ntest % factor) == 0 && (ntest != 1)) {
443  ntest /= factor;
444  m_factors.push_back(factor);
445  }
446 
447  // Deal with any other odd prime factors
448  factor = 3;
449  while (ntest != 1) {
450  while ((ntest % factor) != 0) {
451  factor += 2;
452  }
453  ntest /= factor;
454  m_factors.push_back(factor);
455  }
456 
457  // Check that the factorization is correct
458  int product = 1;
459  for (int i = 0; i < m_factors.size(); ++i) {
460  product *= m_factors[i];
461  }
462  if (product != n) {
463  std::string msg = "Product "+gammalib::str(product)+" of "
464  "factorisation factors differs from array "
465  "length "+gammalib::str(n)+".";
467  }
468 
469  } // endelse: array length is larger than 1
470 
471  // Return
472  return;
473 }
void free_members(void)
Delete class members.
void clear(void)
Clear lookup table for Fast Fourier Transform.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
Lookup table class for Fast Fourier Transformation.
#define G_FACTOR
std::vector< int > m_twiddle
Start index of factors.
Gammalib tools definition.
GFftWavetable * clone(void) const
Clone lookup table for Fast Fourier Transform.
int factor(const int &index) const
Return factorisation factor.
void copy_members(const GFftWavetable &wavetable)
Copy class members.
std::vector< std::complex< double > > m_trig
Trigonometric coefficients.
Lookup table class interface definition for Fast Fourier transformation.
int factors(void) const
Return number of factorisation factors.
GChatter
Definition: GTypemaps.hpp:33
void init_members(void)
Initialise class members.
#define G_SET_MEMBERS
GFftWavetable & operator=(const GFftWavetable &wavetable)
Assignment operator.
#define G_SET_FACTORS
virtual ~GFftWavetable(void)
Destructor.
std::string print(const GChatter &chatter=NORMAL) const
Print lookup table for Fast Fourier Transform information.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Exception handler interface definition.
void set_factors(const int &n)
Compute FFT factorisation.
int size(void) const
Return number of trigonometric coefficients.
const double twopi
Definition: GMath.hpp:36
void set_members(const int &n)
Set wavetable.
GFftWavetable(void)
Void constructor.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
std::vector< int > m_factors
Wavetable factors.
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413