1
$\begingroup$

I am using the Numpy fft2, ifft2, and related functions and I am sometimes running into a strange situation where the output after performing the inverse Fourier transform the results are shifted. I know about fftshift and ifftshift, however my understanding of those is that they are used after fft2 and before ifft2 if it is needed for the DC component to be in the middle of the result instead of at one end. In particular, they are not supposed to be needed in real-space.

One example is in a demonstration of the convolution theorem. The following two things should be identical:

$$F^{-1}(F(im) \cdot circle)$$ $$im \otimes F^{-1}(circle)$$

where $F(x)$ is the Fourier transform, $F^{-1}(x)$ is the inverse transform, $\cdot$ is pointwise multiplication and $\otimes$ is convolution, $im$ is the image, and $circle$ is a simple, ideal, low-pass filter. The code I use to compare these is as follows:

# The fourier duck is a 256x256 8-bit grayscale image
im = skimage.io.imread('fourier_duck.png') / 255

# Create the ideal low-pass filter
x,y = numpy.meshgrid(range(-128,128), range(-128,128))
circle  = numpy.fft.ifftshift(x*x + y*y < 15*15)

# Perform filtering in Fourier space
f_im = numpy.fft.fft2(im)
f_im_circle = f_im * circle
finv_f_im_circle = numpy.real_if_close(numpy.fft.ifft2(f_im_circle), 10000)

# Perform filtering in real space
finv_circle = numpy.real_if_close(numpy.fft.ifft2(circle))
conv_im_finv_circle = scipy.ndimage.convolve(im, finv_circle, mode='wrap') # WARNING: takes lots of memory

The two outputs, finv_f_im_circle and conv_im_finv_circle should be quite similar but are quite different:

Filtering in Fourier spaceFiltering in real space

If the final, real-space, version is put through numpy.fft.fftshift it looks how I would expect (even though the plots don't show it, the value ranges are essentially the same as well, within FP error). The problem seems to crop up with the inverse Fourier transform of the circle which produces a kernel that also looks shifted.

I have tried doing shifts in other places and looking at the imaginary parts all to no available.

Besides the simple example I show above this has also occured with other filters designed in Fourier space, sometimes when doing inverse filtering, and even with filters designed in real space but applied in Fourier space. I can't figure out the pattern of when to shift the final result or when not to.

Ultimately my question is two fold:

  1. Why is this happening?
  2. How do I fix it? (which may just be how do I predict when I need to use the extra real-space numpy.fft.fftshift)

For reference the original duck image is available here

$\endgroup$
1
$\begingroup$

If you have a filter in the time/space domain, the location that is associated with zero linear phase in the frequency domain is the top left pixel of the filter. If the filter is centered in time/space, then you need to ifftshift prior to taking the FFT in order to remove the linear phase that centers it in the time domain. The following example Octave code (sorry, I don't use python) illustrates this:

x = zeros(10,10);   % Make an impulse in time, so that its 
x(6,6) = 1.0;       % magnitude in frequency is all ones.

% Display the angle of the FFT both with
% and without the time-domain quad-swapping.
figure;
imagesc( angle( fft2(x) ) );
title( 'Phase Without Quad-Swapping In Time' );

figure;
imagesc( angle( fft2( ifftshift(x) ) ) );
title( 'Phase With Quad-Swapping In Time' );

The same thing is true in the frequency domain, which is what is happening with your filter. Specifically, you have defined the magnitude of the filter in the frequency domain, but you have not given it the phase necessary to center it in the time/space domain.

I tried to answer a related question here.

$\endgroup$
  • $\begingroup$ Thanks for the explanation. A general fix (at least for various situations I have run into) was to reproduce the functions psf2otf and otf2psf from MATLAB and Octave. This shifts the middle of the PSF and unshifts during FFT and it appropriately unshifts it during inverse. $\endgroup$ – thaimin Feb 24 '18 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.