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
#include
#include
#include
#include

typedef std::complex 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 >inPlot;
std::vector >outPlot;
for (int idx=0; idx halfPoints)
outTemp = std::abs(out[idx-halfPoints][0] +
i*out[idx-halfPoints][1]);
if (idx < 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(st

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 
#include 
#include 
#include 
#include 
#include 

using doublePairs = std::vector >; // using (instead of typedef) for long named vector of pairs 

void makeGaussian(std::vector& v, double tau) { // pass-by-reference indicates v is going to be altered
    for (size_t idx=0; idx& in,
                std::vector& out) {

    fftw_plan p;

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

    std::vector pArr1 (nPoints);
    std::vector 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& n : out) {                           // nicer c++11 range based for loop
        n[0] *= (1/static_cast (nPoints));  // static_cast is preferred to c-style casts
        n[1] *= (1/static_cast (nPoints));
    }
}

void makePlotPairArrays(const std::vector& in,
                        const std::vector& out,
                        doublePairs& inPlot, 
                        doublePairs& outPlot) {

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

    for (size_t idx=0; idx (idx); // for signed/unsigned compiler warnings
        if (idx_ > 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_  in (nPoints);
    std::vector out (nPoints);

    makeGaussian(in, tau);

    performFft(in, out);

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

    // ... do the GnuPlotting

    return 0;
}

StackExchange Code Review Q#154657, answer score: 5