为什么在高斯上进行 FFT 时会出现负数？答案

【问题标题】：Why does negative number show up when I do FFT on a Gaussian?为什么在高斯上进行 FFT 时会出现负数？
【发布时间】：2010-06-24 22:44:28
【问题描述】：

我正在使用这个fftw 库。

目前我正在尝试以 e^(-(x^2+y^2)/a^2) 的形式绘制二维高斯曲线。

代码如下：

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

原来出现了负数。我认为高斯的傅里叶变换是另一个高斯，它不应该包含任何负数。

另外，当前原点在 [0]

【问题讨论】：

两个可能的问题：(i) 你正在做一个离散 FT，而不是 FT 和 (ii) 你的高斯函数已经有效地乘以一个矩形窗口（因为它被截断），相当于通过窗函数的FT在频域进行卷积。
(iii) 浮点数学是出了名的不精确，至少从幼稚的标准来看是这样。 -1E-30 将是一个舍入错误，而不是一个真正的负数。
我猜这就是原因。我当然不能查，因为我不能在我的电脑上做连续的FT。

标签： c++ fftw fft

【解决方案1】：

编辑：先前的答案是错误的，移动高斯中心无济于事，因为它引入了另一个相移。正确的解决方案是将高指数包装成负指数：

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

这允许输入覆盖整个高斯而不是四分之一。完整的代码附在下面。

#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;
}

【讨论】：