Thread Subject:

Correct use of fftshift and ifftshift at inpu tto fft and ifft

I've seen discussion here of the correct use of MATLAB's fftshift() function.

I'd like to verify my interpretation, and to ask what are the effects when people do not use it correctly.

Basically, MATLAB implements an fft (as do other sets of functions) whose results AND inputs are ordered (for an array of N elements), from 0 to (N/2-1) and then from –N/2 to -1. (I will call this 'swapped' from now on.) But in what I might call a 'natural' ordering, for a spectrum or time function centred at zero and extending equal time (or frequency) either side of zero, the arrays of input data and results would be ordered (2’s complement style) from –N/2 to N/2-1 where N is the number of elements. That is, time zero (or frequency zero) are at the centre of such an array.

It is well known that the results of MATLAB's fft() function must be re-ordered with the MATLAB function fftshift() so that they are 'naturally' ordered (for instance when plotting).

   H = fftshift( fft ( h ) );

Similarly, the inverse fft also requires this fftshift:

  h = fftshift ( ifft ( H ) );

It seems to be less well known that both the fft() and ifft() functions also expects their INPUTS to be swapped. You can see this if you do the inverse fft of an fft (yielding the original input):

  h == ifft( fft( h ) );

The fft() results are swapped, so clearly the ifft() expects its input to be swapped.

Thus, if you start with a naturally ordered spectrum (say, a desired filter frequency response centred at DC) then you must make it swapped before you use the ifft(). The function fftshift 'unswaps' an array, and the inverse function ifftshift() swaps one. So we use ifftchift BEFORE an ifft() in the case where the input (spectrum) is NOT the result of a MATLAB fft():

  h = ifft( ifftshift( H ) );

(where H is an unswapped array).

So to do an inverse FFT where both input and results are to be naturally ordered, using MATLAB, we must do:

  h = fftshift( ifft( ifftshift( H ) ) );

(which, by the way, most people seem not to do..).

The fft() also requires a similar idiom. You can see this if you consider the fft of an inverse fft (yielding the original input):

   H == fft( ifft( H ) );

Since the result of ifft() is swapped, clearly fft() must also expect its input to be swapped, so if we start with an unswapped array we must swap it:

  H = fftshift( fft( ifftshift( h ) ) );

(which almost nobody seems to do..).

(The reason things are like this can be explained in two ways, I think. First, that the original Cooley-Tukey FFT results were swapped and so all FFTs since return swapped results. Second, that the fft() implements a particular DFT equation that is a sum up from zero (no negative time or frequency) and so the inputs must be shifted from a DC-centred array.)

I verified this for myself by generating an input whose analytic transform I knew and could code, and comparing the results of the two idioms to the expected anlytic result:

   H = fftshift( fft( h ) ); % most common usage, wrong result
   H = fftshift( fft( ifftshift( h ) ) ); % uncommon usage, correct result

So far so good. My first questions:

1) is the above idiom for using MATLAB's fft() and ifft() correct?
2) is my explanation of why this is so, valid?

Now to the real question. What is the effect if people do not do it right?

The common way to use MATLAB's fft() is:

  H = fftshift( fft( h ) );

The lacking ifftshift is a circular shift by half the array length (there is a tweak depending if the array length is odd or even). If we look only at magnitude of the result, then the result is unchanged by any circular shift. And since probably most people do look only at the magnitude of the result of an FFT, perhaps this is why the incorrect idiom is most common? But if we look at complex numbers, or phase, then there is a clear difference and the most common idiom is wrong.

I came upon this while using MATLAB to design an equalization filter, from a DC-centred spectrum, by an FFT to get the filter coefficients (before windowing etc). The results were incorrect (at least, differd from my expectation) until I adopted the correct idiom I have outlined above.

So my next questions are:

3) are Mathworks aware of this and is it documented?
4) why do so many (most?) people not do it right?
5) what are the effects, and ramifications, of doing it wrong?

Chris
==================
Chris Bore
BORES Signal Processing
www.bores.com

On Jun 23, 5:10 am, "Chris Bore" <chris_b...@yahoo.co.uk> wrote:
> I've seen discussion here of the correct use of MATLAB's fftshift() function.
>
> I'd like to verify my interpretation, and to ask what are the effects when people do not use it correctly.
>
------SNIP
>
> So my next questions are:
>
> 3) are Mathworks aware of this and is it documented?

It's probably known.
It is documented about as poorly as one could imagine.
I cannot guess why ( Maybe it's too complicated?)

It gets more confusing when N is odd because
then fftshift and ifftshift yield different results
and cannot be used interchangeably.

> 4) why do so many (most?) people not do it right?

Many are not even aware.
Most are not interested in phase.
It is not easy to remember (especially when N is odd.)

> 5) what are the effects, and ramifications, of doing it wrong?

phase shifts of 2*pi*Fs/N

Hope this helps.

Greg

On Jun 23, 5:10 am, "Chris Bore" <chris_b...@yahoo.co.uk> wrote:
> I've seen discussion here of the correct use of MATLAB's fftshift() function.
>
> I'd like to verify my interpretation, and to ask what are the effects when people do not use it correctly.
>
> Basically, MATLAB implements an fft (as do other sets of functions) whose results AND inputs are ordered (for an array of N elements), from 0 to (N/2-1) and then from –N/2 to -1. (I will call this 'swapped' from now on.) But in what I might call a 'natural' ordering, for a spectrum or time function centred at zero and extending equal time (or frequency) either side of zero, the arrays of input data and results would be ordered (2’s complement style) from –N/2 to N/2-1 where N is the number of elements. That is, >time zero (or frequency zero) are at the centre of such an array.

You are assuming N is even. If so, neither the bipolar time waveform
xb
defined on dt*(-N/2 : N/2-1) nor the bipolar frequency waveform Xb
defined on df*(-N/2 : N/2-1) is "centered" at zero: there are only
N/2-1 positive points compared to N/2 negative points. The "center"
of the waveform is between 0 and the 1st negative point (-1/2?).

In addition, when N is even, ifftshift(:) = fftshift(:). Therfore,
it it is easy to fall into habits that are invalid when N is odd.

Therefore, if you make any rules about using the shift
functions, make sure they hold when N is odd.

Hope this helps.

Greg

Hello Greg,

Thanks for your reply (and earlier posts on a similar topic that I have read).

> You are assuming N is even. If so, neither the bipolar time waveform
> xb defined on dt*(-N/2 : N/2-1) nor the bipolar frequency waveform Xb
> defined on df*(-N/2 : N/2-1) is "centered" at zero: there are only
> N/2-1 positive points compared to N/2 negative points. The "center"
> of the waveform is between 0 and the 1st negative point (-1/2?).

Indeed. I think in a more generic way what is needed is a circular shift that results in the 'zero' frequency ending up at the 'zero' of the fft's assumed unipolar interval. For inputs that extend the same either side of zero the ifftshift works, but for other cases an explicit circular shift may be needed by the required amount.

In fact, this is all exploiting the property of the DFT in respect of circular shifts, to avoid an explicit demodulation (which is another way to look at it).

It is a pity, I think, because it is an example where MATLAB's abstract and generic way to express things becomes narrowed and stuck to the implementation of a particular function - and a particular way to implement that function too.

Greg Heath <heath@alumni.brown.edu> wrote in message <7ed0576b-2636-4017-95fa-0f9a452d4880@8g2000vbg.googlegroups.com>...
> On Jun 23, 5:10 am, "Chris Bore" <chris_b...@yahoo.co.uk> wrote:
> > I've seen discussion here of the correct use of MATLAB's fftshift() function.
> >
> > I'd like to verify my interpretation, and to ask what are the effects when people do not use it correctly.
> >
> > Basically, MATLAB implements an fft (as do other sets of functions) whose results AND inputs are ordered (for an array of N elements), from 0 to (N/2-1) and then from –N/2 to -1. (I will call this 'swapped' from now on.) But in what I might call a 'natural' ordering, for a spectrum or time function centred at zero and extending equal time (or frequency) either side of zero, the arrays of input data and results would be ordered (2’s complement style) from –N/2 to N/2-1 where N is the number of elements. That is, >time zero (or frequency zero) are at the centre of such an array.
>
> You are assuming N is even. If so, neither the bipolar time waveform
> xb
> defined on dt*(-N/2 : N/2-1) nor the bipolar frequency waveform Xb
> defined on df*(-N/2 : N/2-1) is "centered" at zero: there are only
> N/2-1 positive points compared to N/2 negative points. The "center"
> of the waveform is between 0 and the 1st negative point (-1/2?).
>
> In addition, when N is even, ifftshift(:) = fftshift(:). Therfore,
> it it is easy to fall into habits that are invalid when N is odd.
>
> Therefore, if you make any rules about using the shift
> functions, make sure they hold when N is odd.
>
> Hope this helps.
>
> Greg
>

On Jun 23, 1:29 pm, "Chris Bore" <chris_b...@yahoo.co.uk> wrote:
> Hello Greg,
>
> Thanks for your reply (and earlier posts on a similar topic that I have read).
>
> > You are assuming N is even. If so, neither the bipolar time waveform
> > xb defined on dt*(-N/2 : N/2-1) nor the bipolar frequency waveform Xb
> > defined on df*(-N/2 :  N/2-1) is "centered" at zero: there are only
> > N/2-1 positive points compared to N/2 negative points. The "center"
> > of the waveform is between 0 and the 1st negative point (-1/2?).
>
> Indeed. I think in a more generic way what is needed is a circular shift that results in the 'zero' frequency ending up at the 'zero' of the fft's assumed unipolar interval. For inputs that extend the same either side of zero the ifftshift works, but for other cases an explicit circular shift may be needed by the required amount.
>
> In fact, this is all exploiting the property of the DFT in respect of circular shifts, to avoid an explicit demodulation (which is another way to look at it).
>
> It is a pity, I think, because it is an example where MATLAB's abstract and generic way to express things becomes narrowed and stuck to the implementation of a particular function - and a particular way to implement that function too.
>
>
>
> Greg Heath <he...@alumni.brown.edu> wrote in message <7ed0576b-2636-4017-95fa-0f9a452d4...@8g2000vbg.googlegroups.com>...
> > On Jun 23, 5:10 am, "Chris Bore" <chris_b...@yahoo.co.uk> wrote:
> > > I've seen discussion here of the correct use of MATLAB's fftshift() function.
>
> > > I'd like to verify my interpretation, and to ask what are the effects when people do not use it correctly.
>
> > > Basically, MATLAB implements an fft (as do other sets of functions) whose results AND inputs are ordered (for an array of N elements), from 0 to (N/2-1) and then from –N/2 to -1. (I will call this 'swapped' from now on.) But in what I might call a 'natural' ordering, for a spectrum or time function centred at zero and extending equal time (or frequency) either side of zero, the arrays of input data and results would be ordered (2’s complement style) from –N/2 to N/2-1 where N is the number of elements. That is, >time zero (or frequency zero) are at the centre of such an array.
>
> > You are assuming N is even. If so, neither the bipolar time waveform
> > xb
> > defined on dt*(-N/2 : N/2-1) nor the bipolar frequency waveform Xb
> > defined on df*(-N/2 :  N/2-1) is "centered" at zero: there are only
> > N/2-1 positive points compared to N/2 negative points. The "center"
> > of the waveform is between 0 and the 1st negative point (-1/2?).
>
> > In addition, when N is even, ifftshift(:) = fftshift(:). Therfore,
> > it it is easy to fall into habits that are invalid when N is odd.
>
> > Therefore, if you make any rules about using the shift
> > functions, make sure they hold when N is odd.

After further thought I came up with this and just posted a
similar message on another thread ("Inverse Fourier Transforms"):


The MATLAB fft and ifft functions assume the nonnegative intervals
t = dt*(0:N-1) and f = df*(0:N-1) corresponding to the Fourier pair
x(t) and X(f).

If x and X are considered periodic, fftshift yields the translation
to
the "zero centered" bipolar intervals with indexing

tb = dt * [-ceil((N-1)/2) : floor((N-1)/2)]
fb = df * [-ceil((N-1)/2) : floor((N-1)/2)]

corresponding to the Fourier pair xb(tb) and Xb(fb).

ifftshift is the inverse of fftshift.

When N is even, fftshift and ifftshift yield
the same result.

However, when N is odd,

x = ifftshift(xb)
X = fft(x)
Xb = fftshift(X)

and to return

X = ifftshift(Xb)
x = ifft(X)
xb = fftshift(x)

The trick is to remember that fftshift will
"center" the waveform about the coordinate zero.

If in doubt type fftshift([-2:2]) into the command
line.

Hope this helps.

Greg

I always avoid to use fftshift and ifftshift as much as I can. If you pass a 2D image to the fft2 function, then the output 2D array always has the DC term at the first element position (i.e., the top-left corner of the 2D array for the row-column coordinate system). If all the frequency-domain filters that you design have the DC terms at the same position, then you don't need to use fftshift or ifftshift. I would use fftshift only if I want to visualize the Fourier spectrum of my image.

I have verified the following from my frequent uses of these functions:
* fftshift = ifftshift if your image has even numbers of rows and columns.
* always use fftshift if you want to move the DC term from the top-left corner to the centre of the 2D array, regardless of the size of your input image (odd or even number of rows and colums).
* always use ifftshift if you want to move the DC term from the centre of the 2D array to the top-left corner, regardless of the image size.
* label your frequency axes correctly. If the number of elements is odd (i.e. N is odd below), then use (assuming that your sampling rate is 1)
    w = linspace(-0.5, 0.5, N)
as your frequency array.
If the number of elements is even (i.e. N is even below), then use
   w = linspace(-0.5, 0.5, N); w(end) = [ ]
as your frequency array.

Very often we found that it is easier to directly design a filter in the frequency domain with the DC term in the centre of the array. If this is the case, then I would suggest that, after you have got your filter, you call ifftshift to swap the quadrants so that the DC term is at the top-left corner of the image. From then on, everything you operate on has the DC terms at the top-left corner. To move back to the spatial domain, you don't need to do any quadrant shifting, all you need to do is call ifft2.

The same applies to 1D signals -- just replace ifft2 by ifft and fft2 by fft in the above paragraphs.

Du

> H = fftshift( fft( h ) ); % most common usage, wrong result
> H = fftshift( fft( ifftshift( h ) ) ); % uncommon usage, correct result


Does which is right/wrong depend on whether h is centered at t=0 (namely the time is 0 ~ T-dt) or not (namely the time is -T/2 ~ T/2 )?

Nice discussion~