queso-0.53.0
Public Member Functions | Private Attributes | List of all members
QUESO::Fft< T > Class Template Reference

Class for a Fast Fourier Transform (FFT) algorithm. More...

#include <Fft.h>

Collaboration diagram for QUESO::Fft< T >:
Collaboration graph
[legend]

Public Member Functions

template<>
void forward (const std::vector< std::complex< double > > &data, unsigned int fftSize, std::vector< std::complex< double > > &forwardResult)
 
template<>
void inverse (const std::vector< std::complex< double > > &data, unsigned int fftSize, std::vector< std::complex< double > > &inverseResult)
 
template<>
void forward (const std::vector< double > &data, unsigned int fftSize, std::vector< std::complex< double > > &forwardResult)
 
template<>
void inverse (const std::vector< double > &data, unsigned int fftSize, std::vector< std::complex< double > > &inverseResult)
 
Constructor/Destructor methods
 Fft (const BaseEnvironment &env)
 Default Constructor. More...
 
 ~Fft ()
 Destructor. More...
 
Mathematical methods
void forward (const std::vector< T > &data, unsigned int fftSize, std::vector< std::complex< double > > &result)
 Calculates the forward Fourier transform (for real data. TODO: complex data). More...
 
void inverse (const std::vector< T > &data, unsigned int fftSize, std::vector< std::complex< double > > &result)
 Calculates the inverse Fourier transform for real and complex data. More...
 

Private Attributes

const BaseEnvironment & m_env
 

Detailed Description

template<class T>
class QUESO::Fft< T >

Class for a Fast Fourier Transform (FFT) algorithm.

This class implements a Fast Fourier Transform (FFT) algorithm. Fast Fourier Transforms are efficient algorithms for calculating the discrete Fourier transform (DFT),

\[ x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N) \]

and its inverse. The DFT usually arises as an approximation to the continuous Fourier transform when functions are sampled at discrete intervals in space or time. The naive evaluation of the discrete Fourier transform is a matrix-vector multiplication $ Ab $. A general matrix-vector multiplication takes $ O(N^2)$ operations for $ N $ data-points. Fast Fourier transform algorithms use a divide-and-conquer strategy to factorize the matrix $ A $ into smaller sub-matrices, corresponding to the integer factors of the length $ N $. If $ N $ can be factorized into a product of integers $ f_1 f_2 ... f_n $ then the DFT can be computed in $ O(N \sum f_i) $ operations. For a radix-2 FFT this gives an operation count of $ O(N \log_2 N)$.

Todo:
: Implement Forward Fourier Transform for Complex data.

Definition at line 66 of file Fft.h.

Constructor & Destructor Documentation

template<class T >
QUESO::Fft< T >::Fft ( const BaseEnvironment &  env)

Default Constructor.

Definition at line 31 of file Fft.C.

32  :
33  m_env(env)
34 {
35 }
QUESO::Fft::m_env
const BaseEnvironment & m_env
Definition: Fft.h:113
template<class T >
QUESO::Fft< T >::~Fft ( )

Destructor.

Definition at line 38 of file Fft.C.

39 {
40 }

Member Function Documentation

template<>
void QUESO::Fft< std::complex< double > >::forward ( const std::vector< std::complex< double > > &  data,
unsigned int  fftSize,
std::vector< std::complex< double > > &  forwardResult 
)

Definition at line 32 of file ComplexFft.C.

References queso_not_implemented.

36 {
37  queso_not_implemented();
38 
39  std::complex<double> z = data[0]; z += 0.; // just to avoid icpc warnings
40  unsigned int f = fftSize; f += 1; // just to avoid icpc warnings
41  forwardResult[0] = 0.; // just to avoid icpc warnings
42 #if 0
43  if (forwardResult.size() != fftSize) {
44  forwardResult.resize(fftSize,std::complex<double>(0.,0.));
45  std::vector<std::complex<double> >(forwardResult).swap(forwardResult);
46  }
47 
48  std::vector<double> internalData(fftSize,0.);
49  unsigned int minSize = std::min((unsigned int) data.size(),fftSize);
50  for (unsigned int j = 0; j < minSize; ++j) {
51  internalData[j] = data[j];
52  }
53 
54  gsl_fft_real_workspace* realWkSpace = gsl_fft_real_workspace_alloc(fftSize);
55  gsl_fft_real_wavetable* realWvTable = gsl_fft_real_wavetable_alloc(fftSize);
56 
57  gsl_fft_real_transform(&internalData[0],
58  1,
59  fftSize,
60  realWvTable,
61  realWkSpace);
62 
63  gsl_fft_real_wavetable_free(realWvTable);
64  gsl_fft_real_workspace_free(realWkSpace);
65 
66  unsigned int halfFFTSize = fftSize/2;
67  double realPartOfFFT = 0.;
68  double imagPartOfFFT = 0.;
69  for (unsigned int j = 0; j < internalData.size(); ++j) {
70  if (j == 0) {
71  realPartOfFFT = internalData[j];
72  imagPartOfFFT = 0.;
73  }
74  else if (j < halfFFTSize) {
75  realPartOfFFT = internalData[2*j-1];
76  imagPartOfFFT = internalData[2*j ];
77  }
78  else if (j == halfFFTSize) {
79  realPartOfFFT = internalData[2*j-1];
80  imagPartOfFFT = 0.;
81  }
82  else {
83  realPartOfFFT = internalData[2*(fftSize-j)-1];
84  imagPartOfFFT = -internalData[2*(fftSize-j) ];
85  }
86  forwardResult[j] = std::complex<double>(realPartOfFFT,imagPartOfFFT);
87  }
88 #endif
89  return;
90 }
queso_not_implemented
#define queso_not_implemented()
Definition: asserts.h:56
template<>
void QUESO::Fft< double >::forward ( const std::vector< double > &  data,
unsigned int  fftSize,
std::vector< std::complex< double > > &  forwardResult 
)

Definition at line 34 of file RealFft.C.

38 {
39  if (forwardResult.size() != fftSize) {
40  forwardResult.resize(fftSize,std::complex<double>(0.,0.));
41  std::vector<std::complex<double> >(forwardResult).swap(forwardResult);
42  }
43 
44  std::vector<double> internalData(fftSize,0.);
45  unsigned int minSize = std::min((unsigned int) data.size(),fftSize);
46  for (unsigned int j = 0; j < minSize; ++j) {
47  internalData[j] = data[j];
48  }
49 
50  //double sumOfAllTerms = 0.;
51  //for (unsigned int j = 0; j < fftSize; ++j) {
52  // sumOfAllTerms += internalData[j];
53  //}
54 
55  //allocTables(fftSize);
56  gsl_fft_real_workspace* realWkSpace = gsl_fft_real_workspace_alloc(fftSize);
57  gsl_fft_real_wavetable* realWvTable = gsl_fft_real_wavetable_alloc(fftSize);
58 
59  gsl_fft_real_transform(&internalData[0],
60  1,
61  fftSize,
62  realWvTable,
63  realWkSpace);
64 
65  gsl_fft_real_wavetable_free(realWvTable);
66  gsl_fft_real_workspace_free(realWkSpace);
67  //freeTables();
68 
69  //std::cout << "After FFT"
70  // << ", sumOfAllTerms = " << sumOfAllTerms
71  // << ", sumOfAllTerms - dft[0] = " << sumOfAllTerms - internalData[0]
72  // << std::endl;
73 
74  unsigned int halfFFTSize = fftSize/2;
75  double realPartOfFFT = 0.;
76  double imagPartOfFFT = 0.;
77  for (unsigned int j = 0; j < internalData.size(); ++j) {
78  if (j == 0) {
79  realPartOfFFT = internalData[j];
80  imagPartOfFFT = 0.;
81  }
82  else if (j < halfFFTSize) {
83  realPartOfFFT = internalData[2*j-1];
84  imagPartOfFFT = internalData[2*j ];
85  }
86  else if (j == halfFFTSize) {
87  realPartOfFFT = internalData[2*j-1];
88  imagPartOfFFT = 0.;
89  }
90  else {
91  realPartOfFFT = internalData[2*(fftSize-j)-1];
92  imagPartOfFFT = -internalData[2*(fftSize-j) ];
93  }
94  forwardResult[j] = std::complex<double>(realPartOfFFT,imagPartOfFFT);
95  }
96 
97  return;
98 }
template<class T>
void QUESO::Fft< T >::forward ( const std::vector< T > &  data,
unsigned int  fftSize,
std::vector< std::complex< double > > &  result 
)

Calculates the forward Fourier transform (for real data. TODO: complex data).

This function uses GSL function 'gsl_fft_real_transform' to compute the FFT of data, a real array (of time-ordered real data) of length fftSize, using a mixed radix decimation-in-frequency algorithm. There is no restriction on the length fftSize. Efficient modules are provided for subtransforms of length 2, 3, 4 and 5. Any remaining factors are computed with a slow, $ O(\c fftSize ^2) $, general- fftSize module. The caller must supply a wavetable containing trigonometric lookup tables and a workspace work.
The definition of the forward Fourier transform, $ x = FFT(z) $ of size $ N $ is:

\[ x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N). \]

Referenced by QUESO::ScalarSequence< T >::autoCorrViaFft().

template<>
void QUESO::Fft< std::complex< double > >::inverse ( const std::vector< std::complex< double > > &  data,
unsigned int  fftSize,
std::vector< std::complex< double > > &  inverseResult 
)

Definition at line 94 of file ComplexFft.C.

98 {
99  if (inverseResult.size() != fftSize) {
100  inverseResult.resize(fftSize,std::complex<double>(0.,0.));
101  std::vector<std::complex<double> >(inverseResult).swap(inverseResult);
102  }
103 
104  std::vector<double> internalData(2*fftSize,0.); // Yes, twice the fftSize
105  unsigned int minSize = 2 * std::min((unsigned int) data.size(),fftSize); // Yes, 2*
106  for (unsigned int j = 0; j < minSize; ++j) {
107  internalData[2*j ] = data[j].real();
108  internalData[2*j+1] = data[j].imag();
109  }
110 
111  gsl_fft_complex_workspace* complexWkSpace = gsl_fft_complex_workspace_alloc(fftSize);
112  gsl_fft_complex_wavetable* complexWvTable = gsl_fft_complex_wavetable_alloc(fftSize);
113 
114  gsl_fft_complex_inverse(&internalData[0],
115  1,
116  fftSize,
117  complexWvTable,
118  complexWkSpace);
119 
120  gsl_fft_complex_wavetable_free(complexWvTable);
121  gsl_fft_complex_workspace_free(complexWkSpace);
122 
123  for (unsigned int j = 0; j < fftSize; ++j) {
124  inverseResult[j] = std::complex<double>(internalData[2*j],internalData[2*j+1]);
125  }
126 
127  return;
128 }
template<>
void QUESO::Fft< double >::inverse ( const std::vector< double > &  data,
unsigned int  fftSize,
std::vector< std::complex< double > > &  inverseResult 
)

Definition at line 102 of file RealFft.C.

106 {
107  if (inverseResult.size() != fftSize) {
108  inverseResult.resize(fftSize,std::complex<double>(0.,0.));
109  std::vector<std::complex<double> >(inverseResult).swap(inverseResult);
110  }
111 
112  std::vector<double> internalData(2*fftSize,0.); // Yes, twice the fftSize
113  unsigned int minSize = std::min((unsigned int) data.size(),fftSize);
114  for (unsigned int j = 0; j < minSize; ++j) {
115  internalData[2*j] = data[j];
116  }
117 
118  //if (m_subDisplayFile()) {
119  // *m_subDisplayFile() << "In Fft<double>::inverse()"
120  // << ": about to call gsl_fft_complex_inverse()"
121  // << " with fftSize = " << fftSize
122  // << "; internalData.size() = " << internalData.size()
123  // << std::endl;
124  //}
125 
126  gsl_fft_complex_workspace* complexWkSpace = gsl_fft_complex_workspace_alloc(fftSize);
127  gsl_fft_complex_wavetable* complexWvTable = gsl_fft_complex_wavetable_alloc(fftSize);
128 
129  gsl_fft_complex_inverse(&internalData[0],
130  1,
131  fftSize,
132  complexWvTable,
133  complexWkSpace);
134 
135  gsl_fft_complex_wavetable_free(complexWvTable);
136  gsl_fft_complex_workspace_free(complexWkSpace);
137 
138  //if (m_subDisplayFile()) {
139  // *m_subDisplayFile() << "In Fft<double>::inverse()"
140  // << ": returned from gsl_fft_complex_inverse()"
141  // << " with fftSize = " << fftSize
142  // << "; inverseResult.size() = " << inverseResult.size()
143  // << std::endl;
144  //}
145 
146  for (unsigned int j = 0; j < fftSize; ++j) {
147  inverseResult[j] = std::complex<double>(internalData[2*j],internalData[2*j+1]);
148  }
149 
150  return;
151 }
template<class T>
void QUESO::Fft< T >::inverse ( const std::vector< T > &  data,
unsigned int  fftSize,
std::vector< std::complex< double > > &  result 
)

Calculates the inverse Fourier transform for real and complex data.

This function uses GSL function 'gsl_fft_complex_inverse' to compute the FFT of data, a real array (of time-ordered real data) of length fftSize, using a mixed radix decimation-in-frequency algorithm. There is no restriction on the length fftSize. Efficient modules are provided for subtransforms of length 2, 3, 4 and 5. Any remaining factors are computed with a slow, $ O(fftSize^2) $, general- fftSize module. The caller must supply a wavetable containing trigonometric lookup tables and a workspace work.
The definition of the inverse Fourier transform, $ x = IFFT(z)$ of size $ N $ is:

\[ z_j = {1 \over N} \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N).\]

The factor of $ 1/N $ makes this a true inverse.

Referenced by QUESO::ScalarSequence< T >::autoCorrViaFft().

Member Data Documentation

template<class T>
const BaseEnvironment& QUESO::Fft< T >::m_env
private

Definition at line 113 of file Fft.h.

The documentation for this class was generated from the following files:
Generated on Thu Jun 11 2015 13:52:35 for queso-0.53.0 by  doxygen 1.8.5