GSL Fast-Fourier Transform - Nonsense Output

Question

The Fourier transform of a gaussian is a gaussian, but for some reason the fast Fourier transform library from GSL (GNU scientific library) doesn't give this at all. I've included the code I've used to generate the (attempted) Fourier transform, and two relevant plots right after it. Could help me identify what I've messed up?

#include <gsl/gsl_fft_complex.h>
#include <fstream>

#define REAL(z,i) ((z)[2*(i)]) //complex arrays stored as    
#define IMAG(z,i) ((z)[2*(i)+1])

using namespace std;

int main(){

double N = pow(2,9); //power of 2 for Cooley-Tukey algorithm
int n = (int) N;

double f[2*n];
double dx = 10./N;
double x = -5.;
ofstream fileo("out.txt");

for (int i=0; i<n; ++i){      //initialize gaussian
    REAL(f,i)=exp(-0.5*x*x);  
    IMAG(f,i)=0.;
    x+=dx;
   }

   gsl_fft_complex_radix2_forward(f, 1, n);  //Fourier transform

   for (int i=0; i<n; ++i){
        fileo<<i<<" "<<REAL(f,i)<<'\n';  //plot frequency distribution
   }

   fileo.close();
}

EDIT: Solved!

As stated in @roadrunner66's answer, the width of the original Gaussian was very broad, leading to a ridiculously narrow Gaussian in Fourier space. Moreover, my plot looked funky because, as suggested in @n.m's (now removed) comment, the Fourier transform returns the DFT with projections onto k-values indexed as k=0,1,...,N/2,-N/2,...-2,-1.

roadrunner66 roadrunner66 · Accepted Answer · 2016-05-02T05:07:59

It looks good to me. Shift the output vector by N/2 and plot the absolute value of the output, rather than the real part.

Also note that your input Gaussian is rather wide, which makes it's spectrum very narrow. Check the analytical solution for that case for comparison.

GSL Fast-Fourier Transform - Nonsense Output

1 Answers