Correct order of implementing fftshift and ifftshift (in python)

Question

I want to Fourier transform a function psi(x), multiply it by a k-space function exp(-kx^2-ky^2), and then inverse Fourier transform the product back to x-space.

But my x-space and k-space grids are centred, and I know that I need fftshift and ifftshift to implement my k-space multiplication properly. But I don't understand how they work, so I don't know in which order to implement them. Could someone please tell me if I have done it correctly here?

import scipy.fftpack as spfft
import numpy as np

#Create a centred k-space grid]

kxmax, kymax = 10,10
kxgrid = np.linspace(-kxmax/2, kxmax/2, NX)
kygrid = np.linspace(-kymax/2, kymax/2, NY)
KX, KY = np.meshgrid(kxgrid, kygrid, indexing='xy')

psi = spfft.ifft2(spfft.fftshift(np.exp(-(KX**2 + KY**2)) * spfft.fftshift(spfft.fft2(psi))))

Answer 1

The need for shifting either the signal or the spectrum arises from constraints imposed by DFT/IDFT. These functions have expectations on a reference in their input, the sample corresponding to time t=0 (for the DFT), and to frequency f=0 (for the IDFT). Both functions also output their result using the same convention. Just before presenting shift / inverse shift functions, let's cover what DFT/IDFT expects.

DFT constraints

Periodicity constraint

The input sequence is assumed to be a complete period of the signal and can be juxtaposed to create a periodic signal, e.g. with a 100-sample:

Any 100-sample section of this signal could be used to take the DFT. Example with sections starting at magenta marks above:

The DFT magnitude doesn't change, but the time difference is taken into account in the phase (last row). For the DFT a delay in the signal results in the corresponding spectrum being multiplied by a complex exponential exp(j.2.π.f0.t). As the magnitude is 1, the product only adds 2.π.f0 to the spectral coefficient angles, hence the phase difference.

Ordering constraint

It means DFT evaluates the delay: It assumes the first sample is used for t=0, and computes phases relative to it. The output of the DFT is itself periodic with period N, N being the number of samples in a period:

(Only the magnitude is shown here, but the phase has the same periodicity.)

IDFT constraints

Like previously we could take any 100-sample section of the previous sequence and take an IDFT:

IDFT has the same input constraints than DFT: The input spectrum is assumed to be periodic, the first sample is for f=0 (DC). Any frequency shift in the spectrum corresponds to a multiplication of the signal by a complex exponential with magnitude 1. As for the DFT, this only changes the phase of the output.

Closer look to the horizontal scales

The time corresponding to sample k (k between 0 and N-1) is t=k/N.Ts, Ts being the sampling period, and frequency sample k corresponds to frequency k/N.Fs, fs being the sampling frequency.

Therefore there are no negative frequencies in the output of the DFT, and no negative frequencies are expected in the input of the IDFT.

However, due to the inherent periodicity of DFT and IDFT, we are allowed to translate the spectrum from range 0 to fs to range -1/2fs to +1/2fs. Now we have a spectrum with half of the frequencies negative, but still located after the positive frequencies in the sequence.

Function fftfreq can be used to create the scale with positive and negative frequencies. This function has no direct connection with DFT/IDFT, this is just a helper which knows how to determine the required frequencies values.

The time range can be shifted from 0 to Ts to -1/2Ts to +1/2Ts the same way due to periodicity N.Ts.

Use of fftshift and ifftshift

Now let's answer your question about shifting.

Scales with with half positive and half negative indices are convenient for plotting symmetrical views. But there is a slight problem, negative indices are placed after 0 and positive indices.

fftshift and ifftshift can be used to reorder elements: fftshift prepares the sequence for plotting purpose, ifftshift restores the native order used/expected by DFT/IDFT and described in the first part. Note these functions perform no other action than reordering elements, they are not directly related to FT in spite of their names contains 'fft', and prefix 'i' has no connection with IDFT.

This reordering operation can be seen as a rotation within the sequence, or as a linear shift along the periodic sequence (signal or spectrum).

One can be surprised there are two functions to do the same change. Sequences with an even length have the same number of positive (DC/0 included) and “negative” frequencies. For odd-length sequences, there is one additional negative frequency. Both “half” are not the same length, this must be taken into account. So fftshift manage the size this way:

for k in axes:
    n = tmp.shape[k]
    p2 = (n+1)//2
    mylist = concatenate((arange(p2, n), arange(p2)))
    y = take(y, mylist, k)
return 

Note how the location of the negative half is computed: p2 = (n+1)//2. And ifftshift is identical except it uses this location: p2 = n-(n+1)//2.

As only the user knows whether the sequence has been already swapped, he/she is responsible for using the appropriate function.

---

Bottom line: There is usually no need to shift or unshift sequences. DFT and IDFT expect and output sequences with index 0 as the first sample. A possible case is when a sequence is created centered, and has to be processed by DFT/IDFT. It must be unshifted with ifftshift before being used as input.

Shifting (fftshift) maybe required to center the sample for 0, either the DC frequency or time t=0, for symmetrical display purpose. If a shifted sequence must be later reused with DFT/IDFT functions, then you can restore the native order (ifftshift). Track carefully whether a sequence is currently shifted or not, as this determines which function to use.

Answer 2

No you haven’t, but that’s ok, it can be very confusing.

First thing: fft and ifft require the origin to be in the beginning of the vector (or in your 2D case, in the top-left of the array). Is the input psi’s origin centered like KX? If so, its origin must be shifted to the beginning with ifftshift. (If not, then just leave it alone.)

Second: since KX and KY have origins in their centers, you have to unshift them: you need spfft.ifftshift(np.exp(-(KX**2 + KY**2)) (note the i).

Finally: your output psi will therefore have its origin in the beginning. If you want its origin to be centered like KX, fftshift it.

In summary:

inputOriginStart = # ...
inputOriginStartFFT = spfft.fft2(psiOriginStart)
filterOriginStartFFT = spfft.ifftshift(np.exp(-(KX**2 + KY**2)))
outputOriginStart = spfft.ifft2(filterOriginStartFFT * inputOriginStartFFT)

where inputOriginStart is the input psi assuming it’s origin is in the beginning, and where outputOriginStart is the output psi—renamed for clarity. (I always go for clarity. If it does not work, you can more easily figure it out.)

Edit fixed the error pointed out by the asker—yes, I had a mistake, leave psiOriginStart’s origin at the start; then ifftshift the centered-origin function of KX and KY. (If you want to unshift outputOriginStart’s origin to the center then use fftshift.)

Edit 2 separated the filter (function of KX and KY) from the data to make the correct parentheses obvious.

---

How to keep these straight? A few tricks to remember:

• fft and ifft always need inputs and give outputs whose origins are in the beginning. This should be easy to remember from experience.
• fftshift takes the beginning-origin that fft needs/makes and shifts origin to the center. Again, I tend to readily remember this because of muscle-memory from typing fftshift(fft(...)) a thousand times.
• Finally, the only remaining thing is to deduce that ifftshift is the inverse of fftshift: it takes centered-origin vectors/arrays and shifts the origin to the beginning.