Why does negative number show up when I do FFT on a Gaussian?

Question

I'm using this fftw library.

Currently I'm trying to plot a 2D Gaussian in the form e^(-(x^2+y^2)/a^2).

Here is the code:

using namespace std;
int main(int argc, char** argv ){
    fftw_complex *in, *out, *data;
    fftw_plan p;
    int i,j;
    int w=16;
    int h=16;
    double a = 2;
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*w*h);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*w*h);
    for(i=0;i<w;i++){
        for(j=0;j<h;j++){
            in[i*h+j][0] = exp(- (i*i+j*j)/(a*a));
            in[i*h+j][1] = 0;
        }
    }
    p = fftw_plan_dft_2d(w, h, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);
    //This is something that print what's in the matrix
    print_2d(out,w,h);

    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

Turns out negative numbers shows up. I thought Fourier transform of a Gaussian is another Gaussian, which shouldn't include any negative numbers.

Also, the current origin is at in[0]

Two possible problems: (i) you're doing a **discrete** FT, not an FT and (ii) your Gaussian function has effectively been multiplied by a rectangular window (since it is truncated), which is equivalent to a convolution in the frequency domain by the FT of the window function. — Paul R, Jun 24 '10 at 23:01
(iii) Floating-point math is notoriously imprecise, at least by naive standards. `-1E-30` would be a rounding error, not a true negative number. — MSalters, Jun 25 '10 at 08:50
I guess these are the reasons. Of course I can't check it, since I can't do continuous FT on my computer. — Chao Xu, Jun 25 '10 at 20:08

hbp · Answer 1 · 2021-12-29T04:25:07.273

EDIT: the previous answer is wrong, shifting the center of the Gaussian won't help as it introduces another phase shift. The right solution is to wrap high indices to negative ones:

double x = (i < w*0.5) ? i : (i - w);
double y = (j < h*0.5) ? j : (j - h);
in[i*h+j][0] = exp(-(x*x+y*y)/(a*a));

This allows the input to cover the entire Gaussian instead of a quarter of it. The entire code is attached below.

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

int main(int argc, char** argv)
{
    fftw_complex *in, *out;
    fftw_plan p;
    int i, j, w = 16, h = 16;
    double a = 2, x, y;

    in = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * w * h);
    out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * w * h);
    for (i = 0; i < w; i++) {
        x = (i < w*0.5) ? i : (i - w);
        for (j = 0; j < h; j++) {
            y = (j < h*0.5) ? j : (j - h);
            in[i*h+j][0] = exp(-1.*(x*x+y*y)/(a*a));
            in[i*h+j][1] = 0;
        }
    }
    p = fftw_plan_dft_2d(w, h, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);
    for (i = 0; i < w; i++) {
        for (j = 0; j < h; j++) {
            printf("%4d %4d %+9.4f %+9.4f\n", i, j, out[i*h+j][0], out[i*h+j][1]);
        }
    }

    fftw_destroy_plan(p);
    fftw_cleanup();
    fftw_free(in);
    fftw_free(out);
    return 0;
}

Why does negative number show up when I do FFT on a Gaussian?

1 Answers1