Matlab: for even real functions, FFT complex result, IFFT real result

Question

I am testing the validity of the FFT and IFFT functions in Matlab.

I can compare the outputs of these functions to a well-known mathematical fact: the Fourier transform of an even, real function (like a Gaussian centered at 0), is another even, real function (FFT[real, 0-centered Gaussian] = real, 0-centered Gaussian). This fact should hold both for FFT and IFFT.

First I make my grid:

nx = 256; % grid total pixel count
X = 500; % grid size (um)
dx = X/nx; %  grid spacing (um)
x = linspace(-nx/2,nx/2-1,nx)*dx; % x grid (um)

df = 1/(nx*dx); % spectral grid spacing (1/um)
f = linspace(-nx/2,nx/2-1,nx)*df; % f grid (1/um)

And I make my Gaussian:

A = 1; % magnitude (arbitrary units) 
x_fwhm = 7; % Full width at half maximum diameter (um)

x0 = x_fwhm/sqrt(2*log(2)); % 1/e^2 radius (um)
y = A*exp(-2*x.^2./(x0)^2); % Gaussian (arbitrary units) 

And apply the Fourier transform, with FFT:

yFFT = fftshift(fft(fftshift(y))); 

Alternatively, with IFFT:

yIFFT = fftshift(ifft(fftshift(y))); 

Plotting the results:

The IFFT does a perfect job: yIFFT is a purely real Gaussian. However, FFT yields a complex number: a very small imaginary part exists. This is fine, since an error should be expected in the fourier transform algorithms, and it's negligible anyway. What confuses me is why there is no error at all in IFFT? Are the FFT and IFFT algorithms that much different?

*** Note: fftshift and ifftshift are equivalent here, since my array has an even number of elements.

Answer 1

Handling of real-valued time-domain signal is a fairly common occurrence. So much so that the ifft function has built-in handling for the corresponding symmetry that arises in the frequency domain, as described under the "Algorithm" section of the documentation:

The ifft function tests whether the vectors in Y are conjugate symmetric. A vector v is conjugate symmetric when the i^th element satisfies v(i) = conj(v([1,end:-1:2])). If the vectors in Y are conjugate symmetric, then the inverse transform computation is faster and the output is real.

In other words, ifft constructs the imaginary part of yIFFT to be exactly 0 since it detects that your input has conjugate symmetry.

On the other hand, even time-domain signals are relatively less common and Mathworks didn't deem it necessary to perform a similar test in the fft function. That said, you could still take advantage of the conjugate symmetry test by computing the FFT using the ifft function through

% compute fft(x,[],dim) using ifft:
size(x,dim) * conj(ifft(conj(x),[],dim))

Answer 2

There is an error in your code:

yFFT = fftshift(fft(fftshift(y)));

should read

yFFT = fftshift(fft(ifftshift(y)));

The ifftshift function shifts the origin from the middle to the left-most bin. The fftshift shifts the origin from the left-most bin to the middle. These two operations are similar but not identical for even-sized odd-sized arrays. Note that your inverse transform suffers from the same issue, but shows no error as explained by SleuthEye. [never mind that bit, I think I was confused for a bit there, there should be no difference in your case.]

Also,

linspace(-nx/2,nx/2-1,nx)

would be better written as

-nx/2 : nx/2-1

given that nx is even. It's shorter, but it could also have a smaller numerical error due to the way things are computed (probably there's no difference, but I try to avoid linspace when the step size is 1.