Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

GSL Fast-Fourier Transform - Nonsense Output

Tags:

c++

fft

gsl

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();
}

enter image description here

enter image description here


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.

enter image description here

like image 218
Arturo don Juan Avatar asked May 02 '16 04:05

Arturo don Juan


Video Answer


1 Answers

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.

like image 171
roadrunner66 Avatar answered Oct 14 '22 00:10

roadrunner66