Computing FFT and IFFT with FFTW library C++

Question

I am trying to compute the FFT and then the IFFT just to try out if I can get the same signal back but I am not really sure how to accomplish it. This is how I do the FFT:

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

Consistence

It did but I am not really sure how to Interpret the result, how do I access the frequencies? – Consistence

Answer 1

Have you at least tried to read the more than decent documentation?

They have a whole tutorial thought out for you to get to know FFTW:

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

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

This page has the info you need for FFT vs IFFT (see Arguments->sign). The docs also say that input->FFT->IFFT->n*input. So you'll have to be careful to scale the data correctly.

Congdon

I have but I am just unable to interpret the output and then how to come back within the time domain. – Consistence

Answer 2

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

Second, it inverse Fourier transforms the array out[N] back to in2[N]. And after a 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 output:

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

Predictory

Answer 3

Here's what I made. Mine is designed specifically to give the same output as Matlab. This only works on a column matrix (you can replace AMatrix with a 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's 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);
}

Groundsel

Your answer is the best among all except that we need to implement AMatrix and FastFourierTransform classes. It would be great if you can add those classes. (I will try to implement those but for other people, it would be very helpful/useful) – Amling

Answer 4

Very useful for magnify2x line of image

As in the following example magnify by factor 2 rect 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;
}

Macbeth