fftn, fftshift and interp3

In short, I'm trying to take the FT (FT_shape) of a 3D binary shape (shape), use interp3 on the resulting complex array (distorted_FT_shape), and apply an inverse FT to recover a "distorted" version of the shape (distorted_shape).

This, in theory, should work:

FT_shape = fftshift(fftn(shape));
distorted_FT_shape = interp3(Xgrid, Ygrid, Zgrid, FT_shape, rel_x, rel_y, rel_z,'linear',0);
distorted_shape = ifftn(ifftshift(distorted_FT_shape));

However, distorted_shape doesn't give me a "shape" when plotted. It appears I need an additional fftshift and ifftshift when taking the ffn and ifftn as follows:

FT_shape = fftshift(fftn(fftshift(shape)));
distorted_FT_shape = interp3(Xgrid, Ygrid, Zgrid, FT_shape, rel_x, rel_y, rel_z,'linear',0);
distorted_shape = ifftshift(ifftn(ifftshift(distorted_FT_shape)));

The second set of lines gives me a real, distorted shape as expected, but the additional fftshift and ifftshift do not make sense according to the MATLAB documentation. Can someone please explain why?

Additional info:

shape: a 256x256x256 zero array with a "blob" of ones size 128x128x128 about the centre of the array

Xgrid, Ygrid, Zgrid: meshgrid(-(256-1)/2:(256-1)/2, -(256-1)/2:(256-1)/2, -(256-1)/2:(256-1)/2)

rel_x, rel_y, rel_: rotated Xgrid, Ygrid, Zgrid

I'm doing X-ray diffraction research, so I need to "distort" things in the reciprocal space, for those wondering

Answer 1

[David Goodmanson]

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.

[David Yang]

Cheers for the response David. It's much more clear now.

nu = 1/lambda is called the wavenumber. Normally I see it expressed as k=2pi/lambda (the angular wavenumber or magnitude of the wave vector)

[David Goodmanson]

Hi David,

It would be more accurate I suppose to call k the angular wave number and 1/lambda the wave number, in parallel to f as frequency and w as angular frequency. But people just refer to w as 'omega'. And I have always heard k referred to as the wave number, hardly ever as the angular wave number. It doesn't really work to refer to both k and 1/lambda as the wave number. We seem to lack a symbol for 1/lambda, which would make it possible to remove ambiguity as with f and w.