GSL Fast-Fourier Transform - Double-Valued Gaussia

2019-09-06 12:46发布

All I want to do is evaluate the Fourier transform of a Gaussian (as a test of course), but I'm getting what looks like a double-valued function, as can be seen in the plot below. When you look closely at the tails of the Gaussian, you can see that the FFT returned values are actually oscillating between positive and negative, making the DFT look "double-valued". Why is this happening? Is there a way to fix this (without just plotting the complex modulus of course)?

#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_errno.h>
#include <fstream>
#include <iostream>
#include <iomanip> 

#define REAL(z,i) ((z)[2*(i)]) //complex arrays stored as    [Re(z0),Im(z0),Re(z1),Im(z1),...]
#define IMAG(z,i) ((z)[2*(i)+1])
#define MODU(z,i) ((z)[2*(i)])*((z)[2*(i)])+((z)[2*(i)+1])*((z)[2*(i)+1])
#define PI 3.14159265359

using namespace std;

int main(){

    int n = pow(2,9);
    double data[2*n];
    double N = (double) n;

    ofstream file_out("out.txt");

    double xmin=-10.;
    double xmax=10.;
    double dx=(xmax-xmin)/N;
    double x=xmin;

    for (int i=0; i<n; ++i){
        REAL(data,i)=exp(-100.*x*x);
        IMAG(data,i)=0.;
        x+=dx;
    }

    gsl_fft_complex_radix2_forward(data, 1, n); 

    for (int i=0; i<n; ++i){
        file_out<<(i-n/2)<<"    "<<REAL(data,((i+n/2)%n))<<'\n';
    }

    file_out.close();
}

enter image description here

标签: fft gsl dft
1条回答
唯我独甜
2楼-- · 2019-09-06 13:31

You are plotting just the real component, e.g. the even or cosine component. Note that a cosine wave of integer frequency N toggles between being -1 in the middle and to being 1 in the middle, as N increases from an odd number to an even number. Thus any noise in the middle of the input to a DFT window can cause various real components in the DFT result to toggle (unless that noise is exactly orthogonal to all those DFT basis vectors) .

查看更多
登录 后发表回答