FFTW and fft with MatLab

Question

I have a weird problem with the discrete fft. I know that the Fourier Transform of a Gauss function exp(-x^2/2) is again the same Gauss function exp(-k^2/2). I tried to test that with some simple code in MatLab and FFTW but I get strange results.

First, the imaginary part of the result is not zero (in MatLab) as it should be.

Second, the absolute value of the real part is a Gauss curve but without the absolute value half of the modes have a negative coefficient. More precisely, every second mode has a coefficient that is the negative of that what it should be.

Third, the peak of the resulting Gauss curve (after taking the absolute value of the real part) is not at one but much higher. Its height is proportional to the number of points on the x-axis. However, the proportionality factor is not 1 but nearly 1/20.

Could anyone explain me what I am doing wrong?

Here is the MatLab code that I used:

    function [nooutput,M] = fourier_test

    Nx = 512;      % number of points in x direction

    Lx = 50;        % width of the window containing the Gauss curve

    x = linspace(-Lx/2,Lx/2,Nx);     % creating an equidistant grid on the x-axis

    input_1d = exp(-x.^2/2);                 % Gauss function as an input
    input_1d_hat = fft(input_1d);            % computing the discrete FFT
    input_1d_hat = fftshift(input_1d_hat);   % ordering the modes such that the peak is centred

    plot(real(input_1d_hat), '-')
    hold on
    plot(imag(input_1d_hat), 'r-')

You may be mixing up the continuous FT with the discrete FT. The continuous FT of a gaussian is a gaussian, but I'm not sure that's true for a DFT (FFT). — Paul R
Of course, discreteness make things imprecise but the discrete FT is created as an approximation of the continuous FT. When the number of points and the interval are large enough, which is the case here, the approximation should be good. There will be some deviations, but they should definitely not be so large and so systematic. Also, when I increase Nx and Lx, the result does not change which means that convergence is achieved. — Fortran_user
I think also there may be an implicit time shift in your DFT case, which translates to a phase rotation in the frequency domain - I would expect the magnitude to be correct and the phase to be a straight line (modulo 2π if you're not unwrapping it, i.e. a sawtooth). — Paul R

Buck Thorn Buck Thorn · Accepted Answer · 2013-10-10T14:46:59

The answer is basically what Paul R suggests in his second comment, you introduce a phase shift (linearly dependent on the frequency) because the center of the Gaussian described by input_1d_hat is effectively at k>0, where k+1 is the index into input_1d_hat. Instead if you center your data (such that input_1d_hat(1) corresponds to the center) as follows you get a phase-corrected Gaussian in the frequency domain:

    Nx = 512;      % number of points in x direction
    Lx = 50;        % width of the window containing the Gauss curve

    x = linspace(-Lx/2,Lx/2,Nx);     % creating an equidistant grid on the x-axis

    %%%%%%%%%%%%%%%%
    x=fftshift(x);   % <-- center
    %%%%%%%%%%%%%%%%

    input_1d = exp(-x.^2/2);                 % Gauss function as an input
    input_1d_hat = fft(input_1d);            % computing the discrete FFT
    input_1d_hat = fftshift(input_1d_hat);   % ordering the modes such that the peak is centered

    plot(real(input_1d_hat), '-')
    hold on
    plot(imag(input_1d_hat), 'r-')

From the definition of the DFT, if the Gaussian is not centered such that maximum occurs at k=0, you will see a phase twist. The effect off fftshift is to perform a circular shift or swapping of left and right sides of the dataset, which is equivalent to shifting the center of the peak to k=0.

As for the amplitude scaling, that is an issue with the definition of the DFT implemented in Matlab. From the documentation for the FFT:

For length N input vector x, the DFT is a length N vector X,
with elements
                 N
   X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                n=1
The inverse DFT (computed by IFFT) is given by
                 N
   x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
                k=1

Note that in the forward step the summation is not normalized by N. Therefore if you increase the number of points Nx in the summation while keeping the width Lx of the Gaussian function constant you will increase X(k) proportionately.

As for signal leaking into the imaginary frequency dimension, that is due to the discrete form of the DFT, which results in truncation and other effects, as noted again by Paul R. If you reduce Lx while keeping Nx constant, you should see a reduction in the amount of signal in the imaginary dimension relative to the real dimension (compare the spectra while keeping peak intensities in the real dimension equal).

You'll find additional answers to similar questions here and here.

FFTW and fft with MatLab

1 Answers