25 #include <queso/Fft.h>
26 #include <gsl/gsl_fft_complex.h>
33 const std::vector<std::complex<double> >& data,
35 std::vector<std::complex<double> >& forwardResult)
39 "Fft<complex>::forward()",
40 "not implemented yet");
42 std::complex<double> z = data[0]; z += 0.;
43 unsigned int f = fftSize; f += 1;
44 forwardResult[0] = 0.;
46 if (forwardResult.size() != fftSize) {
47 forwardResult.resize(fftSize,std::complex<double>(0.,0.));
48 std::vector<std::complex<double> >(forwardResult).swap(forwardResult);
51 std::vector<double> internalData(fftSize,0.);
52 unsigned int minSize = std::min((
unsigned int) data.size(),fftSize);
53 for (
unsigned int j = 0; j < minSize; ++j) {
54 internalData[j] = data[j];
57 gsl_fft_real_workspace* realWkSpace = gsl_fft_real_workspace_alloc(fftSize);
58 gsl_fft_real_wavetable* realWvTable = gsl_fft_real_wavetable_alloc(fftSize);
60 gsl_fft_real_transform(&internalData[0],
66 gsl_fft_real_wavetable_free(realWvTable);
67 gsl_fft_real_workspace_free(realWkSpace);
69 unsigned int halfFFTSize = fftSize/2;
70 double realPartOfFFT = 0.;
71 double imagPartOfFFT = 0.;
72 for (
unsigned int j = 0; j < internalData.size(); ++j) {
74 realPartOfFFT = internalData[j];
77 else if (j < halfFFTSize) {
78 realPartOfFFT = internalData[2*j-1];
79 imagPartOfFFT = internalData[2*j ];
81 else if (j == halfFFTSize) {
82 realPartOfFFT = internalData[2*j-1];
86 realPartOfFFT = internalData[2*(fftSize-j)-1];
87 imagPartOfFFT = -internalData[2*(fftSize-j) ];
89 forwardResult[j] = std::complex<double>(realPartOfFFT,imagPartOfFFT);
98 const std::vector<std::complex<double> >& data,
100 std::vector<std::complex<double> >& inverseResult)
102 if (inverseResult.size() != fftSize) {
103 inverseResult.resize(fftSize,std::complex<double>(0.,0.));
104 std::vector<std::complex<double> >(inverseResult).swap(inverseResult);
107 std::vector<double> internalData(2*fftSize,0.);
108 unsigned int minSize = 2 * std::min((
unsigned int) data.size(),fftSize);
109 for (
unsigned int j = 0; j < minSize; ++j) {
110 internalData[2*j ] = data[j].real();
111 internalData[2*j+1] = data[j].imag();
114 gsl_fft_complex_workspace* complexWkSpace = gsl_fft_complex_workspace_alloc(fftSize);
115 gsl_fft_complex_wavetable* complexWvTable = gsl_fft_complex_wavetable_alloc(fftSize);
117 gsl_fft_complex_inverse(&internalData[0],
123 gsl_fft_complex_wavetable_free(complexWvTable);
124 gsl_fft_complex_workspace_free(complexWkSpace);
126 for (
unsigned int j = 0; j < fftSize; ++j) {
127 inverseResult[j] = std::complex<double>(internalData[2*j],internalData[2*j+1]);
const int UQ_UNAVAILABLE_RANK
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Class for a Fast Fourier Transform (FFT) algorithm.