Using FFT and rescaling properly on a Gaussian in MATLAB

0
votes

I have an issue getting the correct output for the Fourier transform of Gaussian exp(-2(x^2 + y^2). I know analytically that the Fourier transform should be 0.25 * exp(-0.125(k^2 + j^2)), where k and j are the Fourier variables, but the output of the absolute value of the FFT output doesn't match this. This is the analytic solution:

versus my output:

This is the relevant code:

Fs = 10;                                        %space frequency
Dx = 1 / Fs;                                    %sampling period
x = -10: Dx : 10;                               %space vector on domain [-10, 10], one dimensional
L = length(x);                                  %length of signal
[X, Y] = meshgrid(x, x);                        %2D domain defined by the space vector x

U =  exp(-2 *( X.^2 + Y.^2));                                      
n = 2^nextpow2(L);                              %padding for FFT to optimize performance

U_hat = fft2(U, n, n);                          %FFT of the Gaussian

Dk = Fs*((-n/2):((n/2) - 1)) / n;               %space frequency domain (ie, the fourier domain) in one dimension, shifted, rescaled by n
P = abs(fftshift(U_hat) / n);                   %power spectrum (ie | X_hat | = sqrt(X_hat * complex_conjugate(X_hat)), shifted
                                                                                                                 
[K1, K2] = meshgrid(Dk, Dk);                    %Fourier domain

mesh(K1, K2, P);
title('Gaussian Pulse in Frequency Domain');
xlabel('Frequency (k_1)');
ylabel('Frequency (k_2)');
zlabel('|P(f)|');

Did I mess up making the frequency domain?

2
matlabfftgaussian
You should try computing the analytical solution again. The numerical one is correct. – Cris Luengo
I just used Mathematica to get the solution symbolically, and it agrees with my claim for the analytic solution. Is there anything I’m missing or misinterpreting? Particularly, I used FourierTransform[Exp[-2(x^2 + y^2)], {x, y}, {k, l}]. I know the Fourier transform of a Gaussian is another Gaussian, which is why I didn’t question the result. – dweeby
That’s interesting, because the same company owns this website that says something different: mathworld.wolfram.com/FourierTransformGaussian.html – Cris Luengo
Ooh, that's weird. I'll compare my numerical answer to the solution on the link you sent me, which I think should be easy to compute analytically since I can just split up the 2D gaussian into two 1D Gaussians – dweeby
Thank you so much for all the help! I'm new to stackoverflow - how do I close this out? Should I follow up with "answer your question" saying you helped me? I don't think I can checkmark comments. – dweeby

2 Answers

1
votes

Credit to Cris Luengo for the answer in the comments - the analytic solution present in the question was incorrect. Mathematica reported the Fourier transform to be 0.25 * exp(-0.125(k^2 + j^2)) but the actual transform is (pi / 2) * exp(-( (pi^2) / 2) * (k^2 + j^2)). The code is correct as is.

0
votes

I analyzed your implementation for quite some time. Your frequency-domain implementation looks good. Why?

  1. You got gaussian like FFT O/P for gaussian like input. (More on this later)
  2. The frequency spectrum bin increments i.e df = Fs/n. Which says there is nothing wrong in the X-axis of your Frequency Transform.
  3. Your Y-axis of transform is abs of fftshift of fft divided by n, to normalize the transform. You did not multiply by 2 which most people would have taken as you are interested also in the negative frequencies.

So implementation is correct, but the final result is wrong... That can mean only one thing. Your gaussian signal only looks gaussian, but is actually not!

A gaussian wave should have the following equation.

Click for Gaussian signal equation

If your means are 0 and variance and sd are 1. It simply becomes

U = (1/(2pi))exp(-( X.^2 + Y.^2)/2) --- (1)

Even if you are to ignore the first part of (1) i.e 1/2*pi as it is a simple amplitude scaling, your equation *exp(-2 ( X.^2 + Y.^2)) still does not follow the basic equation of a gaussian wave.

Please recheck your input.