1 2 #include "pffft.hpp" 3 4 #include <complex> 5 #include <iostream> 6 7 cxx11_forward_real_double(const int transformLen)8void cxx11_forward_real_double(const int transformLen) 9 { 10 std::cout << "running " << __FUNCTION__ << "()" << std::endl; 11 12 // first check - might be skipped 13 using FFT_T = pffft::Fft<double>; 14 if (transformLen < FFT_T::minFFtsize()) 15 { 16 std::cerr << "Error: minimum FFT transformation length is " << FFT_T::minFFtsize() << std::endl; 17 return; 18 } 19 20 // instantiate FFT and prepare transformation for length N 21 pffft::Fft<double> fft { transformLen }; 22 23 // one more check 24 if (!fft.isValid()) 25 { 26 std::cerr << "Error: transformation length " << transformLen << " is not decomposable into small prime factors. " 27 << "Next valid transform size is: " << FFT_T::nearestTransformSize(transformLen) 28 << "; next power of 2 is: " << FFT_T::nextPowerOfTwo(transformLen) << std::endl; 29 return; 30 } 31 32 // allocate aligned vectors for (real) input X and (complex) output Y 33 auto X = fft.valueVector(); // input vector; type is AlignedVector<double> 34 auto Y = fft.spectrumVector(); // output vector; type is AlignedVector< std::complex<double> > 35 36 // alternative access: get raw pointers to aligned vectors 37 double *Xs = X.data(); 38 std::complex<double> *Ys = Y.data(); 39 40 // prepare some input data 41 for (int k = 0; k < transformLen; k += 2) 42 { 43 X[k] = k; // access through AlignedVector<double> 44 Xs[k+1] = -1-k; // access through raw pointer 45 } 46 47 // do the forward transform; write complex spectrum result into Y 48 fft.forward(X, Y); 49 50 // print spectral output 51 std::cout << "output should be complex spectrum with " << fft.getSpectrumSize() << " bins" << std::endl; 52 std::cout << "output vector has size " << Y.size() << " (complex bins):" << std::endl; 53 for (unsigned k = 0; k < Y.size(); k += 2) 54 { 55 std::cout << "Y[" << k << "] = " << Y[k] << std::endl; 56 std::cout << "Y[" << k+1 << "] = " << Ys[k+1] << std::endl; 57 } 58 } 59 60 main(int argc,char * argv[])61int main(int argc, char *argv[]) 62 { 63 int N = (1 < argc) ? atoi(argv[1]) : 32; 64 cxx11_forward_real_double(N); 65 return 0; 66 } 67