C Programming: The reconstruction of the Fast Fourier transform of a real signal using fftw libraries

This post is to help any one who is implementing the Fast Fourier transform on a real sequence of numbers or one dimensional signal in the C programming language. The code provided performs a one dimensional FFT and calculates the magnitude and phase from the real and imaginary components, it then performs the inverse Fourier transform to extract the original sequence of numbers from the magnitude and phase. This process is often used for signal processing applications such as speech enhancement or any kind of audio modification. More information about the Fast Fourier Transform can be found at Wikipedia.

To use this code, ensure you download the FFTW libraries

Create a new c program called test_fftw.c and ensure you compile with a link to the FFTW libraries by using -lfftw3, for example:

gcc -g -lfftw3 -lm fftw_test.c -o fftw_test

Here is the source code for using the FFTW libraries to calculate the fast Fourier transform:

#include 
#include 
#include 
#include 
#include 
#include 

int main(void)
{
    // Variable Declaration
    double array[] = {0.1, 0.6, 0.1, 0.4, 0.5, 0, 0.8, 0.7, 0.8, 0.6, 0.1,0};
    double *out,*mag,*phase;
    double real,imag;
    int i,j;
    int size = 12;
    fftw_complex *out_cpx, *mod_cpx;
    fftw_plan fft; 
    fftw_plan ifft; 
    
    //Allocate Memory
    out_cpx = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(size/2+1));
    mod_cpx = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(size/2+1));
    out = (double *) malloc(size*sizeof(double));
    mag = (double *) malloc(size*sizeof(double));
    phase = (double *) malloc(size*sizeof(double));
   
    fft = fftw_plan_dft_r2c_1d(size, array, out_cpx, FFTW_ESTIMATE);  //Setup fftw plan for fft 1 dimensional, real signal
    ifft = fftw_plan_dft_c2r_1d(size, mod_cpx, out, FFTW_ESTIMATE);   //Setup fftw plan for ifft 1 dimensional, complex signal

    fftw_execute(fft);	//perform fft
	for(j=0;j<size/2+1;j++)
	{
        	real = out_cpx[j][0];	//Extract real component
		imag = out_cpx[j][1];   //Extract imaginary component
		mag[j] = sqrt((real*real)+(imag*imag));  // Calculate the Magnitude
		phase[j] = atan2(imag,real); // Calculate the Phase
	}
	
	/***********MODIFICATION****************************/
	//You can perform frequency domain modification here  	
	/***************************************************/	

	for(j=0;j<size/2+1;j++)
	{
		mod_cpx[j][0] = (mag[j]*cos(phase[j]));  //Construct new real component
		mod_cpx[j][1] = (mag[j]*sin(phase[j]));  //Construct new imaginary  component
	}

    fftw_execute(ifft); //perform ifft
 
    // Print input and output
    printf("Input:    tOutput:n");
    for(i=0;i<size;i++)
    {
	out[i] = out[i]/size;	
	printf("%ft%fn",(array[i]),out[i]);
    }

    // Free all memory
    fftw_destroy_plan(fft);
    fftw_destroy_plan(ifft);
    fftw_free(out_cpx);
    fftw_free(mod_cpx);
    free(out);
    free(mag);
    free(phase);
    return 0;
}

Which outputs:

Input:    	Output:
0.100000	0.100000
0.600000	0.600000
0.100000	0.100000
0.400000	0.400000
0.500000	0.500000
0.000000	0.000000
0.800000	0.800000
0.700000	0.700000
0.800000	0.800000
0.600000	0.600000
0.100000	0.100000
0.000000	0.000000

You can also calculate FFT using the light weight kiss fft libraries, however it doesn’t have the same support as the FFTW libraries and is not optimised.