Plotting an input curve and its FFT using gnuplot in C++ I'm currently working on a project1 which will result in millions of FFTs being solved time and time again. I have a finished prototype in Python (using numpy, scipy, and the pyfftw3 package), but there are significant bottlenecks in calculation speed. I have been writing a GPU version using arrayfire, but there's still a large amount of number crunching which can't be done in parallel. I'm hoping that translating the project to C++, either the entire thing or just the compute-intensive parts, will provide an extra boost over Python where the GPU isn't appropriate.

As a starting point, I have written a very basic C++ script which defines a Gaussian curve, takes the FFT using the FFTW3 library, and plots the input curve with its FFT using gnuplot. As this is my first script in C++ (neglecting 'hello world') I'd like to make sure I'm not wandering severely off track when it comes to optimizing efficiency and speed. I haven't benchmarked this yet, but was hoping that someone with more experience in using these libraries and the language could scan over it and point out any howlers.

The code:

#include "gnuplot-iostream.h"
#include &lt;iostream&gt;
#include &lt;vector&gt;
#include &lt;cmath&gt;
#include &lt;fftw3.h&gt;
#include &lt;complex&gt;

typedef std::complex&lt;double&gt; doubleComp;

int main()
{
    Gnuplot gp;

    doubleComp i;
    i = -1;
    i = sqrt(i);

    int gridOrder = 9;
    int nPoints = pow(2, gridOrder);
    int halfPoints = pow(2, gridOrder-1);
    float tau = 18;
    float inTemp, outTemp;

    fftw_complex in[nPoints], out[nPoints], pArr1[nPoints], pArr2[nPoints];
    fftw_plan p;

    for (int idx=0; idx&lt;nPoints; idx++) {
        in[idx][0] = exp(-pow((idx-halfPoints)/tau, 2));
        in[idx][1] = 0;
    }

    p = fftw_plan_dft_1d(nPoints, pArr1, pArr2, FFTW_FORWARD, FFTW_MEASURE);
    fftw_execute_dft(p, in, out);

    for (int idx=0; idx&lt;nPoints; idx++) {
        out[idx][0] *= (1/double(nPoints));
        out[idx][1] *= (1/double(nPoints));
    }

    std::vector&lt;std::pair&lt;double, double&gt; &gt;inPlot;
    std::vector&lt;std::pair&lt;double, double&gt; &gt;outPlot;
    for (int idx=0; idx&lt;nPoints; idx++) {
        inTemp = std::abs(in[idx][0] + i*in[idx][1]);

        // Deal with fftshifts needed for visuals:
        if (idx &gt; halfPoints)
            outTemp = std::abs(out[idx-halfPoints][0] +
                                i*out[idx-halfPoints][1]);
        if (idx &lt; halfPoints)
            outTemp = std::abs(out[idx+halfPoints][0] +
                                i*out[idx+halfPoints][1]);

        inPlot.push_back(std::make_pair(idx-halfPoints, inTemp));
        outPlot.push_back(std::make_pair(idx-halfPoints, outTemp));
    }

    gp &lt;&lt; "set xlabel 'Grid points' font 'Helvetica, 12'\n" &lt;&lt; 
        "set ylabel 'Intensity' font 'Helvetica, 12'" &lt;&lt; std::endl;
    gp &lt;&lt; "set tics font 'Helvetica, 10'" &lt;&lt; std::endl;
    gp &lt;&lt; "set title 'Gaussian pulse example' font 'Helvetica, 13'" &lt;&lt; std::endl;
    gp &lt;&lt; "set grid xtics lt 0 lw 1 lc 'gray'\n" &lt;&lt; 
        "set grid ytics lt 0 lw 2 lc 'gray'" &lt;&lt; std::endl;

    gp &lt;&lt; "set term wxt '0'" &lt;&lt; std::endl;
    gp &lt;&lt; "plot '-' lt 1 lw 2 linecolor 'blue' with line title 'Time domain'\n";
    gp.send1d(inPlot);

    gp &lt;&lt; "set term wxt '1'" &lt;&lt; std::endl;
    gp &lt;&lt; "plot '-' lt 11 lw 2 with line title"
        "'Spectral domain'\n";
    gp.send1d(outPlot);

    fftw_destroy_plan(p);

    return 0;
}


To compile (Ubuntu 16.04):

g++ -O3 GaussianPulseDefinition.cpp -o GaussianPulseDefinition -lboost_system -lboost_iostreams -lboost_filesystem -lfftw3 -lm

1 For those who like to know the application: The project is basically a coupled PDE solver (with Qt GUI) aimed at simulating the nonlinear dynamics of optical pulse propagation in optical fibre, fibre laser systems, optical amplifiers, etc.
Firstly, does the program give correct results?

If so, how do you know?

Without going into extensive unit testing, checking your results against a Matlab prototype/your Python implementation is probably good enough for now.

I can't comment on the Gnuplot functions as I don't have it installed, also but some general C++ tips are...


std::vectors are your friend! Try use them instead of C-style arrays. They 'know' their size and accessing out-of-bounds values are a major problem with C-style arrays and definitely for fiddly calculations like FFT. (You can use my_vector.data() to pass it into C-style functions).
Reduce reliance on a 'global' nPoints constant by using the vectors .size() member. (Unfortunately this outputs an unsigned integer size_t so there are some casts in the code to silence the compiler warnings). This is so that if nPoints goes out of sync with the size of the array/vector, you're not going to access an out-of-bounds value.
Always try to define a variable's value as soon as you declare it. Undefined variables are a major problem as C++ doesn't (usually) auto-define say, a double to zero. (double x; would refer to some garbage value in memory).
Always try to define variables when you need them i.e. as late as possible (see code for examples)
Try get into the habit of spitting out little functions right from the get-go, instead of one big block of code. It will break the problem up and be easier to debug. Also this will limit some variables scope so its easier to see what parts of the program change what variables.
For FFTs probably best to stick to doubles all round for the extra precision (removed the floats)
Use C++11 if you can and turn up the compiler warning level g++ -std=c++11 -Wall ...


Hopefully a bit to get you started!

I whole heartedly recommend this book if you're coming from Python, it's a short in length but really comprehensive and does a much better job than I could at explaining how to write good C++, and by the man himself - Stroustrup's "A Tour Of C++"

//#include "gnuplot-iostream.h"
#include &lt;iostream&gt;
#include &lt;vector&gt;
#include &lt;cmath&gt;
#include &lt;fftw3.h&gt;
#include &lt;complex&gt;
#include &lt;cassert&gt;

using doublePairs = std::vector&lt;std::pair&lt;double, double&gt; &gt;; // using (instead of typedef) for long named vector of pairs 

void makeGaussian(std::vector&lt;fftw_complex&gt;&amp; v, double tau) { // pass-by-reference indicates v is going to be altered
    for (size_t idx=0; idx&lt;v.size(); idx++) {   // in.size() so we don't go out of range
        const int halfPoints = v.size() / 2;
        v[idx][0] = exp(-pow((idx-halfPoints)/tau, 2));
        v[idx][1] = 0;
    }
}

void performFft(std::vector&lt;fftw_complex&gt;&amp; in,
                std::vector&lt;fftw_complex&gt;&amp; out) {

    fftw_plan p;

    assert(in.size() == out.size()); // test same sizes
    const int nPoints = static_cast&lt;int&gt; (in.size());

    std::vector&lt;fftw_complex&gt; pArr1 (nPoints);
    std::vector&lt;fftw_complex&gt; pArr2 (nPoints);
    p = fftw_plan_dft_1d(nPoints, pArr1.data(), pArr2.data(), FFTW_FORWARD, FFTW_MEASURE);
    fftw_execute_dft(p, in.data(), out.data()); // access raw pointers to vectors for c-style function call

    fftw_destroy_plan(p);

    for (auto&amp; n : out) {                           // nicer c++11 range based for loop
        n[0] *= (1/static_cast&lt;double&gt; (nPoints));  // static_cast&lt;T&gt; is preferred to c-style casts
        n[1] *= (1/static_cast&lt;double&gt; (nPoints));
    }
}

void makePlotPairArrays(const std::vector&lt;fftw_complex&gt;&amp; in,
                        const std::vector&lt;fftw_complex&gt;&amp; out,
                        doublePairs&amp; inPlot, 
                        doublePairs&amp; outPlot) {

    assert(in.size() == out.size()); // test same sizes
    const std::complex&lt;double&gt; i {std::sqrt(-1)}; // typedef not necessary? (it's only used once)
    const int halfPoints = in.size() / 2;

    for (size_t idx=0; idx&lt;in.size(); idx++) {
        // these Temps should be doubles?? since it's a pair of doubles not floats?
        // would lead to less precise FFTs
        double inTemp = std::abs(in[idx][0] + i*in[idx][1]); // keep these local to the loop scope

        // Deal with fftshifts needed for visuals:
        double outTemp {0};

        const int idx_ = static_cast&lt;int&gt; (idx); // for signed/unsigned compiler warnings
        if (idx_ &gt; halfPoints)                              // what about when idx == halfPoints ???
            outTemp = std::abs(out[idx-halfPoints][0] +     // try make this an if/else statement so
                                i*out[idx-halfPoints][1]);  // at least one branch is always touched
        if (idx_ &lt; halfPoints)
            outTemp = std::abs(out[idx+halfPoints][0] +
                                i*out[idx+halfPoints][1]);

        inPlot.push_back(std::make_pair(idx-halfPoints, inTemp));
        outPlot.push_back(std::make_pair(idx-halfPoints, outTemp));
    }
}

int main()
{
    constexpr int gridOrder = 9;
    constexpr int nPoints = std::pow(2, gridOrder);
    constexpr double tau = 18;
    // 'halfPoints' not needed here by declaring in function's scope
    // std::pow(2, gridOrder-1); BTW pow is heavy why not nPoints / 2?

    // our fft arrays
    std::vector&lt;fftw_complex&gt; in (nPoints);
    std::vector&lt;fftw_complex&gt; out (nPoints);

    makeGaussian(in, tau);

    performFft(in, out);

    doublePairs inPlot;
    doublePairs outPlot;
    makePlotPairArrays (in, out, inPlot, outPlot);

    // ... do the GnuPlotting

    return 0;
}