fftshift with even and odd number of data points (scaling the positive and negative frequency axis)

If we wish to draw double sided frequency axis, i.e., positive and negative frequencies with odd and even number of data points respectively, how should we scale the frequency axis in each case? There is perhaps a very subtle point here.

Fs= sampling frequency

1) In the case of an even number of data points (N), the maximum values on the frequency axis are (Fs/2). However, there is one extra value in the output corresponding to the "negative" frequency when N is even as shown in this table after fftshift. Apparently, the extreme FT values are not the same after fftshift (as shown in the Table) when N is even.

2) In the case of an odd number of data points (N), the maximum values on the frequency axis are not (Fs/2), they are always less than (Fs/2). However, FT output are symmetric after fftshift.

What is the proper way to address this scaling issue in MATLAB for a double sided spectrum? Let us say we have the following, with Fs= 20 Hz

A_even=[0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0]

length(A_even)

B=fft(A_even)

C=fftshift(B)

X_odd=[0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 ]

length(X_odd)

Y=fft(X_odd)

Z=fftshift(Y)

3) This Table from an online FFT Primer summarizes the above query. It does not explain why MATLAB prefers a negative frequency value first in case of even number of data points. Please see the case of N=2, 4, 6.

Thanks.

Answer 1

[David Goodmanson]

Hi Farooq,

[1] For N odd, the frequencies before and after using fftshift (or ifftshift) are as you say:

N = 9 example:

for frequencies with spacing deltaf,

frequency array corresponding to the fft output array is
[0 1 2 3 4 -4 -3 -2 -1]*deltaf

if fftshift is then employed,

frequency array corresponding to the shifted fft output array is
[-4 -3 -2 -1 0 1 2 3 4]*deltaf

Fairlly straigthforward.

You can use ifftshift to undo the effect of fftshift, but note that for odd N, fftshift and its inverse function ifftshift are different functions. Consequently fftshift is not its own inverse.

[2] For even N, aside from zero, pos and neg frequencies there is also the Nyquist frequency, corresponding to exactly half an oscillation in each time interval. That frequency is really no more negative than it is positive.

N = 10 example:

frequency array corresponding to the fft output array is
[0 1 2 3 4 ny -4 -3 -2 -1]*deltaf

if fftshift is then employed,

frequency array corresponding to the shifted fft output array is
[ny -4 -3 -2 -1 0 1 2 3 4]*deltaf

Using the 'ny' label removes the unsatisfying asymmetry in the indices that you saw in Table 1.

Here fftshift and ifftshift are identical and both swap the two halves of the array. For even N, fftshift is its own inverse.

[FW]

Thanks for the useful discussion. The point you raise about the N=10, is

frequency array corresponding to the shifted fft output array is

[ny -4 -3 -2 -1 0 1 2 3 4]*deltaf

[1] I think this is would be my question. With even N, the frequency axis extremes should be +/- Fs/2, at least on paper. However in the array we have only one value corresponding to the Nyquist frequency. Following the notation of the table in the first column, N=10, we have

[c0 c1 c2 c3 c4 c-5 c-4 c-3 c-2 c-1] which after fftshift should become

[c-5 c-4 c-3 c-2 c-1 c0 c1 c2 c3 c2 c1]

so is your ny = c-5?

So how is our double sided frequency spectrum symmetric for N=even?

[2] Is it documented anywhere that MATLAB, calculates the correspoding the negative frequency value first? i.e. c-5 is calculated first rather than c5. I could not find this discussed anywhere except in that Table.

Example:

A_even=[ 0 0 0 1 1 1 1 0 0 0]

N=10

B=fft(A_even)'

C=fftshift(B)

Then C is equal to

C =

ny= 0.0000 + 0.0000i

c-4= -0.3090 - 0.9511i

c-3=0.4271 + 0.5878i

c-2= 0.8090 + 0.5878i

c-1=-2.9271 - 0.9511i

c0= 4.0000 + 0.0000i

c1=-2.9271 + 0.9511i

c2=0.8090 - 0.5878i

c3= 0.4271 - 0.5878i

c4= -0.3090 + 0.9511i

This implies c-5 = ny = 0.0000 + 0.0000i ? and c0= 4.0000+0000i

Thanks.

[David Goodmanson]

Hi Farooq,

For even N, the double sided frequency spectrum isn't symmetric about zero. For the fft, the frequencies that correspond to the output array are

0:N-1

multiplied by delf = Fs/N which I will leave out for the moment.

You can split this into two halves,

0:(N/2)-1 and (N/2):N-1 

The Nyquist frequency is (N/2), at the lower end of the upper half, just as 0 frequency is at the lower end of the lower half. So there is a kind of symmetry there.

Because of aliasing you can subtract N off of the upper half, getting

-(N/2):-1

which takes those frequencies just below the first half. With this and the first half, altogether you have (putting delf =Fs/2 back in)

 ( -(N/2):(N/2)-1 ) * Fs/N  = -Fs/2 : Fs/2 -Fs/N

So the frequency array definitely does not go from -Fs/2 to Fs/2.

[FW]

Thank you David! This makes perfect sense. Is there a reason as to why MATLAB prefers a negative coefficient first e.g. why c-5 rather than c5 is calculated first?

Just to confirm, for the N=odd case, it should be

( -(N-1)/2):(N-1)/2) ) * Fs/N

[David Goodmanson]

Hi Farooq,

The odd case is as you say. As is mentioned, for N even If you splint the array in half, 0 frequency is at the lower end of the lower half and the nyquist frequency is at the lower end of the upper half.

For even N there is no real middle point and there are two possible choices to place 0 frequency with an fftshift.

a) highest point of lower half, point N/2

b) lowest point of upper half, point N/2+1.

simply swapping the two halves of the array accomplishes b) so that's what they did. You end up with a frequency array -(N/2) : (N/2)-1