CODEWITHSUNDEEP
matlabfftgaussian
dweeby
dweeby

Reputation: 11

Using FFT and rescaling properly on a Gaussian in MATLAB

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?

Upvotes: 0

Views: 536

Answers (2)

Pratik Purkayastha
Pratik Purkayastha

Reputation: 1

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.

Upvotes: 0

dweeby
dweeby

Reputation: 11

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.

Upvotes: 1

Related Questions