how to do spectrum analyse with FFTW?

Question

I would like to better understand FFTW's API. FFTW is a library for computing the discrete Fourier transform (DFT) in one or more dimensions.

Now, assume I have a sinus waveform x=30*sin(2*M_PI*f*i*T), where f is the frequency (for example f=1000Hz). If I use FFTW's function to analyse my waveform, I expect to get one frequency f=1000Hz.

My question is how can I do this in c++ using FFTW Library? Any help would be appreciated.

Answer 1

You can find more details in FFTW's documentation.

However, for the relatively simple case of one-dimensional real-valued signal, the following is a sum-up of the general steps you have to do.

Typically you would need to allocate input/output buffers, and a data-structure which FFTW uses for its own bookkeeping which the library refers as a plan. This can be done in a number of ways (more details in FFTW's documentation), such as:

  #include "fftw3.h"

  // First choose a buffer size:
  //   Typically best performance with a power of 2
  //   but could be a product of small primes
  int           input_size    = 1024; 
  //   Compute corresponding number of complex output samples
  int           output_size   = (input_size/2 + 1);

  // Allocate input and output buffers
  double*       input_buffer  = static_cast<double*      >(fftw_malloc(input_size  * sizeof(double)));
  fftw_complex* output_buffer = static_cast<fftw_complex*>(fftw_malloc(output_size * sizeof(fftw_complex)));

  // Create plan
  //   Select plan creation flags
  //   see http://www.fftw.org/fftw3_doc/Planner-Flags.html#Planner-Flags
  int           flags = FFTW_ESTIMATE;
  fftw_plan     plan = fftw_plan_dft_r2c_1d(input_size, 
                                            input_buffer, 
                                            output_buffer, 
                                            flags);

Once that is done, you can fill in the input_buffer with the real-valued data samples to analyse, and perform the FFT using:

 fftw_execute(plan);

The complex-valued results will be stored in output_buffer, where output_buffer[0] corresponds to a frequency of 0, and output_buffer[output_size-1] corresponds to half the sampling rate. This plan can be executed multiple times (with updated values in input_buffer, resulting in correspondingly updated values in output_buffer).

Note that typically fftw_complex (which is the datatype used for the output in this example) is implemented as an array of 2 values: index 0 corresponds to the real-part, and index 1 corresponds to the imaginary-part (e.g. output_buffer[i][0] corresponds to the real-part of the i^th frequency component).

Once you are done, you can release allocated resources with:

  fftw_free(input_buffer);
  fftw_free(output_buffer);
  fftw_destroy_plan(plan);

Note that if you may use float, double or long double versions of those functions. Simply link against the corresponding libfftw3f-3.lib, libfftw3-3.lib or libfftw3l-3.lib.

Update: If you want to use complex-valued input samples together with fftw_plan_dft_1d, you will then have to set the real and imaginary parts like so:

for (i = 0; i < N-1; ++i) {
  t[i]=i*T;
  signal[i][0] = 0.7 * sin(2*M_PI*f*t[i]); // real-part
  signal[i][1] = 0.0; // imaginary-part
}

Alternatively change the input sample type to float, double or long double (together with using fftw_plan_dft_r2c_1d).