Почему отрицательное число отображается, когда я делаю БПФ на гауссиане? - PullRequest
Новая
       5

Почему отрицательное число отображается, когда я делаю БПФ на гауссиане?

2 голосов
/ 25 июня 2010

Я использую эту библиотеку fftw .

В настоящее время я пытаюсь построить 2D гауссиан в форме 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]

1 Ответ

3 голосов
/ 26 февраля 2014

Вы забыли сместить центр гауссиана в середину (w / 2, h / 2) 

in[i*h+j][0] = exp(-(i*i+j*j)/(a*a));

следует читать 

in[i*h+j][0] = exp(-1.*((i-w/2)*(i-w/2)+(j-h/2)*(j-h/2))/(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;
    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(-1.*((i-w/2)*(i-w/2)+(j-h/2)*(j-h/2))/(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 i\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;
}
...