Calculating a real-number FFT with FFTW3

Calculating a real-number FFT with FFTW3

In this article I’m going to explain how to calculate a real-number based FFT using FFTW, providing the source code. First, though, I will introduce some important concepts.

Then I will provide code with comments within it. I will also provide a header containing sample data from a single channel sound file whose frames have been converted into double values. The samples are of a 440 Hz sine tone mixed with noise.

I exported the raw sound data from Audacity, and used the od program (GNU coreutils), to convert the raw data to double precision floating point numbers, and outputted that into a file that I could reformat into the header, “a_with_noise.h”.

Acknowledgements

Thanks Jake for pointing me in the right direction, providing sample code for computing a complex number FFT.

Assumptions

  1. Programming explanations and examples are for the C programming language.
  2. You understand what a FFT is used for.
  3. You have installed FFTW3 on your system.
  4. You are using a GNU/Linux distribution.

Precursory Information

Sampling Rate

This value is used to help determine the frequency bin size. For the sample data, the audio was originally sampled at 44100 times per second or 44100 Hz.

N

In the source code I’ve set ‘N’ to 4096, because it’s evenly divisible by two.

Frequency Bin Size

This is equal to the sampling rate divided by ‘N’, or 44100 divided by 4096.

Complex Numbers

These numbers consist of two parts: a real and imaginary component. By splitting them into two, it allows us to perform algebraic operations on it that normally couldn’t work with real numbers.

In C, we have a standard header, “complex.h”, which defines a complex type.

Sample Code

fft.c

// fft.c
// Written by Ben Cottrell
// This is licensed under CC-BY-SA-4.0

#include <complex.h>
#include <stdio.h>
#include <math.h>
#include <fftw3.h>
#include "a_with_noise.h"

// https://github.com/audacity/audacity/blob/c5ebc396eb06857b4509101fdd2b0620dc0658b3/src/FFT.cpp#L344
void hamming(double* in, int num_samples)
{
    const double multiplier = 2 * M_PI / num_samples;
    static const double coeff0 = 0.54, coeff1 = -0.46;
    for (int ii=0; ii < num_samples; ++ii) {
        in[ii] *= coeff0 + coeff1 * cos(ii * multiplier);
    }
}

int main() {
    int N;
    fftw_complex *out;
    // the optimal plan that fftw will use
    fftw_plan my_plan;
    // number of frames in awithnoise
    N=4096;
    // perform a hamming windowing function on our input data, smoothing extraneous frequencies
    hamming(awithnoise, N);
    out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*N);
    // create a plan to execute a DFT with real input and complex output, one dimension
    my_plan = fftw_plan_dft_r2c_1d(N, awithnoise, out, FFTW_ESTIMATE);
    // execute the plan
    fftw_execute(my_plan);
    double p = (double)N / 44100;
    double freq_bin = 0;
    double freq_bin_size = (double)44100 / N;
    double magnitude = 0;
    for (int i=0; i < N; i++) {
        // get the magnitude by calculating the square root of x^2 + y^2
        magnitude = sqrt(pow(creal(out[i]),2.0)+pow(cimag(out[i]),2.0));
        printf("%f %f\n", freq_bin, magnitude);
        freq_bin += freq_bin_size;
    }
    // destroy the plan, free memory
    fftw_destroy_plan(my_plan);
    fftw_free(out);
    return 0;
}

a_with_noise.h

// a_with_noise.h
// Written by Ben Cottrell
// This is licensed under CC-BY-SA-4.0

#ifndef A_WITH_NOISE_H
#define A_WITH_NOISE_H

double awithnoise[] =
{
        -0.10978938,
        0.06689477,
        -0.00797962,
        0.08235065,
        0.17361747,
        0.08064409,
        0.057826433,
        0.10899856,
        0.26787728,
        0.19625948,
        0.17439666,
        0.1664341,
        0.4255249,
        0.27694514,
        0.41154468,
        0.31954277,
        0.46462804,
        0.33654302,
        0.43988577,
        0.31898364,
        0.3791127,
...

Get the rest of the file from here: https://gist.githubusercontent.com/ben-cottrell-nz/d2a55a105fa49842da00291ddfa22625/raw/62e30668928ac5666da0510a67586c8f48d0bd49/a_with_noise.h

Building

Make sure you have built FFTW3 and installed it. When I tested this code, I used the CMake build system generator, but for the sake of brevity you could simply compile and test this code by executing the following text in your terminal emulator:

mkdir build 
cd build
gcc -o fft_test ../fft.c -I.. -lfftw3 -lm

If the compiler finished with no errors, run:

./fft_test

Here’s a CMakeLists.txt file as well:

cmake_minimum_required(VERSION 3.5)
project(fft_test LANGUAGES C)
set(CMAKE_C_STANDARD 11)
set(CMAKE_INCLUDE_CURRENT_DIR ON)
add_executable(fft_test
    fft.c)
target_link_libraries(fft_test
    m
    fftw3)

In Conclusion

I hope this gets you closer to being able to compute real-number FFT’s with FFTW3. If you have any questions about this, please leave them below.

Leave a Reply

Your email address will not be published. Required fields are marked *