FFTW3 : C++ Implementation on MATLAB/OCTAVE Perspective

I used to code using MATLAB and OCTAVE for my signal processing research. But, when it comes to real implementation and performance, I always stop and wonder how to make my concept coded in C/C++. Moreover, my MATLAB license is expired. All I need is Fourier Transform because it is the basic operation for signal processing. Then I made research on how Fourier Transform could be done in C/C++. I ended up on Fast Fourier Transform in the West (FFTW). It is considered the best and fastest FFT implementation in C/C++.

In this article, I would like to demonstrate basic forward and inverse transform of monotone signal. Because of we use audo signal in assumption, so this is 1D FFT problem. Well I think this simple example could be an entrance gate for MATLAB and OCTAVE programmers to implement their concepts on C/C++. I hope this article can be helpful for us.

MATLAB/Octave Point of View

Assume we have 2 seconds duration of 440 Hz waveform. We are going to do transformation on this signal, do some operation, and transform it back to time-domain signal. For this example, I am intentionally left the do some operation blank. The implementaion could be like this :

% Generating waveform
fs = 8000;
t = linspace(0, 1, 2 * fs);
x = sin(2 * pi * 440 * t);

% Doing forward transformation
spectrum = fft(x);

% do some operation
% left blank

% Doing inverse transform
y = ifft(spectrum);

These codes are typical way to perform FFT in MATLAB. We generate 16000 samples sinusoidal signal in 8000 Hz which mean 2 seconds length. Then we transform it into frequency domain, perform some operation, and transform it back to time domain.

C++ Implementation

• Step 0 : Installing the library
  If you don’t care about serial or parallel computation, let’s make it serial. You can install FFTW3 through terminal as follows.

 sudo apt-get install libfftw3-dev libfftw3-doc

As far as I know, this repository only support serial computation. If parallelization realy matters, you will need to follow steps in this link.

• Step 1 : The code
  Create our source code, for instance fourier.cc.

#include <cstdio>
#include <cmath>
#include <fftw3.h>

#define PI 3.14159265359
#define REAL 0
#define IMAG 1 

int main(int argc, char const *argv[]) {
    // Generating waveform
    int fs = 8000;  // sampling frequency
    size_t wave_duration = 2 * fs;  // 2 seconds wave
    fftw_complex x[wave_duration];
    for (size_t i = 0; i < wave_duration; i++) {
        x[i][REAL] = sin (2.0 * PI * 440.0 * (float)i / (float)fs );
        x[i][IMAG] = 0;
    }

    // Doing forward transform
    fftw_complex spectrum[wave_duration];
    fftw_plan fft = fftw_plan_dft_1d(wave_duration, x, spectrum,
        FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(fft);

    // Doing inverse tranform
    fftw_complex y[wave_duration];
    fftw_plan ifft = fftw_plan_dft_1d(wave_duration, spectrum, y,
        FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(ifft);
    for (size_t i = 0; i < wave_duration; i++) {
        y[i][REAL] /= 2 * fs;
    }

    // Evaluating result
    printf("time\torigin\tinverse\tdelta\n");
    for (size_t i = 0; i < wave_duration; i++) {
        double t = (float) i / (float) fs; 
        double delta = x[i][REAL] - y[i][REAL];
        printf("%.4f\t%.4f\t%.4f\t%.12f\n", t, x[i][REAL], y[i][REAL], delta);
    }

    // Cleaning memory
    fftw_destroy_plan(fft);
    fftw_destroy_plan(ifft);
    return 0;
}

• Step 2 : header and definition

#include <cstdio>
#include <cmath>
#include <fftw3.h>

#define PI 3.14159265359
#define REAL 0
#define IMAG 1

We are including two standard libraries, cstdio for printing into console and cmath for mathematic operation; and fftw3.h library. We also define our constants for PI, array index for REAL number, and array index IMAG number.

• Step 3 : Generating waveform

    // Generating waveform
    int fs = 8000;  // sampling frequency
    size_t wave_duration = 2 * fs;  // 2 seconds wave
    fftw_complex x[wave_duration];
    for (size_t i = 0; i < wave_duration; i++) {
        x[i][REAL] = sin (2.0 * PI * 440.0 * (float)i / (float)fs );
        x[i][IMAG] = 0;
    }

We declare the container of our signal using fft_complex. Both for forward transform and inverse transform container use fft_complex. In fftw, complex array are declared as two dimension array. In other words, fftw_complex is simply the same as c++ code below

double complex[N][2];

where N is the array length, real number is in index 0, and imaginary number is in index 1.

Our signal is obtained according this equation
$$x(index) = sin(2\cdot\pi\cdot f\cdot\frac{index}{f_{s}})$$

We work on real component only and the imaginary component is left blank.

• Step 4 : Forward transform

    // Doing forward transform
    fftw_complex spectrum[wave_duration];
    fftw_plan fft = fftw_plan_dft_1d(wave_duration, x, spectrum,
        FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(fft);

To do transformation, we should make a plan. With 1D Fourier problem, we have to use fftw_plan_dft_1d function. This function requires 5 inputs. The first input is the array length, the second is input complex array, the third is output complex array, transform sign, and flag.

• Step 5 : Inverse transform

    // Doing inverse tranform
    fftw_complex y[wave_duration];
    fftw_plan ifft = fftw_plan_dft_1d(wave_duration, spectrum, y,
        FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(ifft);
    for (size_t i = 0; i < wave_duration; i++) {
        y[i][REAL] /= 2 * fs;
    }

We need to declare inverse transform result container, for instance we give it name y. The only difference in fftw_plan for inverse transform is the transformation sign which is FFTW_BACKWARD. Then we execute the plan.

One that must be our attention is the result must be normalized first. The amplitude of inverted signal is magnified as big as two times sampling frequency sampling. Thus, the normalization is done by dividing all values by 2 * fs. This operation is only conducted for real component because if our spectrum is perfectly symmetric, the imaginary component will be zero altogether.

• Step 6 : Evaluation

    // Evaluating result
    printf("time\torigin\tinverse\tdelta\n");
    for (size_t i = 0; i < wave_duration; i++) {
        double t = (float) i / (float) fs; 
        double delta = x[i][REAL] - y[i][REAL];
        printf("%.4f\t%.4f\t%.4f\t%.12f\n", t, x[i][REAL], y[i][REAL], delta);
    }

In this loop, we compare between the original monotone waveform and inverse transform wave. This evaluation is designated to find out if the result correct or not by subtracting these two signal. If the deltas are all zero, the results are correct.

• Step 7 : Compiling and executing
  Compile your code using g++ and add fftw3 library

g++ -std=c++11 -o fourier fourier.cc -lfftw3

Sample result (skipped) :

$ ./fourier
time    origin  inverse delta
0.0000  0.0000  0.0000  -0.000000000000
0.0001  0.0785  0.0785  -0.000000000000
0.0003  0.1564  0.1564  0.000000000000
0.0004  0.2334  0.2334  0.000000000000
0.0005  0.3090  0.3090  0.000000000000
0.0006  0.3827  0.3827  -0.000000000000
0.0008  0.4540  0.4540  -0.000000000000
0.0009  0.5225  0.5225  -0.000000000000
0.0010  0.5878  0.5878  0.000000000000
.
.
.
1.9990  -0.5878 -0.5878 -0.000000000000
1.9991  -0.5225 -0.5225 0.000000000000
1.9993  -0.4540 -0.4540 0.000000000000
1.9994  -0.3827 -0.3827 0.000000000000
1.9995  -0.3090 -0.3090 -0.000000000000
1.9996  -0.2334 -0.2334 -0.000000000000
1.9998  -0.1564 -0.1564 0.000000000000
1.9999  -0.0785 -0.0785 0.000000000000

As we see, the original signal and inverted signal are completely the same.

That’s all. We could use this pattern for more complex implementation of signal processing ideas. Although the C++ implementation is not as handy as MATLAB, it give us the insight of low level language. Hope you could easily translate your concept from MATLAB/Octave into C++.

---

Correct me if I am wrong
Sincerely,

Tirtadwipa Manunggal
tirtadwipa.manunggal@gmail.com