Question concerning fftshift in FFT operations.

[veritas9872]

Hello, I have a question concerning the ifftshift and fftshift operations in the current implementation.
I am aware that applying fftshift(fft(ifftshift(x))) for FFT and fftshift(ifft(ifftshift(X))) for iFFT is standard in MATLAB.
However, it is my understanding that the FFT implementation is different in different languages and I found that the FFT in NumPy (and hence PyTorch) does not require shifting. In fact, applying shifting produces the wrong answer.
Is there a reason for applying the shifting operations before and after the FFT operations? Possibly because the original k-space was produced with MATLAB compatibility in mind. Also, I think that having tests to check if the FFT operations are correct would be beneficial.

I include my test code below.

import numpy as np
from scipy import signal

def test_fft2_noshift():
    a = np.random.random((17, 17))
    b = np.random.random((17, 17))
    c1 = signal.convolve2d(a, b, mode='full')
    ar = np.fft.fft2(a, s=(34, 34))
    br = np.fft.fft2(b, s=(34, 34), norm='forward')
    c2 = np.fft.ifft2(ar * br, s=(34, 34), norm='forward')[:-1, :-1]
    assert np.allclose(c1, c2)

def test_fft2_shifted():
    a = np.random.random((17, 17))
    b = np.random.random((17, 17))
    c1 = signal.convolve2d(a, b, mode='full')
    ar = np.fft.fft2(np.fft.ifftshift(a), s=(34, 34))
    br = np.fft.fft2(np.fft.ifftshift(b), s=(34, 34), norm='forward')
    cr = np.fft.fftshift(ar * br)
    c2 = np.fft.ifft2(cr, s=(34, 34), norm='forward')[:-1, :-1]
    assert not np.allclose(c1, c2)

[mmuckley]

Hello @veritas9872, I think there should be some clarifications of what is going on with FFT. The FFT is well-defined to apply the Discrete Fourier Transform, which by convention always starts at 0 and goes to N-1. Effectively, this means the 0-phase is the first element and the middle element is the highest phase offset (or frequency in Fourier space). Conventionally with images, we don't usually think of the left side of the image as being the 0-phase position, so we use FFT shifts because they're more intuitive and they match our point-of-reference for how complex phase works across MRIs.

I think the reason you're getting an incorrect result here is due to the s=(34, 34), which applies zero-padding. Conventionally this might be done at the end of the sequence, and with the ifftshift operation before you are moving the location of the end. If you alter this to apply zero-padding at the correct location after ifftshift, you should get the right result.

Practically we don't use FFTs for convolutions like this in the repository. We use FFTs because they represent the physics of an MRI scanner, and the code in the repository is configured (and has been extensively tested) to reconstruct real-world MRI data.

[roflmaostc]

It just depends where you shift! Shifting before the FFT and ittshifting at the very end produces the same results since the convolution acts everywhere the same.

You always need to remember: fft Fourier transforms a signal and outputs the 0 frequency at the very first pixel. ifft does the inverse Fourier transform and expects the 0 frequency at the very first pixel in the input. fftshift shifts the 0 frequency (the first pixel) to the center pixel and ifftshift shifts the center pixel to the very first pixel.