Calculation of FFT and IFFT with the FFTW C ++ library

Question

I am trying to calculate FFT and then IFFT just try if I can get the same signal, but I'm not quite sure how to do this. Here's how I do FFT:

plan = fftw_plan_r2r_1d(blockSize, datas, out, FFTW_R2HC, FFTW_ESTIMATE); fftw_execute(plan); 

Answer 1

Have you at least tried to read more decent documentation?

They have a whole tutorial in which you can learn FFTW:

http://fftw.org/fftw3_doc/Tutorial.html#Tutorial

UPDATE: I assume that you know how to work with C-arrays, like what is used as input and output.

This page contains the information needed for FFT and IFFT (see Arguments-> Sign). The docs also say that input-> FFT-> IFFT-> n * input. Therefore, you need to be careful to properly scale the data.

Answer 2

Here is an example. It does two things. First, it prepares the input array in[N] as a cosine wave, whose frequency is 3, and the magnitude is 1.0, and Fourier transforms it. So, at the output you should see a peak in out[3] , and the other in out[N-3] . Since the magnitude of the cosine wave is 1.0, you get N / 2 with out[3] and out[N-3] .

Secondly, it reverses the Fourier transform of the array out[N] back to in2[N] . And after proper normalization, you can see that in2[N] is identical to in[N] .

 #include <stdlib.h> #include <math.h> #include <fftw3.h> #define N 16 int main(void) { fftw_complex in[N], out[N], in2[N]; /* double [2] */ fftw_plan p, q; int i; /* prepare a cosine wave */ for (i = 0; i < N; i++) { in[i][0] = cos(3 * 2*M_PI*i/N); in[i][1] = 0; } /* forward Fourier transform, save the result in 'out' */ p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute(p); for (i = 0; i < N; i++) printf("freq: %3d %+9.5f %+9.5f I\n", i, out[i][0], out[i][1]); fftw_destroy_plan(p); /* backward Fourier transform, save the result in 'in2' */ printf("\nInverse transform:\n"); q = fftw_plan_dft_1d(N, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE); fftw_execute(q); /* normalize */ for (i = 0; i < N; i++) { in2[i][0] *= 1./N; in2[i][1] *= 1./N; } for (i = 0; i < N; i++) printf("recover: %3d %+9.5f %+9.5f I vs. %+9.5f %+9.5f I\n", i, in[i][0], in[i][1], in2[i][0], in2[i][1]); fftw_destroy_plan(q); fftw_cleanup(); return 0; } 

This is the conclusion:

 freq: 0 -0.00000 +0.00000 I freq: 1 +0.00000 +0.00000 I freq: 2 -0.00000 +0.00000 I freq: 3 +8.00000 -0.00000 I freq: 4 +0.00000 +0.00000 I freq: 5 -0.00000 +0.00000 I freq: 6 +0.00000 -0.00000 I freq: 7 -0.00000 +0.00000 I freq: 8 +0.00000 +0.00000 I freq: 9 -0.00000 -0.00000 I freq: 10 +0.00000 +0.00000 I freq: 11 -0.00000 -0.00000 I freq: 12 +0.00000 -0.00000 I freq: 13 +8.00000 +0.00000 I freq: 14 -0.00000 -0.00000 I freq: 15 +0.00000 -0.00000 I Inverse transform: recover: 0 +1.00000 +0.00000 I vs. +1.00000 +0.00000 I recover: 1 +0.38268 +0.00000 I vs. +0.38268 +0.00000 I recover: 2 -0.70711 +0.00000 I vs. -0.70711 +0.00000 I recover: 3 -0.92388 +0.00000 I vs. -0.92388 +0.00000 I recover: 4 -0.00000 +0.00000 I vs. -0.00000 +0.00000 I recover: 5 +0.92388 +0.00000 I vs. +0.92388 +0.00000 I recover: 6 +0.70711 +0.00000 I vs. +0.70711 +0.00000 I recover: 7 -0.38268 +0.00000 I vs. -0.38268 +0.00000 I recover: 8 -1.00000 +0.00000 I vs. -1.00000 +0.00000 I recover: 9 -0.38268 +0.00000 I vs. -0.38268 +0.00000 I recover: 10 +0.70711 +0.00000 I vs. +0.70711 +0.00000 I recover: 11 +0.92388 +0.00000 I vs. +0.92388 +0.00000 I recover: 12 +0.00000 +0.00000 I vs. +0.00000 +0.00000 I recover: 13 -0.92388 +0.00000 I vs. -0.92388 +0.00000 I recover: 14 -0.70711 +0.00000 I vs. -0.70711 +0.00000 I recover: 15 +0.38268 +0.00000 I vs. +0.38268 +0.00000 I 

Answer 3

Very useful for magnify2x image lines.

As in the following example, increase signal 2 of direct signal

 int main (void) { fftw_complex in[N], out[N], in2[N2], out2[N2]; /* double [2] */ fftw_plan p, q; int i; int half; half=(N/2+1); /* prepare a cosine wave */ for (i = 0; i < N; i++) { //in[i][0] = cos ( 3 * 2 * M_PI * i / N); in[i][0] = (i > 3 && i< 12 )?1:0; in[i][1] = (i > 3 && i< 12 )?1:0; //in[i][1] = 0; } /* forward Fourier transform, save the result in 'out' */ p = fftw_plan_dft_1d (N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute (p); for (i = 0; i < N; i++) printf ("input: %3d %+9.5f %+9.5f I %+9.5f %+9.5f I\n", i, in[i][0], in[i][1],out[i][0],out[i][1]); fftw_destroy_plan (p); for (i = 0; i<N2; i++) {out2[i][0]=0.;out2[i][1]=0.;} for (i = 0; i<half; i++) { out2[i][0]=2.*out[i][0]; out2[i][1]=2.*out[i][1]; } for (i = half;i<N;i++) { out2[N+i][0]=2.*out[i][0]; out2[N+i][1]=2.*out[i][1]; } /* backward Fourier transform, save the result in 'in2' */ printf ("\nInverse transform:\n"); q = fftw_plan_dft_1d (N2, out2, in2, FFTW_BACKWARD, FFTW_ESTIMATE); fftw_execute (q); /* normalize */ for (i = 0; i < N2; i++) { in2[i][0] /= N2; in2[i][1] /= N2; } for (i = 0; i < N2; i++) printf ("recover: %3d %+9.1f %+9.1f I\n", i, in2[i][0], in2[i][1]); fftw_destroy_plan (q); fftw_cleanup (); return 0; } 

Answer 4

Here is what I did. Mine is specifically designed to get the same output as Matlab. This only works with a column matrix (you can replace AMatrix with std::vector ).

 AMatrix<complex<double>> FastFourierTransform::Forward1d(const AMatrix<double> &inMat) { AMatrix<complex<double>> outMat{ inMat.NumRows(), inMat.NumCols() }; size_t N = inMat.NumElements(); bool isOdd = N % 2 == 1; size_t outSize = (isOdd) ? ceil(N / 2 + 1) : N / 2; fftw_plan plan = fftw_plan_dft_r2c_1d( inMat.NumRows(), reinterpret_cast<double*>(&inMat.mutData()[0]), // mutData here is as same v.data() for std::vector reinterpret_cast<fftw_complex*>(&outMat.mutData()[0]), FFTW_ESTIMATE); fftw_execute(plan); // mirror auto halfWayPoint = outMat.begin() + (outMat.NumRows() / 2) + 1; auto startCopyDest = (isOdd) ? halfWayPoint : halfWayPoint - 1; std::reverse_copy(outMat.begin() + 1, halfWayPoint, startCopyDest); // skip DC (+1) std::for_each(halfWayPoint, outMat.end(), [](auto &c) { return c = std::conj(c);}); return std::move(outMat); } AMatrix<complex<double>> FastFourierTransform::Forward1d(const AMatrix<double> &inMat, size_t points) { // append zeros to matrix until it the same size as points AMatrix<double> sig = inMat; sig.Resize(points, sig.NumCols()); for (size_t i = inMat.NumRows(); i < points; i++) { sig(i, 0) = 0; } return Forward1d(sig); } AMatrix<complex<double>> FastFourierTransform::Inverse1d(const AMatrix<complex<double>> &inMat) { AMatrix<complex<double>> outMat{ inMat.NumRows(), inMat.NumCols() }; size_t N = inMat.NumElements(); fftw_plan plan = fftw_plan_dft_1d( inMat.NumRows(), reinterpret_cast<fftw_complex*>(&inMat.mutData()[0]), reinterpret_cast<fftw_complex*>(&outMat.mutData()[0]), FFTW_BACKWARD, FFTW_ESTIMATE); fftw_execute(plan); // Matlab normalizes auto normalize = [=](auto &c) { return c *= 1. / N; }; std::for_each(outMat.begin(), outMat.end(), normalize); fftw_destroy_plan(plan); return std::move(outMat); } // From Matlab documentation: "ifft tests X to see whether vectors in X along the active dimension // are conjugate symmetric. If so, the computation is faster and the output is real. // An N-element vector x is conjugate symmetric if x(i) = conj(x(mod(N-i+1,N)+1)) for each element of x." // http://uk.mathworks.com/help/matlab/ref/ifft.html AMatrix<double> FastFourierTransform::Inverse1dConjSym(const AMatrix<complex<double>> &inMat) { AMatrix<double> outMat{ inMat.NumRows(), inMat.NumCols() }; size_t N = inMat.NumElements(); fftw_plan plan = fftw_plan_dft_c2r_1d( inMat.NumRows(), reinterpret_cast<fftw_complex*>(&inMat.mutData()[0]), reinterpret_cast<double*>(&outMat.mutData()[0]), FFTW_BACKWARD | FFTW_ESTIMATE); fftw_execute(plan); auto normalize = [=](auto &c) { return c *= 1. / N; }; std::for_each(outMat.begin(), outMat.end(), normalize); fftw_destroy_plan(plan); return std::move(outMat); }