GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GFft.cpp
Go to the documentation of this file.
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 #define G_OP_ADD "GFft::operator+=(GFft&)"
39 #define G_OP_SUB "GFft::operator-=(GFft&)"
40 #define G_OP_MUL "GFft::operator*=(GFft&)"
41 #define G_OP_DIV "GFft::operator/=(GFft&)"
42 #define G_FORWARD "GFft::forward(GNdarray&)"
43 #define G_BACKWARD "GFft::backward()"
44 
45 
46 /*==========================================================================
47  = =
48  = Constructors/destructors =
49  = =
50  ==========================================================================*/
51 
52 /***********************************************************************//**
53  * @brief Void constructor
54  ***************************************************************************/
56 {
57  // Initialise class members
58  init_members();
59 
60  // Return
61  return;
62 }
63 
64 
65 /***********************************************************************//**
66  * @brief N-dimensional array constructor
67  *
68  * @param[in] array N-dimensional array.
69  *
70  * Constructs a Fast Fourier Transformation from a n-dimensional array.
71  ***************************************************************************/
72 GFft::GFft(const GNdarray& array)
73 {
74  // Initialise class members
75  init_members();
76 
77  // Perform foward transformation
78  forward(array);
79 
80  // Return
81  return;
82 }
83 
84 
85 /***********************************************************************//**
86  * @brief Copy constructor
87  *
88  * @param[in] fft Fast Fourier Transform.
89  ***************************************************************************/
90 GFft::GFft(const GFft& fft)
91 {
92  // Initialise class members
93  init_members();
94 
95  // Copy members
96  copy_members(fft);
97 
98  // Return
99  return;
100 }
101 
102 
103 /***********************************************************************//**
104  * @brief Destructor
105  ***************************************************************************/
107 {
108  // Free members
109  free_members();
110 
111  // Return
112  return;
113 }
114 
115 
116 /*==========================================================================
117  = =
118  = Operators =
119  = =
120  ==========================================================================*/
121 
122 /***********************************************************************//**
123  * @brief Assignment operator
124  *
125  * @param[in] fft Fast Fourier Transform.
126  * @return Fast Fourier Transform.
127  ***************************************************************************/
129 {
130  // Execute only if object is not identical
131  if (this != &fft) {
132 
133  // Free members
134  free_members();
135 
136  // Initialise private members
137  init_members();
138 
139  // Copy members
140  copy_members(fft);
141 
142  } // endif: object was not identical
143 
144  // Return this object
145  return *this;
146 }
147 
148 
149 /***********************************************************************//**
150  * @brief Unary addition operator
151  *
152  * @param[in] fft Fast Fourier Transform.
153  * @return Fast Fourier Transform.
154  *
155  * Adds a Fast Fourier Transform to another Fast Fourier Transform. The
156  * method requires that both FFTs have the same dimension.
157  ***************************************************************************/
159 {
160  // Throw an exception if the FFTs have not the same shape
162 
163  // Add elements
164  for (int i = 0; i < m_size; ++i) {
165  *(m_data+i) += *(fft.m_data+i);
166  }
167 
168  // Return Fast Fourier Transform
169  return *this;
170 }
171 
172 
173 /***********************************************************************//**
174  * @brief Unary subtraction operator
175  *
176  * @param[in] fft Fast Fourier transform.
177  * @return Fast Fourier transform.
178  *
179  * Subtracts a Fast Fourier Transform from another Fast Fourier Transform.
180  * The method requires that both FFTs have the same dimension.
181  ***************************************************************************/
183 {
184  // Throw an exception if the FFTs have not the same shape
186 
187  // Subtract elements
188  for (int i = 0; i < m_size; ++i) {
189  *(m_data+i) -= *(fft.m_data+i);
190  }
191 
192  // Return Fast Fourier Transform
193  return *this;
194 }
195 
196 
197 /***********************************************************************//**
198  * @brief Unary multiplication operator
199  *
200  * @param[in] fft Fast Fourier transform.
201  * @return Fast Fourier transform.
202  *
203  * Multiplies a Fast Fourier Transform to another Fast Fourier Transform.
204  * The method requires that both FFTs have the same dimension.
205  ***************************************************************************/
207 {
208  // Throw an exception if the FFTs have not the same shape
210 
211  // Multiply elements
212  for (int i = 0; i < m_size; ++i) {
213  *(m_data+i) *= *(fft.m_data+i);
214  }
215 
216  // Return Fast Fourier Transform
217  return *this;
218 }
219 
220 
221 /***********************************************************************//**
222  * @brief Unary division operator
223  *
224  * @param[in] fft Fast Fourier transform.
225  * @return Fast Fourier transform.
226  *
227  * Divides a Fast Fourier Transform by another Fast Fourier Transform.
228  * The method requires that both FFTs have the same dimension.
229  ***************************************************************************/
231 {
232  // Throw an exception if the FFTs have not the same shape
234 
235  // Divide elements
236  for (int i = 0; i < m_size; ++i) {
237  *(m_data+i) /= *(fft.m_data+i);
238  }
239 
240  // Return Fast Fourier Transform
241  return *this;
242 }
243 
244 
245 /***********************************************************************//**
246  * @brief Unary minus operator
247  *
248  * @return Fast Fourier transformation.
249  *
250  * Negates all elements of the Fast Fourier Transformation.
251  ***************************************************************************/
253 {
254  // Copy FFT
255  GFft fft = *this;
256 
257  // Negate all elements
258  for (int i = 0; i < fft.size(); ++i) {
259  fft(i) = -fft(i);
260  }
261 
262  // Return FFT
263  return fft;
264 }
265 
266 
267 /*==========================================================================
268  = =
269  = Public methods =
270  = =
271  ==========================================================================*/
272 
273 /***********************************************************************//**
274  * @brief Clear Fast Fourier Transform
275  ***************************************************************************/
276 void GFft::clear(void)
277 {
278  // Free members
279  free_members();
280 
281  // Initialise private members
282  init_members();
283 
284  // Return
285  return;
286 }
287 
288 
289 /***********************************************************************//**
290  * @brief Clone Fast Fourier Transform
291  *
292  * @return Pointer to deep copy of Fast Fourier Transform.
293  ***************************************************************************/
294 GFft* GFft::clone(void) const
295 {
296  // Clone FFT
297  return new GFft(*this);
298 }
299 
300 
301 /***********************************************************************//**
302  * @brief Forward Fast Fourier Transform
303  *
304  * @param[in] array N-dimensional array.
305  *
306  * GException::feature_not_implemented
307  * Array dimension is not supported.
308  *
309  * Performs discret Fast Fourier forward transformation for 1- and
310  * 2-dimensional arrays.
311  *
312  * For one dimensional arrays of length \f$N\f$ the transformation is given
313  * by
314  *
315  * \f[
316  * X_k = \sum_{n=0}^{N-1} z_n e^{-i 2 \pi \frac{k n}{N}}
317  * \f]
318  *
319  * and for two dimensional arrays of length \f$ M \times N\f$ the
320  * transformation is given by
321  *
322  * \f[
323  * X_{k,l} = \sum_{n=0}^{N-1} \left(
324  * \sum_{m=0}^{M-1} z_{m,n} e^{-i 2 \pi \frac{m k}{M}} \right)
325  * e^{-i 2 \pi \frac{n l}{N}}
326  * \f]
327  ***************************************************************************/
328 void GFft::forward(const GNdarray& array)
329 {
330  // Set data
331  set_data(array);
332 
333  // Continue only if there are elements in the data member
334  if (m_size > 0) {
335 
336  // Case A: array is 1-dimensional
337  if (array.dim() == 1) {
338  transform(m_data, 1, m_shape[0], m_wavetable[0], true);
339  }
340 
341  // Case B: array is 2-dimensional
342  else if (array.dim() == 2) {
343 
344  // Loop over rows
345  for (int row = 0; row < m_shape[1]; ++row) {
346  transform(m_data + row*m_strides[1], m_strides[0], m_shape[0],
347  m_wavetable[0], true);
348  }
349 
350  // Loop over columns
351  for (int col = 0; col < m_shape[0]; ++col) {
352  transform(m_data + col*m_strides[0], m_strides[1], m_shape[1],
353  m_wavetable[1], true);
354  }
355 
356  }
357 
358  // ... otherwise throw an exception
359  else {
360  std::string msg = "Fast Fourier Transformation for "+
361  gammalib::str(array.dim())+"-dimensional array "
362  "not implemented.";
364  }
365 
366  } // endif: there were elements in the FFT data member
367 
368  // Return
369  return;
370 }
371 
372 
373 /***********************************************************************//**
374  * @brief Backward Fast Fourier Transform
375  *
376  * @return N-dimensional array.
377  *
378  * GException::feature_not_implemented
379  * Array dimension is not supported.
380  *
381  * Performs the inverse (or backward) discret Fast Fourier forward
382  * transformation for 1- and 2-dimensional arrays.
383  *
384  * For one dimensional arrays of length \f$N\f$ the transformation is given
385  * by
386  *
387  * \f[
388  * z_k = \frac{1}{N} \sum_{n=0}^{N-1} X_n e^{i 2 \pi \frac{k n}{N}}
389  * \f]
390  *
391  * and for two dimensional arrays of length \f$ M \times N\f$ the
392  * transformation is given by
393  *
394  * \f[
395  * z_{k,l} = \frac{1}{N} \sum_{n=0}^{N-1} \left(
396  * \frac{1}{M} \sum_{m=0}^{M-1} X_{m,n}
397  * e^{i 2 \pi \frac{m k}{M}} \right)
398  * e^{i 2 \pi \frac{n l}{N}}
399  * \f]
400  ***************************************************************************/
402 {
403  // Initialise n-dimensional array
404  GNdarray array(shape());
405 
406  // Continue only if there are elements in the array
407  if (m_size > 0) {
408 
409  // Create copy of complex array
410  std::complex<double>* data = new std::complex<double>[m_size];
411  for (int i = 0; i < m_size; ++i) {
412  *(data+i) = *(m_data+i);
413  }
414 
415  // Case A: array is 1-dimensional
416  if (array.dim() == 1) {
417 
418  // Perform backward transformation (circumvent const correctness)
419  const_cast<GFft*>(this)->transform(data, 1, m_shape[0], m_wavetable[0], false);
420 
421  // Get normalisation factor
422  const double norm = 1.0 / double(m_shape[0]);
423 
424  // Extract 1-dimensional array and normalise with 1/n
425  for (int i = 0; i < m_size; ++i) {
426  array(i) = (data+i)->real() * norm;
427  }
428 
429  }
430 
431  // Case B: array is 2-dimensional
432  else if (array.dim() == 2) {
433 
434  // Loop over rows
435  for (int row = 0; row < m_shape[1]; ++row) {
436  const_cast<GFft*>(this)->transform(data + row*m_strides[1], m_strides[0], m_shape[0],
437  m_wavetable[0], false);
438  }
439 
440  // Loop over columns
441  for (int col = 0; col < m_shape[0]; ++col) {
442  const_cast<GFft*>(this)->transform(data + col*m_strides[0], m_strides[1], m_shape[1],
443  m_wavetable[1], false);
444  }
445 
446  // Get normalisation factor
447  const double norm = 1.0 / double(m_shape[0]*m_shape[1]);
448 
449  // Extract 2-dimensional array and normalise with 1/(n*m)
450  for (int i = 0; i < m_size; ++i) {
451  array(i) = (data+i)->real() * norm;
452  }
453 
454  }
455 
456  // ... otherwise throw an exception
457  else {
458 
459  // Delete copy of complex array
460  delete [] data;
461 
462  // Throw exception
463  std::string msg = "Fast Fourier Transformation for "+
464  gammalib::str(array.dim())+"-dimensional array "
465  "not implemented.";
467  }
468 
469  // Delete copy of complex array
470  delete [] data;
471 
472  } // endif: there were elements in the array
473 
474  // Return N-dimensional array
475  return array;
476 }
477 
478 
479 /***********************************************************************//**
480  * @brief Print Fast Fourier Transform information
481  *
482  * @param[in] chatter Chattiness.
483  * @return String containing Fast Fourier Transform information.
484  ***************************************************************************/
485 std::string GFft::print(const GChatter& chatter) const
486 {
487  // Initialise result string
488  std::string result;
489 
490  // Continue only if chatter is not silent
491  if (chatter != SILENT) {
492 
493  // Append header
494  result.append("=== GFft ===");
495 
496  // Append array information
497  result.append("\n"+gammalib::parformat("Dimension"));
498  result.append(gammalib::str(dim()));
499  result.append("\n"+gammalib::parformat("Shape"));
500  result.append("(");
501  for (int i = 0; i < dim(); ++i) {
502  if (i > 0) {
503  result.append(",");
504  }
505  result.append(gammalib::str(m_shape[i]));
506  }
507  result.append(")");
508  result.append("\n"+gammalib::parformat("Size"));
509  result.append(gammalib::str(m_size));
510 
511  // VERBOSE: Put all FFT elements in string
512  if (chatter >= VERBOSE) {
513  result.append("\n"+gammalib::parformat("Elements"));
514  for (int i = 0; i < size(); ++i) {
515  if (i > 0) {
516  result.append(",");
517  }
518  result.append(gammalib::str((*this)(i)));
519  }
520  result.append(")");
521  }
522 
523  } // endif: chatter was not silent
524 
525  // Return result
526  return result;
527 }
528 
529 
530 /*==========================================================================
531  = =
532  = Private methods =
533  = =
534  ==========================================================================*/
535 
536 /***********************************************************************//**
537  * @brief Initialise class members
538  ***************************************************************************/
540 {
541  // Initialise members
542  m_size = 0;
543  m_data = NULL;
544  m_shape.clear();
545  m_strides.clear();
546  m_wavetable.clear();
547 
548  // Return
549  return;
550 }
551 
552 
553 /***********************************************************************//**
554  * @brief Copy class members
555  *
556  * @param[in] fft Fast Fourier Transform.
557  ***************************************************************************/
558 void GFft::copy_members(const GFft& fft)
559 {
560  // Copy members
561  m_size = fft.m_size;
562  m_shape = fft.m_shape;
563  m_strides = fft.m_strides;
564  m_wavetable = fft.m_wavetable;
565 
566  // Copy data
567  if (m_size > 0) {
568  m_data = new std::complex<double>[m_size];
569  for (int i = 0; i < m_size; ++i) {
570  *(m_data+i) = *(fft.m_data+i);
571  }
572  }
573 
574  // Return
575  return;
576 }
577 
578 
579 /***********************************************************************//**
580  * @brief Delete class members
581  ***************************************************************************/
583 {
584  // Delete data
585  if (m_data != NULL) {
586  delete [] m_data;
587  }
588 
589  // Reset size and pointer
590  m_size = 0;
591  m_data = NULL;
592 
593  // Return
594  return;
595 }
596 
597 
598 /***********************************************************************//**
599  * @brief Set data from n-dimensional array
600  *
601  * @param[in] array N-dimensional array
602  ***************************************************************************/
603 void GFft::set_data(const GNdarray& array)
604 {
605  // Clear information
606  free_members();
607  init_members();
608 
609  // Continue only if there are elements in the array
610  if (array.size() > 0) {
611 
612  // Set size of array
613  m_size = array.size();
614 
615  // Allocate data and initialise data to zero
616  m_data = new std::complex<double>[m_size];
617 
618  // Set real part of data from n-dimensional array
619  for (int i = 0; i < m_size; ++i) {
620  *(m_data+i) = std::complex<double>(array(i), 0.0);
621  }
622 
623  // Save shape and strides
624  m_shape = array.shape();
625  m_strides = array.strides();
626 
627  // Set trigonometric coefficients
628  for (int i = 0; i < m_shape.size(); ++i) {
629  GFftWavetable wavetable(m_shape[i]);
630  m_wavetable.push_back(wavetable);
631  }
632 
633  } // endif: there were elements in the array
634 
635  // Return
636  return;
637 }
638 
639 
640 /***********************************************************************//**
641  * @brief Check if FFT has the same shape
642  *
643  * @param[in] fft Fast Fourier Transform.
644  * @return True if the FFT has the same shape.
645  ***************************************************************************/
646 bool GFft::has_same_shape(const GFft& fft) const
647 {
648  // Check if the array has the same dimension
649  bool identity = m_shape.size() == fft.m_shape.size();
650 
651  // Check if the shape is identical. Break as soon as one shape value
652  // differs.
653  if (identity) {
654  for (int i = 0; i < m_shape.size(); ++i) {
655  if (m_shape[i] != fft.m_shape[i]) {
656  identity = false;
657  break;
658  }
659  }
660  }
661 
662  // Return result
663  return identity;
664 }
665 
666 
667 /***********************************************************************//**
668  * @brief Throw exception if FFT shapes differ
669  *
670  * @param[in] method Method that throws exception.
671  * @param[in] fft Fast Fourier Transform.
672  *
673  * @exception GException::invalid_argument
674  * FFT shapes differ.
675  ***************************************************************************/
676 void GFft::require_same_shape(const std::string& method,
677  const GFft& fft) const
678 {
679  // If the shape differs then throw an exception
680  if (!has_same_shape(fft)) {
681 
682  // Compose exception message
683  std::string msg = "Incompatible FFT dimensions (";
684  for (int i = 0; i < m_shape.size(); ++i) {
685  if (i > 0) {
686  msg += ", ";
687  }
688  msg += gammalib::str(m_shape[i]);
689  }
690  msg += ") <=> (";
691  for (int i = 0; i < fft.m_shape.size(); ++i) {
692  if (i > 0) {
693  msg += ", ";
694  }
695  msg += gammalib::str(fft.m_shape[i]);
696  }
697  msg += ").";
698 
699  // Throw exception
700  throw GException::invalid_argument(method, msg);
701 
702  } // endif: array shapes differed
703 
704  // Return
705  return;
706 }
707 
708 
709 /***********************************************************************//**
710  * @brief Perform Fast Fourier Transform
711  *
712  * @param[in] data Pointer to complex array to be transformed.
713  * @param[in] stride Step size when traversing complex array.
714  * @param[in] n Number of elements in complex array.
715  * @param[in] wavetable Trigonometric coefficients.
716  * @param[in] forward Forward transform if true, otherwise backward transform.
717  ***************************************************************************/
718 void GFft::transform(std::complex<double>* data,
719  const int& stride,
720  const int& n,
721  const GFftWavetable& wavetable,
722  const bool& forward)
723 {
724  // Allocate scratch space and initialize to zero
725  std::complex<double>* scratch = new std::complex<double>[n];
726  std::complex<double> zero(0.0, 0.0);
727  for (int i = 0; i < n; ++i) {
728  *(scratch+i) = zero;
729  }
730 
731  // Initialise transformation state
732  int state = 0;
733 
734  // Initialise factorisation product
735  int product = 1;
736 
737  // Set sign
738  int sign = (forward) ? -1 : +1;
739 
740  // Loop over all factors
741  for (int i = 0; i < wavetable.factors(); ++i) {
742 
743  // Get factorisation factor and start index
744  int factor = wavetable.factor(i);
745  int index = wavetable.index(i);
746 
747  // Update factorisation product
748  product *= factor;
749 
750  // Compute number of coefficients per factor
751  int q = n / product;
752 
753  // Set state dependent stuff
754  std::complex<double>* in;
755  std::complex<double>* out;
756  int istride;
757  int ostride;
758  if (state == 0) {
759  in = data;
760  istride = stride;
761  out = scratch;
762  ostride = 1;
763  state = 1;
764  }
765  else {
766  in = scratch;
767  istride = 1;
768  out = data;
769  ostride = stride;
770  state = 0;
771  }
772 
773  // Call factor dependent method
774  if (factor == 2) {
775  factor2(in, istride, out, ostride, wavetable, sign, product, n, index);
776  }
777  else if (factor == 3) {
778  factor3(in, istride, out, ostride, wavetable, sign, product, n, index);
779  }
780  else if (factor == 4) {
781  factor4(in, istride, out, ostride, wavetable, sign, product, n, index);
782  }
783  else if (factor == 5) {
784  factor5(in, istride, out, ostride, wavetable, sign, product, n, index);
785  }
786  else if (factor == 6) {
787  factor6(in, istride, out, ostride, wavetable, sign, product, n, index);
788  }
789  else if (factor == 7) {
790  factor7(in, istride, out, ostride, wavetable, sign, product, n, index);
791  }
792  else {
793  factorn(in, istride, out, ostride, wavetable, sign, factor, product, n, index);
794  }
795 
796  } // endfor: looped over all factors
797 
798  // In case the loop was exited in state 1 the results are in the scratch
799  // array and we need to copy the results back from the scratch array to
800  // the data array
801  if (state == 1) {
802  for (int i = 0; i < n; ++i) {
803  *(data+stride*i) = *(scratch+i);
804  }
805  }
806 
807  // Free scratch space
808  delete [] scratch;
809 
810  // Return
811  return;
812 }
813 
814 
815 /***********************************************************************//**
816  * @brief Compute FFT for factor 2
817  *
818  * @param[in] in Pointer to input array.
819  * @param[in] istride Step size when traversing input array.
820  * @param[in] out Pointer to output array.
821  * @param[in] ostride Step size when traversing output array.
822  * @param[in] wavetable Trigonometric coefficients.
823  * @param[in] sign Forward (-1) or backward (+1).
824  * @param[in] product ...
825  * @param[in] n Logical array length.
826  * @param[in] index Start index of trigonometric coefficients.
827  ***************************************************************************/
828 void GFft::factor2(const std::complex<double>* in,
829  const int& istride,
830  std::complex<double>* out,
831  const int& ostride,
832  const GFftWavetable& wavetable,
833  const int& sign,
834  const int& product,
835  const int& n,
836  const int& index)
837 {
838  // Compute ...
839  const int factor = 2;
840  const int m = n / factor;
841  const int q = n / product;
842  const int p_1 = product / factor;
843  const int jump = (factor - 1) * p_1;
844 
845  // Loop over ...
846  for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
847 
848  // Extract coefficients from wavetable
849  std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 1, sign);
850 
851  // Compute x = W(2) z
852  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
853 
854  // Get z
855  const std::complex<double> z0 = *(in+istride*i);
856  const std::complex<double> z1 = *(in+istride*(i+m));
857 
858  // Compute x
859  const std::complex<double> x0 = z0 + z1;
860  const std::complex<double> x1 = z0 - z1;
861 
862  // out = w * x
863  *(out+ostride*j) = x0;
864  *(out+ostride*(j+p_1)) = w[0] * x1;
865 
866  } // endfor: k1
867 
868  } // endfor: k
869 
870  // Return
871  return;
872 }
873 
874 
875 /***********************************************************************//**
876  * @brief Compute FFT for factor 3
877  *
878  * @param[in] in Pointer to input array.
879  * @param[in] istride Step size when traversing input array.
880  * @param[in] out Pointer to output array.
881  * @param[in] ostride Step size when traversing output array.
882  * @param[in] wavetable Trigonometric coefficients.
883  * @param[in] sign Forward (-1) or backward (+1).
884  * @param[in] product ...
885  * @param[in] n Logical array length.
886  * @param[in] index Start index of trigonometric coefficients.
887  ***************************************************************************/
888 void GFft::factor3(const std::complex<double>* in,
889  const int& istride,
890  std::complex<double>* out,
891  const int& ostride,
892  const GFftWavetable& wavetable,
893  const int& sign,
894  const int& product,
895  const int& n,
896  const int& index)
897 {
898  // Compute ...
899  const int factor = 3;
900  const int m = n / factor;
901  const int q = n / product;
902  const int p_1 = product / factor;
903  const int jump = (factor - 1) * p_1;
904 
905 
906  // Precompute some factors
907  const double tau = std::sqrt(3.0) / 2.0;
908 
909  // Loop over ...
910  for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
911 
912  // Extract coefficients from wavetable
913  std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 2, sign);
914 
915  // Compute x = W(3) z
916  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
917 
918  // Get z
919  const std::complex<double> z0 = *(in+istride*i);
920  const std::complex<double> z1 = *(in+istride*(i+m));
921  const std::complex<double> z2 = *(in+istride*(i+2*m));
922 
923  // Compute t
924  const std::complex<double> t1 = z1 + z2;
925  const std::complex<double> t2 = z0 - t1/2.0;
926  const std::complex<double> t3 = double((int)sign) * tau * (z1 - z2);
927 
928  // Compute x
929  const std::complex<double> x0 = z0 + t1;
930  const std::complex<double> x1 = t2 + timesi(t3);
931  const std::complex<double> x2 = t2 - timesi(t3);
932 
933  // out = w * x
934  *(out+ostride*j) = x0;
935  *(out+ostride*(j+p_1)) = w[0] * x1;
936  *(out+ostride*(j+2*p_1)) = w[1] * x2;
937 
938  } // endfor: k1
939 
940  } // endfor: k
941 
942  // Return
943  return;
944 }
945 
946 
947 /***********************************************************************//**
948  * @brief Compute FFT for factor 4
949  *
950  * @param[in] in Pointer to input array.
951  * @param[in] istride Step size when traversing input array.
952  * @param[in] out Pointer to output array.
953  * @param[in] ostride Step size when traversing output array.
954  * @param[in] wavetable Trigonometric coefficients.
955  * @param[in] sign Forward (-1) or backward (+1).
956  * @param[in] product ...
957  * @param[in] n Logical array length.
958  * @param[in] index Start index of trigonometric coefficients.
959  ***************************************************************************/
960 void GFft::factor4(const std::complex<double>* in,
961  const int& istride,
962  std::complex<double>* out,
963  const int& ostride,
964  const GFftWavetable& wavetable,
965  const int& sign,
966  const int& product,
967  const int& n,
968  const int& index)
969 {
970  // Compute ...
971  const int factor = 4;
972  const int m = n / factor;
973  const int q = n / product;
974  const int p_1 = product / factor;
975  const int jump = (factor - 1) * p_1;
976 
977  // Loop over ...
978  for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
979 
980  // Extract coefficients from wavetable
981  std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 3, sign);
982 
983  // Compute x = W(4) z
984  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
985 
986  // Get z
987  const std::complex<double> z0 = *(in+istride*i);
988  const std::complex<double> z1 = *(in+istride*(i+m));
989  const std::complex<double> z2 = *(in+istride*(i+2*m));
990  const std::complex<double> z3 = *(in+istride*(i+3*m));
991 
992  // Compute t
993  const std::complex<double> t1 = z0 + z2;
994  const std::complex<double> t2 = z1 + z3;
995  const std::complex<double> t3 = z0 - z2;
996  const std::complex<double> t4 = double((int)sign) * (z1 - z3);
997 
998  // Compute x
999  const std::complex<double> x0 = t1 + t2;
1000  const std::complex<double> x1 = t3 + timesi(t4);
1001  const std::complex<double> x2 = t1 - t2;
1002  const std::complex<double> x3 = t3 - timesi(t4);
1003 
1004  // out = w * x
1005  *(out+ostride*j) = x0;
1006  *(out+ostride*(j+p_1)) = w[0] * x1;
1007  *(out+ostride*(j+2*p_1)) = w[1] * x2;
1008  *(out+ostride*(j+3*p_1)) = w[2] * x3;
1009 
1010  } // endfor: k1
1011 
1012  } // endfor: k
1013 
1014  // Return
1015  return;
1016 }
1017 
1018 
1019 /***********************************************************************//**
1020  * @brief Compute FFT for factor 5
1021  *
1022  * @param[in] in Pointer to input array.
1023  * @param[in] istride Step size when traversing input array.
1024  * @param[in] out Pointer to output array.
1025  * @param[in] ostride Step size when traversing output array.
1026  * @param[in] wavetable Trigonometric coefficients.
1027  * @param[in] sign Forward (-1) or backward (+1).
1028  * @param[in] product ...
1029  * @param[in] n Logical array length.
1030  * @param[in] index Start index of trigonometric coefficients.
1031  ***************************************************************************/
1032 void GFft::factor5(const std::complex<double>* in,
1033  const int& istride,
1034  std::complex<double>* out,
1035  const int& ostride,
1036  const GFftWavetable& wavetable,
1037  const int& sign,
1038  const int& product,
1039  const int& n,
1040  const int& index)
1041 {
1042  // Compute ...
1043  const int factor = 5;
1044  const int m = n / factor;
1045  const int q = n / product;
1046  const int p_1 = product / factor;
1047  const int jump = (factor - 1) * p_1;
1048 
1049  // Precompute some factors
1050  const double sin_2pi_by_5 = std::sin(gammalib::twopi / 5.0);
1051  const double sin_2pi_by_10 = std::sin(gammalib::twopi / 10.0);
1052  const double tau = std::sqrt(5.0) / 4.0;
1053 
1054  // Loop over ...
1055  for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
1056 
1057  // Extract coefficients from wavetable
1058  std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 4, sign);
1059 
1060  // Compute x = W(5) z
1061  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1062 
1063  // Get z
1064  const std::complex<double> z0 = *(in+istride*i);
1065  const std::complex<double> z1 = *(in+istride*(i+m));
1066  const std::complex<double> z2 = *(in+istride*(i+2*m));
1067  const std::complex<double> z3 = *(in+istride*(i+3*m));
1068  const std::complex<double> z4 = *(in+istride*(i+4*m));
1069 
1070  // Compute t
1071  const std::complex<double> t1 = z1 + z4;
1072  const std::complex<double> t2 = z2 + z3;
1073  const std::complex<double> t3 = z1 - z4;
1074  const std::complex<double> t4 = z2 - z3;
1075  const std::complex<double> t5 = t1 + t2;
1076  const std::complex<double> t6 = tau * (t1 - t2);
1077  const std::complex<double> t7 = z0 - t5/4.0;
1078  const std::complex<double> t8 = t7 + t6;
1079  const std::complex<double> t9 = t7 - t6;
1080  const std::complex<double> t10 = double((int)sign) *
1081  (sin_2pi_by_5*t3 + sin_2pi_by_10*t4);
1082  const std::complex<double> t11 = double((int)sign) *
1083  (sin_2pi_by_10*t3 - sin_2pi_by_5*t4);
1084 
1085  // Compute x
1086  const std::complex<double> x0 = z0 + t5;
1087  const std::complex<double> x1 = t8 + timesi(t10);
1088  const std::complex<double> x2 = t9 + timesi(t11);
1089  const std::complex<double> x3 = t9 - timesi(t11);
1090  const std::complex<double> x4 = t8 - timesi(t10);
1091 
1092  // Compute out = w * x
1093  *(out+ostride*j) = x0;
1094  *(out+ostride*(j+p_1)) = w[0] * x1;
1095  *(out+ostride*(j+2*p_1)) = w[1] * x2;
1096  *(out+ostride*(j+3*p_1)) = w[2] * x3;
1097  *(out+ostride*(j+4*p_1)) = w[3] * x4;
1098 
1099  } // endfor: k1
1100 
1101  } // endfor: k
1102 
1103  // Return
1104  return;
1105 }
1106 
1107 
1108 /***********************************************************************//**
1109  * @brief Compute FFT for factor 6
1110  *
1111  * @param[in] in Pointer to input array.
1112  * @param[in] istride Step size when traversing input array.
1113  * @param[in] out Pointer to output array.
1114  * @param[in] ostride Step size when traversing output array.
1115  * @param[in] wavetable Trigonometric coefficients.
1116  * @param[in] sign Forward (-1) or backward (+1).
1117  * @param[in] product ...
1118  * @param[in] n Logical array length.
1119  * @param[in] index Start index of trigonometric coefficients.
1120  ***************************************************************************/
1121 void GFft::factor6(const std::complex<double>* in,
1122  const int& istride,
1123  std::complex<double>* out,
1124  const int& ostride,
1125  const GFftWavetable& wavetable,
1126  const int& sign,
1127  const int& product,
1128  const int& n,
1129  const int& index)
1130 {
1131  // Compute ...
1132  const int factor = 6;
1133  const int m = n / factor;
1134  const int q = n / product;
1135  const int p_1 = product / factor;
1136  const int jump = (factor - 1) * p_1;
1137 
1138  // Precompute some factors
1139  const double tau = std::sqrt(3.0) / 2.0;
1140 
1141  // Loop over ...
1142  for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
1143 
1144  // Extract coefficients from wavetable
1145  std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 5, sign);
1146 
1147  // Compute x = W(6) z. W(6) is a combination of sums and differences of
1148  // W(3) acting on the even and odd elements of z
1149  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1150 
1151  // Get z
1152  const std::complex<double> z0 = *(in+istride*i);
1153  const std::complex<double> z1 = *(in+istride*(i+m));
1154  const std::complex<double> z2 = *(in+istride*(i+2*m));
1155  const std::complex<double> z3 = *(in+istride*(i+3*m));
1156  const std::complex<double> z4 = *(in+istride*(i+4*m));
1157  const std::complex<double> z5 = *(in+istride*(i+5*m));
1158 
1159  // Compute ta
1160  const std::complex<double> ta1 = z2 + z4;
1161  const std::complex<double> ta2 = z0 - ta1/2.0;
1162  const std::complex<double> ta3 = double((int)sign) * tau * (z2 - z4);
1163 
1164  // Compute a
1165  const std::complex<double> a0 = z0 + ta1;
1166  const std::complex<double> a1 = ta2 + timesi(ta3);
1167  const std::complex<double> a2 = ta2 - timesi(ta3);
1168 
1169  // Compute tb
1170  const std::complex<double> tb1 = z5 + z1;
1171  const std::complex<double> tb2 = z3 - tb1/2.0;
1172  const std::complex<double> tb3 = double((int)sign) * tau * (z5 - z1);
1173 
1174  // Compute b
1175  const std::complex<double> b0 = z3 + tb1;
1176  const std::complex<double> b1 = tb2 + timesi(tb3);
1177  const std::complex<double> b2 = tb2 - timesi(tb3);
1178 
1179  // Compute x
1180  const std::complex<double> x0 = a0 + b0;
1181  const std::complex<double> x1 = a1 - b1;
1182  const std::complex<double> x2 = a2 + b2;
1183  const std::complex<double> x3 = a0 - b0;
1184  const std::complex<double> x4 = a1 + b1;
1185  const std::complex<double> x5 = a2 - b2;
1186 
1187  // Compute out = w * x
1188  *(out+ostride*j) = x0;
1189  *(out+ostride*(j+p_1)) = w[0] * x1;
1190  *(out+ostride*(j+2*p_1)) = w[1] * x2;
1191  *(out+ostride*(j+3*p_1)) = w[2] * x3;
1192  *(out+ostride*(j+4*p_1)) = w[3] * x4;
1193  *(out+ostride*(j+5*p_1)) = w[4] * x5;
1194 
1195  } // endfor: k1
1196 
1197  } // endfor: k
1198 
1199  // Return
1200  return;
1201 }
1202 
1203 
1204 /***********************************************************************//**
1205  * @brief Compute FFT for factor 7
1206  *
1207  * @param[in] in Pointer to input array.
1208  * @param[in] istride Step size when traversing input array.
1209  * @param[in] out Pointer to output array.
1210  * @param[in] ostride Step size when traversing output array.
1211  * @param[in] wavetable Trigonometric coefficients.
1212  * @param[in] sign Forward (-1) or backward (+1).
1213  * @param[in] product ...
1214  * @param[in] n Logical array length.
1215  * @param[in] index Start index of trigonometric coefficients.
1216  ***************************************************************************/
1217 void GFft::factor7(const std::complex<double>* in,
1218  const int& istride,
1219  std::complex<double>* out,
1220  const int& ostride,
1221  const GFftWavetable& wavetable,
1222  const int& sign,
1223  const int& product,
1224  const int& n,
1225  const int& index)
1226 {
1227  // Compute ...
1228  const int factor = 7;
1229  const int m = n / factor;
1230  const int q = n / product;
1231  const int p_1 = product / factor;
1232  const int jump = (factor - 1) * p_1;
1233 
1234  // Precompute some factors
1235  static const double twopi7 = gammalib::twopi / 7.0;
1236  static const double c1 = std::cos(1.0 * twopi7);
1237  static const double c2 = std::cos(2.0 * twopi7);
1238  static const double c3 = std::cos(3.0 * twopi7);
1239  static const double s1 = std::sin(1.0 * twopi7);
1240  static const double s2 = std::sin(2.0 * twopi7);
1241  static const double s3 = std::sin(3.0 * twopi7);
1242  static const double tau1 = (c1 + c2 + c3) / 3.0 - 1.0;
1243  static const double tau2 = (2.0 * c1 - c2 - c3) / 3.0;
1244  static const double tau3 = (c1 - 2.0*c2 + c3) / 3.0;
1245  static const double tau4 = (c1 + c2 - 2.0 * c3) / 3.0;
1246  static const double tau5 = (s1 + s2 - s3) / 3.0;
1247  static const double tau6 = (2.0 * s1 - s2 + s3) / 3.0;
1248  static const double tau7 = (s1 - 2.0 * s2 - s3) / 3.0;
1249  static const double tau8 = (s1 + s2 + 2.0 * s3) / 3.0;
1250 
1251  // Loop over ...
1252  for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
1253 
1254  // Extract coefficients from wavetable
1255  std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 6, sign);
1256 
1257  // Compute x = W(7) z
1258  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1259 
1260  // Get z
1261  const std::complex<double> z0 = *(in+istride*i);
1262  const std::complex<double> z1 = *(in+istride*(i+m));
1263  const std::complex<double> z2 = *(in+istride*(i+2*m));
1264  const std::complex<double> z3 = *(in+istride*(i+3*m));
1265  const std::complex<double> z4 = *(in+istride*(i+4*m));
1266  const std::complex<double> z5 = *(in+istride*(i+5*m));
1267  const std::complex<double> z6 = *(in+istride*(i+6*m));
1268 
1269  // Compute t
1270  const std::complex<double> t0 = z1 + z6;
1271  const std::complex<double> t1 = z1 - z6;
1272  const std::complex<double> t2 = z2 + z5;
1273  const std::complex<double> t3 = z2 - z5;
1274  const std::complex<double> t4 = z4 + z3;
1275  const std::complex<double> t5 = z4 - z3;
1276  const std::complex<double> t6 = t2 + t0;
1277  const std::complex<double> t7 = t5 + t3;
1278 
1279  // Compute b
1280  const std::complex<double> b0 = z0 + t6 + t4;
1281  const std::complex<double> b1 = tau1 * (t6 + t4);
1282  const std::complex<double> b2 = tau2 * (t0 - t4);
1283  const std::complex<double> b3 = tau3 * (t4 - t2);
1284  const std::complex<double> b4 = tau4 * (t2 - t0);
1285  const std::complex<double> b5 = double(-(int)sign) * tau5 * (t7 + t1);
1286  const std::complex<double> b6 = double(-(int)sign) * tau6 * (t1 - t5);
1287  const std::complex<double> b7 = double(-(int)sign) * tau7 * (t5 - t3);
1288  const std::complex<double> b8 = double(-(int)sign) * tau8 * (t3 - t1);
1289 
1290  // Compute T
1291  const std::complex<double> T0 = b0 + b1;
1292  const std::complex<double> T1 = b2 + b3;
1293  const std::complex<double> T2 = b4 - b3;
1294  const std::complex<double> T3 = -b2 - b4;
1295  const std::complex<double> T4 = b6 + b7;
1296  const std::complex<double> T5 = b8 - b7;
1297  const std::complex<double> T6 = -b8 - b6;
1298  const std::complex<double> T7 = T0 + T1;
1299  const std::complex<double> T8 = T0 + T2;
1300  const std::complex<double> T9 = T0 + T3;
1301  const std::complex<double> T10 = T4 + b5;
1302  const std::complex<double> T11 = T5 + b5;
1303  const std::complex<double> T12 = T6 + b5;
1304 
1305  // Compute x
1306  const std::complex<double> x0 = b0;
1307  const std::complex<double> x1 = T7 - timesi(T10);
1308  const std::complex<double> x2 = T9 - timesi(T12);
1309  const std::complex<double> x3 = T8 + timesi(T11);
1310  const std::complex<double> x4 = T8 - timesi(T11);
1311  const std::complex<double> x5 = T9 + timesi(T12);
1312  const std::complex<double> x6 = T7 + timesi(T10);
1313 
1314  // Compute out = w * x
1315  *(out+ostride*j) = x0;
1316  *(out+ostride*(j+p_1)) = w[0] * x1;
1317  *(out+ostride*(j+2*p_1)) = w[1] * x2;
1318  *(out+ostride*(j+3*p_1)) = w[2] * x3;
1319  *(out+ostride*(j+4*p_1)) = w[3] * x4;
1320  *(out+ostride*(j+5*p_1)) = w[4] * x5;
1321  *(out+ostride*(j+6*p_1)) = w[5] * x6;
1322 
1323  } // endfor: k1
1324 
1325  } // endfor: k
1326 
1327  // Return
1328  return;
1329 }
1330 
1331 
1332 /***********************************************************************//**
1333  * @brief Compute FFT for arbitrary factor
1334  *
1335  * @param[in] in Pointer to input array.
1336  * @param[in] istride Step size when traversing input array.
1337  * @param[in] out Pointer to output array.
1338  * @param[in] ostride Step size when traversing output array.
1339  * @param[in] wavetable Trigonometric coefficients.
1340  * @param[in] sign Forward (-1) or backward (+1).
1341  * @param[in] factor Factor.
1342  * @param[in] product ...
1343  * @param[in] n Logical array length.
1344  * @param[in] index Start index of trigonometric coefficients.
1345  ***************************************************************************/
1346 void GFft::factorn(std::complex<double>* in,
1347  const int& istride,
1348  std::complex<double>* out,
1349  const int& ostride,
1350  const GFftWavetable& wavetable,
1351  const int& sign,
1352  const int& factor,
1353  const int& product,
1354  const int& n,
1355  const int& index)
1356 {
1357  // Compute ...
1358  const int m = n / factor;
1359  const int q = n / product;
1360  const int p_1 = product / factor;
1361  const int jump = (factor - 1) * p_1;
1362 
1363  // ...
1364  for (int i = 0; i < m; ++i) {
1365  *(out+ostride*i) = *(in+istride*i);
1366  }
1367 
1368  // ...
1369  for (int e = 1; e < (factor - 1) / 2 + 1; ++e) {
1370  for (int i = 0; i < m; ++i) {
1371  const int idx = i + e * m;
1372  const int idxc = i + (factor - e) * m;
1373  *(out+ostride*idx) = *(in+istride*idx) + *(in+istride*idxc);
1374  *(out+ostride*idxc) = *(in+istride*idx) - *(in+istride*idxc);
1375  }
1376  }
1377 
1378  // e = 0
1379  for (int i = 0; i < m; ++i) {
1380  *(in+istride*i) = *(out+ostride*i);
1381  }
1382  for (int e = 1; e < (factor - 1) / 2 + 1; ++e) {
1383  for (int i = 0; i < m; ++i) {
1384  const int idx = i + e * m;
1385  *(in+istride*i) += *(out+ostride*idx);
1386  }
1387  }
1388 
1389  // ...
1390  for (int e = 1; e < (factor-1)/2 + 1; ++e) {
1391 
1392  //
1393  int idx = e * q;
1394  const int idx_step = e * q;
1395 
1396  //
1397  double w_real, w_imag;
1398 
1399  //
1400  const int em = e * m;
1401  const int ecm = (factor - e) * m;
1402 
1403  // ...
1404  for (int i = 0; i < m; ++i) {
1405  *(in+istride*(i+em)) = *(out+ostride*i);
1406  *(in+istride*(i+ecm)) = *(out+ostride*i);
1407  }
1408 
1409  // ...
1410  for (int e1 = 1; e1 < (factor - 1) / 2 + 1; ++e1) {
1411 
1412  //
1413  if (idx == 0) {
1414  w_real = 1.0;
1415  w_imag = 0.0;
1416  }
1417  else {
1418 
1419  // Compute indices
1420  int twiddle = index + idx - 1;
1421 
1422  // Set trigonometric coefficients for forward transform
1423  if (sign == -1) {
1424  w_real = wavetable[twiddle].real();
1425  w_imag = wavetable[twiddle].imag();
1426  }
1427  // ... otherwise set trigonometric coefficients for backward
1428  // tranform: w -> conjugate(w)
1429  else {
1430  w_real = wavetable[twiddle].real();
1431  w_imag = -wavetable[twiddle].imag();
1432  }
1433  }
1434 
1435  // Loop over ...
1436  for (int i = 0; i < m; ++i) {
1437 
1438  //
1439  const double xp_real = (out+ostride*(i + e1 * m))->real();
1440  const double xp_imag = (out+ostride*(i + e1 * m))->imag();
1441  const double xm_real = (out+ostride*(i + (factor-e1)*m))->real();
1442  const double xm_imag = (out+ostride*(i + (factor-e1)*m))->imag();
1443 
1444  //
1445  const double ap = w_real * xp_real;
1446  const double am = w_imag * xm_imag;
1447 
1448  //
1449  double sum_real = ap - am;
1450  double sumc_real = ap + am;
1451 
1452  //
1453  const double bp = w_real * xp_imag;
1454  const double bm = w_imag * xm_real;
1455 
1456  //
1457  double sum_imag = bp + bm;
1458  double sumc_imag = bp - bm;
1459 
1460  //
1461  *(in+istride*(i+em)) = std::complex<double>((in+istride*(i+em))->real() + sum_real,
1462  (in+istride*(i+em))->imag() + sum_imag);
1463  *(in+istride*(i+ecm)) = std::complex<double>((in+istride*(i+ecm))->real() + sumc_real,
1464  (in+istride*(i+ecm))->imag() + sumc_imag);
1465 
1466 
1467  } // endfor: i
1468 
1469  // Increment
1470  idx += idx_step ;
1471  idx %= factor * q ;
1472 
1473  } // endfor: e1
1474 
1475  } // endfor: e
1476 
1477  // k = 0
1478  for (int k1 = 0; k1 < p_1; ++k1) {
1479  *(out+ostride*k1) = *(in+istride*k1);
1480  }
1481 
1482  // k > 0
1483  for (int e1 = 1; e1 < factor; ++e1) {
1484  for (int k1 = 0; k1 < p_1; ++k1) {
1485  *(out+ostride*(k1 + e1 * p_1)) = *(in+istride*(k1 + e1 * m));
1486  }
1487  }
1488 
1489  // e = 0
1490  for (int k = 1, i = p_1, j = product; k < q; ++k, j += jump) {
1491  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1492  *(out+ostride*j) = *(in+istride*i);
1493  }
1494  }
1495 
1496  // e > 0
1497  for (int k = 1, i = p_1, j = product; k < q; ++k, j += jump) {
1498 
1499  for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1500 
1501  for (int e1 = 1; e1 < factor; ++e1) {
1502 
1503  // Get x
1504  std::complex<double> x = *(in+istride*(i + e1 * m));
1505 
1506  // Compute index
1507  int twiddle = index + (e1-1)*q + k-1;
1508 
1509  // Get w
1510  std::complex<double> w = (sign == -1)
1511  ? wavetable[twiddle]
1512  : std::conj(wavetable[twiddle]);
1513 
1514  // Compute out = w * x
1515  *(out+ostride*(j + e1 * p_1)) = w * x;
1516 
1517  } // endfor: e1
1518 
1519  } // endfor: k1
1520 
1521  } // endfor: k
1522 
1523  // Return
1524  return;
1525 }
1526 
1527 
1528 /***********************************************************************//**
1529  * @brief Extract coefficients from wavetable
1530  *
1531  * @param[in] index Start index of trigonometric coefficients.
1532  * @param[in] k ...
1533  * @param[in] q ...
1534  * @param[in] n Number of coefficients.
1535  * @param[in] sign Forward (-1) or backward (+1) transformation.
1536  * @param[in] wavetable Trigonometric coefficients.
1537  ***************************************************************************/
1538 std::vector<std::complex<double> > GFft::get_w(const GFftWavetable& wavetable,
1539  const int& index,
1540  const int& k,
1541  const int& q,
1542  const int& n,
1543  const int& sign) const
1544 {
1545  // Allocate w vectors
1546  std::vector<std::complex<double> > w(n);
1547 
1548  // Set trigonometric coefficients for k=0 since they are not stored in the
1549  // wavetable object
1550  if (k == 0) {
1551  w.assign(n, std::complex<double>(1.0, 0.0));
1552  }
1553 
1554  // ... otherwise use coefficients stored in wavetable
1555  else {
1556 
1557  // Compute indices
1558  int twiddle = index + k - 1;
1559 
1560  // If sign==-1 then set trigonometric coefficients for forward transform
1561  if (sign == -1) {
1562  for (int i = 0; i < n; ++i, twiddle += q) {
1563  w[i] = wavetable[twiddle];
1564  }
1565  }
1566 
1567  // ... otherwise set trigonometric coefficients for backward tranform
1568  // w -> conjugate(w)
1569  else {
1570  for (int i = 0; i < n; ++i, twiddle += q) {
1571  w[i] = std::conj(wavetable[twiddle]);
1572  }
1573  }
1574 
1575  } // endelse: k != 0
1576 
1577  // Return w vectors
1578  return w;
1579 }
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:888
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:1346
void set_data(const GNdarray &array)
Set data from n-dimensional array.
Definition: GFft.cpp:603
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
const std::vector< int > & strides(void) const
Return strides of array.
Definition: GNdarray.hpp:321
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:1538
void clear(void)
Clear Fast Fourier Transform.
Definition: GFft.cpp:276
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
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:1032
#define G_BACKWARD
Definition: GFft.cpp:43
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
Gammalib tools definition.
const std::vector< int > & shape(void) const
Return shape of array.
Definition: GNdarray.hpp:307
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
int factor(const int &index) const
Return factorisation factor.
GFft & operator-=(const GFft &fft)
Unary subtraction operator.
Definition: GFft.cpp:182
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
#define G_OP_DIV
Definition: GFft.cpp:41
std::vector< GFftWavetable > m_wavetable
Trigonometric coefficients.
Definition: GFft.hpp:184
#define G_OP_SUB
Definition: GFft.cpp:39
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
int factors(void) const
Return number of factorisation factors.
#define G_OP_ADD
Definition: GFft.cpp:38
GFft & operator*=(const GFft &fft)
Unary multiplication operator.
Definition: GFft.cpp:206
N-dimensional array class interface definition.
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:1217
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:1238
int size(void) const
Return number of elements in array.
Definition: GNdarray.hpp:293
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
const int & index(const int &factor) const
Return start index for a given factor.
N-dimensional array class.
Definition: GNdarray.hpp:44
GFft operator-(void) const
Unary minus operator.
Definition: GFft.cpp:252
Fast Fourier transformation class interface definition.
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:828
#define G_OP_MUL
Definition: GFft.cpp:40
#define G_FORWARD
Definition: GFft.cpp:42
int dim(void) const
Return dimension of Fast Fourier Transformation.
Definition: GFft.hpp:268
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Exception handler interface definition.
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:1121
const double twopi
Definition: GMath.hpp:36
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
int dim(void) const
Return dimension of array.
Definition: GNdarray.hpp:279
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:960
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413