Hi David,
I believe the reason for this is the following. Staying in 1d for simplicity and not worrying about constants, the fourier transform of e.g. a gaussian is
exp(-x^2/2) --> exp(-k^2/2)
and the translated version picks up a phase factor.
exp(-(x-a)^2/2) --> exp(-ika) exp(-x^2/2)
The issue is that for the fft, x=0 is the first point in the configuration space array. If the span of x is from 0 to x0, and you create a blob in the middle of the array, to the fft it basically looks like blob(x-x0/2) and picks up a highly oscillatory phase factor of exp(-ik x0/2). The ifft can of course undo that, but in the meantime you are interpolating, and the phase factor plays havoc with the interpolation.
The extra fftshift takes the center of the blob down to x=0 and results in a good transform with moderate values of k.
If the number of fft points is even it doesn't matter, but if the number is odd, then I think it's necessary to use
FT_shape = fftshift(fftn(ifftshift(shape)));
distorted_shape = fftshift(ifftn(ifftshift(distorted_FT_shape)));
That's because for odd n, fftshift and ifftshift are different. for odd n, fftshift is not its own inverse, same for ifftshift. But fftshift and ifftshift are inverses.
>> fftshift(1:9) ans = 6 7 8 9 1 2 3 4 5
>> ifftshift(1:9) ans = 5 6 7 8 9 1 2 3 4
Since you are here and conversant in this whole topic I have a question for you. In time and frequency, everything has a name.
exp(2 pi i f t) exp(i w t) w = 2 pi f
In config space and reciprocal space,
exp(2 pi i [?] x) exp(i k x) k = 2 pi [?]
Is there a standard name and symbol for [?] = 1/lambda, other than the cumbersome 'inverse wavelengh'? Spectroscopists use cm^-1 but that's unit-dependent.
Best Answer