queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
RealFft.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/Fft.h>
26 #include <gsl/gsl_fft_real.h>
27 #include <gsl/gsl_fft_complex.h>
28 
29 namespace QUESO {
30 
31 // Math methods------------------------------------------
32 template <>
33 void
35  const std::vector<double>& data,
36  unsigned int fftSize,
37  std::vector<std::complex<double> >& forwardResult)
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 }
99 //-------------------------------------------------------
100 template <>
101 void
103  const std::vector<double>& data,
104  unsigned int fftSize,
105  std::vector<std::complex<double> >& inverseResult)
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 }
152 #if 0
153 template <class T>
154 void
155 Fft<T>::allocTables(unsigned int fftSize)
156 {
157  if (m_fftSize != fftSize) {
158  if (m_fftSize != 0) freeTables();
159  m_fftSize = fftSize;
160  m_realWkSpace = gsl_fft_real_workspace_alloc (fftSize);
161  m_realWvTable = gsl_fft_real_wavetable_alloc (fftSize);
162  //m_hcWvTable = gsl_fft_halfcomplex_wavetable_alloc(fftSize);
163  }
164 
165  return;
166 }
167 
168 template <class T>
169 void
170 Fft<T>::freeTables()
171 {
172  if (m_fftSize != 0) {
173  //gsl_fft_halfcomplex_wavetable_free(m_hcWvTable);
174  gsl_fft_real_wavetable_free (m_realWvTable);
175  gsl_fft_real_workspace_free (m_realWkSpace);
176  m_fftSize = 0;
177  }
178 
179  return;
180 }
181 #endif
182 
183 } // End namespace QUESO
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).
Class for a Fast Fourier Transform (FFT) algorithm.
Definition: Fft.h:66

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5