# Order of MATLAB FFT frequencies

This wikibook states that the output of MATLAB's FFT corresponds with the wavenumbers ordered as:

$$k=\left\{0,1,...,\frac{n}{2},-\frac{n}{2}+1,-\frac{n}{2}+2,...,-1\right\}$$

However, in the example codes on the same page, the wavenumbers are coded as

k = [0:n/2-1 0 -n/2+1:-1];


which is the same as the first, but with the $n/2$-wavenumber (the "maximum wavemnumber") replaced with $0$. It seems strange that they would include $0$ twice.

It seems the correct order is necessary for taking derivatives via Fourier transforms, as described in the wikibook. Which of these is correct, and does MATLAB document this anywhere?

• look into the function fftshift. It will take the output of fft, and reorder it from [-n/2+1 : n/2-1], which should help with your confusion. Feb 18 '13 at 18:05
• Isn't this something that you can quickly test for yourself? Finding out in which order the output of FFT is seems easy to determine with a couple of experiments. Jan 13 '18 at 20:39

I want to expand on my comment and rework the example you reference in a way that should be more understandable than the original and to explain why fft returns the coefficients the way it does.

For reference, the fft portion of the example is :

Nx = size(x,2);
k = 2*pi/(b-a)*[0:Nx/2-1 0 -Nx/2+1:-1];
dFdx = ifft(1i*k.*fft(f));
d2Fdx2 = ifft(-k.^2.*fft(f));


I added another section of code directly below it:

Nx = size(x,2);
k = 2*pi/(b-a)*(-Nx/2:Nx/2-1);
dFdxp = ifft(ifftshift(1i*k.*fftshift(fft(f))));
d2Fdx2p = ifft(ifftshift(-k.^2.*fftshift(fft(f))));


and wrapped both pieces of code in a tic; toc for rought timings. In a more readable format, the second method uses:

ckf = fftshift(fft(f));
ckdf = 1i*k.*ckf;
df = ifft(ifftshift(ckdf));


The first difference is that the second example has a much more intuitive k. This is the main advantage of the second example, since k is now in the form that we think about them. In the second and third lines I had to add fftshift around the call to fft, then a call to ifftshift directly inside the call to ifft. These additional function calls reorder the coefficients from what is required for the computer to work with them to the way humans usually think about them.

The issue with the second example, is that while k is more intuitive for us, this leaves the internal matrices for solving and inverting fft in forms that aren't as advantageous. So either we have to switch the order with calls to fftswitch and ifftswitch or it has to be hard coded into the fft functions. This is less prone to error from users (assuming they are unfamiliar with the workings of fft, as many people are), but you pay a price in run time.

As I stated before, I added timing calls around the two blocks for comparison and ran for multiple N. The timing results were:

N =     1000,  Ex1 = 0.000222 s,   Ex2 = 0.007072 s
N =    10000,  Ex1 = 0.001576 s,   Ex2 = 0.003506 s
N =   100000,  Ex1 = 0.023857 s,   Ex2 = 0.034051 s
N =  1000000,  Ex1 = 0.213816 s,   Ex2 = 0.406250 s
N = 10000000,  Ex1 = 4.555143 s,   Ex2 = 7.102348 s


As you can see, the act of switching the values back and forth slows the process considerably, especially at low N (where it 30x slower). This is just an example, and your computer may show slightly different trends depending on things such as memory speed, processor cores/speed, etc. but it is illustrative of the point. The reason fft has confusing output is because it is saving you a nontrivial fraction of your computing time.

• This still doesn't answer the original question -- why does it work at all to include zero twice? And why does the documentation suggest that zero should only be included once? Jan 28 '15 at 17:05

Your question concerning wavenumber 'replacement' is rather tricky. In general, wavenumber modification of this sort is not intended to save flops, as some have suggested here, but instead designed to respect the analytic peculiarities of, say, certain differential operators. I'm very surprised that I couldn't find a related discussion in Trefethen's Spectral Methods. To proceed, I will assume that you are concerned about FFT-based spectral differentiation and that you are performing transformations on a domain of even cardinality.

The rule of thumb, for odd derivatives, is to set

k = [0:n/2-1 0 -n/2+1:-1];


and for even derivatives, to set

k = [0:n/2 -n/2+1:-1];


If you are dealing with any padding, the treatment of the Nyquist frequency is similarly tedious. There is an excellent write-up that proves and remarks on these topics here!

It seems that both are correct; here's why. The $n/2$ mode is aliased to the zero mode on an equispaced grid with $n$ points (because that mode passes through zero precisely at all of the grid points). Therefore the two are indistinguishable. I'm not sure why in the code (and in Trefethen's book) they bother to replace $n/2$ with $0$; perhaps because a few flops can be saved in that way.

For a short explanation and demo of aliasing, see my IPython notebook here.

If $N$ is even (which we typically use for the fft to be fast), why set $N/2$ mode to zero? Well if the signal is real $x \in \mathbb R$, taking the inverse fourier transform, the mode $N/2$ may produce an imaginary part since there is no conjugate $-N/2$ component to cancel the $N/2$ mode contribution, so we set it to zero.