Discrete Fourier Transform using FFTW

If we have a time series data and we want to know the frequency content of it, then we need to perform Fourier transformation. This can be done numerically using available mathematical softwares such as Matlab, Mathematica and Maple. For C programming language, I mostly use FFTW since it is free and easy to use on linux OS and in Cygwin on Windows.

To understand how to do DFT, let us a look at Fourier transform of a square function given by f(t) = 1 for |t-2|<=1 and f(t) = 0 for |t-2|>1. We consider the function at interval (0,4). First discretize the data into grids with spacing 0.1.

--- try_fftw.c ---

#include <stdio.h>
#include <math.h>
#include <fftw3.h>

#define N 40

int main(void)
{
int i;
double dt, t;
fftw_complex *f, *g;
fftw_plan p;

dt = 0.1;
/* allocate memory for f and g */
f = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
g = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

/* compute square function f(x) */
for(i=0;i<N;i++){
t = i*dt;
if(fabs(t-2)<=1){
f[i][0] = 1;
f[i][1] = 0;
}
else{
f[i][0] = 0;
f[i][1] = 0;
}
}

/* Print f(t) */
printf(" -- i f(t) -- \n");
for(i=0;i<N;i++){
printf("%d %e \n", i, f[i][0]);
}

/* create plan */
p = fftw_plan_dft_1d(N, f, g, FFTW_FORWARD, FFTW_ESTIMATE);

/* perform FFT */
fftw_execute(p);

/* Print results */
printf(" -- i g(w) -- \n");
for(i=0;i<N;i++){
printf("%d %e %e \n", i, g[i][0], g[i][1]);
}

/* Free memory */
fftw_destroy_plan(p);
fftw_free(f);
fftw_free(g);

return 0;
}


-------

To compile use:

 gcc try_fftw.c  -lm -lfftw3 -o try  

Results are shown in the figures below





Ref: FFTW3 manual