Confusion about fftshift and fft2 in MATLAB

Question

I'm mighty confused about how to correctly apply the FFT transform and its inverse in Matlab. I have a program where I need to

1. apply FFT2 to a matrix of size 4x32 (corresponding to modes m=-1:2, n=-15:16)
2. do some processing, which leads to a coefficient matrix for another function, whose Fourier coefficients relate to the first set of data by a simple algebraic (component-wise) formula
3. Use some properties of the two functions that lets me calculate what I need by just summing up an expression on the form 2*abs(A_n)*cos(phi+n*theta+alpha_n) where A_n is the nth coefficient of the m=1 mode and alpha_n = arg(A_n).

I have experimented a little with the FFT2 function and tried to understand how it arranges its output. From what I understand (from sources in my course literature), the coefficients will be ordered as is illustrated by the following script:

>>m = -1:2; n = -7:8;
>>[N,M] = meshgrid(n,m); 
>>MN = M; MN(:,:,2) = N; 
>>asfft = @(X) [X(2:4,8:16,:) X(2:4,1:7,:); X(1,8:16,:) X(1,1:7,:)];
>>asfft(MN)
ans(:,:,1) =
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2
    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1
ans(:,:,2) =
     0     1     2     3     4     5     6     7     8    -7    -6    -5    -4    -3    -2    -1
     0     1     2     3     4     5     6     7     8    -7    -6    -5    -4    -3    -2    -1
     0     1     2     3     4     5     6     7     8    -7    -6    -5    -4    -3    -2    -1
     0     1     2     3     4     5     6     7     8    -7    -6    -5    -4    -3    -2    -1

where asfft reorders the indices in the same way I believe fft2 does, but does nothing else. In other words, each index is ordered from 0 up to the max, then from the min to -1. According to the documentation, I should be able to rearrange this so I get 0 in the middle by using fftshift, but it's not giving me the output I expect. Instead, I get this:

>> fftshift(asfft(MN))
ans(:,:,1) =
     8    -7    -6    -5    -4    -3    -2    -1     0     1     2     3     4     5     6     7
     8    -7    -6    -5    -4    -3    -2    -1     0     1     2     3     4     5     6     7
     8    -7    -6    -5    -4    -3    -2    -1     0     1     2     3     4     5     6     7
     8    -7    -6    -5    -4    -3    -2    -1     0     1     2     3     4     5     6     7
ans(:,:,2) =
     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2
    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1

As you can see, the max is on the wrong side of the spectrum - instead of -7 -6 ... -1 0 1 ... 8 and -1 0 1 2 I have 8 -7 6 ... and 2 -1 0 1. This is fatal for me, since to be able to do the calculation described in 3 above, I need to know the indices of the respective coefficients. (The two layers are also switched, but that doesn't matter for me since I'm only going to do this on MxN matrices, not on N-d-arrays, later.)

Why is this? What am I doing wrong here?

Answer 1

Considering the simple 1D case first, fft gives you:

[X(0) X(1) X(2) ... X(N/2-1) X(-N/2) X(-N/2+1) ... X(-1)]

where I'm using X(k) to denote elements of the mathematical DFT.

All fftshift does is rotate these by N/2, so you end up with:

[X(-N/2) X(-N/2+1) ... X(-1) X(0) X(1) X(2) ... X(N/2-1)]

i.e. linear order.¹

In the multi-dimensional case, fftshift simply applies this rotation in all dimensions.

---

_{1. Although since by definition, X(k) == X(N+k) for a DFT, the unrotated vector is also in "linear" order!}

Answer 2

Check this newsread entry too http://www.mathworks.com/matlabcentral/newsreader/view_thread/285244

Answer 3

taking a queue from Oli, the problem is that you're not parsing it in the middle:

so for asfft = @(X) [X(2:4,8:16,:) X(2:4,1:7,:); X(1,8:16,:) X(1,1:7,:)];

size(8:16)
ans =
     1     9
size(1:7)
ans =
     1     7

so do:

asfft = @(X) [X(3:4,9:16,:) X(3:4,1:8,:); X(1:2,9:16,:) X(1:2,1:8,:)];