queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ComplexFft.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_complex.h>
27 
28 namespace QUESO {
29 
30 template <>
31 void
33  const std::vector<std::complex<double> >& data,
34  unsigned int fftSize,
35  std::vector<std::complex<double> >& forwardResult)
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 }
91 
92 template <>
93 void
95  const std::vector<std::complex<double> >& data,
96  unsigned int fftSize,
97  std::vector<std::complex<double> >& inverseResult)
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 }
129 
130 } // End namespace QUESO
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