Use of MATLAB fftshift

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 (originators of MATLAB) aware of this and is it
documented? (I'll ask them..)
4) why do so many (most?) people not do it right?
5) what are the effects, and ramifications, of doing it wrong, when
MATLAB is so widely used?

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

On 23 Jun, 11:04, Chris Bore <chris.b...@gmail.com> wrote:

> So my next questions are:
>
> 3) are Mathworks (originators of MATLAB) aware of this and is it
> documented? (I'll ask them..)

You might want to be aware that matlab, in its early days,
had a number of amateurishly implemented signal processing
functions. These things might have improved with the redesign
of the signal processing toolbox of recent years, but don't
haev too high hopes.

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

Because few people are willing or able to make the effort
to learn the basics. There is no need for a FFTSHIFT function,
provided one knows what the FFT is and how it works.

> 5) what are the effects, and ramifications, of doing it wrong, when
> MATLAB is so widely used?

Read the current newspapers, about the Gulf oils spill, and
project to whatever subject you might be interested in.

Rune

On Jun 23, 5:04&#4294967295;am, Chris Bore <chris.b...@gmail.com> 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 &#4294967295;N/2 to -1. (I will call this
> 'swapped' from now on.)

and Chris, you're counting like a normal DSPer, not like MATLAB (who
counts only from 1 to N).

> 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&#4294967295;s complement style) from &#4294967295;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.
>
...
> 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 ifftshift
> BEFORE an ifft() in the case where the input (spectrum) is NOT the
> result of a MATLAB fft():
>
...
>
> &#4294967295; H = fftshift( fft( ifftshift( h ) ) );
>
> (which almost nobody seems to do..).

i do something like that.

recently i was involved in the discussion and had learned of a subtle
difference between fftshift() and ifftshift() (the latter i wasn't
even aware of) that never affected any of my code because i only have
been using the FFT with even N (usually a power of two).  now, setting
that difference aside, the main use for fftshift() (or maybe it should
be ifftshift()) is to swap the two halves when the center of your data
is where all the business is happening.  usually this happens (for me
in audio, at least) when you window off a piece of data and send the
windowed segment to the FFT.

if you look at the phase result of

   n0 = some_big_index_into_array_x;

   X = fft( x(n0+1:n0+N).*hanning(N)  );

you will see a phase term of pi added (or subtracted) to only every
odd (well, in stupid MATLAB it's every even) indexed value of arg(X)
(in MATLAB they would call it "angle(X)").  to appropriately center
the data passed to the FFT at bin 0 (well, in stupid MATLAB, that
would be bin 1), you have to send the first half (chronologically),
namely x(n0+1:n0+N/2), to the "negative time" half of the FFT input
and the second half, x(n0+1+N/2:n0+N), to the "positive time" half of
the input.  so

    X = fft( fftshift( x(n0+1:n0+N).*hanning(N) )  );

works a little better, if you want smoothly connected phase in the
spectrum of that snipped piece of array x().

to me, that's all this fftshift is about.

r b-j

On Jun 23, 5:04&#4294967295;am, Chris Bore <chris.b...@gmail.com> wrote:
> I've seen discussion here of the correct use of MATLAB's fftshift()
> function

Please do not send separate copies of the same post to different
newsgroups. Just post one copy but include all group names,
separated by commas, on the newsgroup line.


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

Greg Heath wrote:
> *SNIP*
>
> Please do not send separate copies of the same post to different
> newsgroups. Just post one copy but include all group names,
> separated by commas, on the newsgroup line.
>

Did you follow your own advice ?   *NOT* lol
Was it possible? ??? ????? [I'm using SeaMonkey 2.0.5 as reader]

May I suggest posting to all 'relevant' groups with "Reply To:" 
set to *ONLY* one group [that group noted in body]


I now decease from "rant"

On Jun 23, 11:30&#4294967295;pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 23, 5:04&#4294967295;am, Chris Bore <chris.b...@gmail.com> wrote:
>
> > I've seen discussion here of the correct use of MATLAB's fftshift()
> > function
>
> Please do not send separate copies of the same post to different
> newsgroups. Just post one copy but include all group names,
> separated by commas, on the newsgroup line.
>
> 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 &#4294967295;= ifftshift(xb)
> X &#4294967295;= fft(x)
> Xb = fftshift(X)
>
> and to return
>
> X &#4294967295;= ifftshift(Xb)
> x &#4294967295;= 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

Thanks to all contributors.

Can I state the issue another way, and try to suggest a more generic
idiom?

The issue, I think, is that the fft() and ifft() in MATLAB are defined
on an interval n = [0 : N-1] where N is the number of samples, and n
is the 'sample number' that has n ==0 at the 'center' (time or
frequency). Whereas, I addressed the problem of a signal or spectrum
that is defined on the interval n = [ -N/2 : (N/2 - 1) ].

First, because the FFT assumes a periodic function, we can use a
circular shift to translate from one interval to another.
Second, the ifftshift() does such a circular shift for the interval I
specified (with appropriate caveats as to correct use)

Second, the issue is partly that it is easy to confuse the matrix
index (call this 'i' where in matlab i = [1 : N ] ) and the 'sample
number' n (which are not the same).

So in matlab's native usage, n( 1 ) == 0 whereas in my interval n( N/2
+ 1 ) == 0. That is, to re-state the obvious , the 'centers' are at
index i == 1 and i == (N/2+1) respectively.

But if I extend this to the case of a signal or spectrum that is
defined on any interval of length N, and where the i corresponding to
the 'zero' (time or frequency) is i0, then this interval is n = [ (1 -
i0 ) : (N - i0 ) ] and what is needed to translate this to matlab's
interval with n( 1 ) == 0 is a circular shift by n0.

Thus we need to do:

circshift( h, [ 0, ( 1 - i0 ) ] );

which is matlab-ese for a circular shift of h left by ( i0 - 1 ) - the
-1 because matlab starts array indices at 1..

Thus, by explicitly stating the matlab array index of the sample that
corresponds to 'zero' time or frequency, any interval of length N may
be translated correctly. And the matlab array index of 'zero' time
could alternatively be implicitly calculated from a vector of the
sample numbers (n) or from a time- or frequency- base.

This is one way I would expect to see this addressed in a slightly
object-oriented DSP architecture, where a 'signal' would also contain
metadata such as its 'zero' time or frequency position (offset if you
like).

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

On Jun 23, 11:30&#4294967295;pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 23, 5:04&#4294967295;am, Chris Bore <chris.b...@gmail.com> wrote:
>
> > I've seen discussion here of the correct use of MATLAB's fftshift()
> > function
>
> Please do not send separate copies of the same post to different
> newsgroups. Just post one copy but include all group names,
> separated by commas, on the newsgroup line.
>
> 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 &#4294967295;= ifftshift(xb)
> X &#4294967295;= fft(x)
> Xb = fftshift(X)
>
> and to return
>
> X &#4294967295;= ifftshift(Xb)
> x &#4294967295;= 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

Thanks to all contributors. And apologies for mis-posting (I have been
away from comp.dsp for a while and I seem to have clicked the wrong
'reply' several times).

Can I state the issue another way, and try to suggest a more generic
idiom?

The issue, I think, is that the fft() and ifft() in MATLAB are defined
on an interval n = [0 : N-1] where N is the number of samples, and n
is the 'sample number' that has n ==0 at the 'center' (time or
frequency). Whereas, I addressed the problem of a signal or spectrum
that is defined on the interval n = [ -N/2 : (N/2 - 1) ].

First, because the FFT assumes a periodic function, we can use a
circular shift to translate from one interval to another.
Second, the ifftshift() does such a circular shift for the interval I
specified (with appropriate caveats as to correct use)

Third, the issue is partly that it is easy to confuse the matrix index
(call this 'i' where in matlab i = [1 : N ] ) and the 'sample number'
n (which are not the same).

So in matlab's native usage, n( 1 ) == 0 whereas in my interval n( N/2
+ 1 ) == 0. That is, to re-state the obvious , the 'centers' are at
index i == 1 and i == (N/2+1) respectively.

But if I extend this to the case of a signal or spectrum that is
defined on any interval of length N, and where the i corresponding to
the 'zero' (time or frequency) is i0, then this interval is n = [ (1 -
i0 ) : (N - i0 ) ] and what is needed to translate this to matlab's
interval with n( 1 ) == 0 is a circular shift by n0.

Thus we need to do:

circshift( h, [ 0, ( 1 - i0 ) ] );

which is matlab-ese for a circular shift of h left by ( i0 - 1 ) - the
-1 because matlab starts array indices at 1..

Thus, by explicitly stating the matlab array index of the sample that
corresponds to 'zero' time or frequency, any interval of length N may
be translated correctly. And the matlab array index of 'zero' time
could alternatively be implicitly calculated from a vector of the
sample numbers (n) or from a time- or frequency- base.

This is one way I would expect to see this addressed in a slightly
object-oriented DSP architecture, where a 'signal' would also contain
metadata such as its 'zero' time or frequency position (offset if you
like). That would permit a 'generic' dft that would automatically
apply the correct shift, based on the offset, and so would save people
some potential difficulites.

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

On Jun 24, 12:56 am, Richard Owlett <rowl...@pcnetinc.com> wrote:
> Greg Heath wrote:
> > *SNIP*
>
> > Please do not send separate copies of the same post to different
> > newsgroups. Just post one copy but include all group names,
> > separated by commas, on the newsgroup line.
>
> Did you follow your own advice ?   *NOT* lol
> Was it possible? ??? ????? [I'm using SeaMonkey 2.0.5 as reader]
>
> May I suggest posting to all 'relevant' groups with "Reply To:"
> set to *ONLY* one group [that group noted in body]
>
> I now decease from "rant"


Whenever I am the OP I try to follow the comma-separated group list
procedure.

However, when I am just a replier, sometimes things get confused.
In particular, I have written one or more lengthy replies in one
group
only to find, later, that the OP has sent identical separate posts to
several other groups.

I post from Google Groups.

When I am the OP, there are 3 lines:

To Newsgroups:
CC:
Subject:

When I am a replier there is only one:

Newsgroups:

So, if I want previous replies to be seen by members of all
of the other groups AND, all subsequent replies to be sent
to all of the groups, I have to repost my original replies with
the appropriate comma-separated newsgroup list.

For the current case I just wanted the comp.dspers to see
my previous reply. I wasn't concerned about replies to my
reply.

>May I suggest posting to all 'relevant' groups with "Reply To:"
>set to *ONLY* one group [that group noted in body]

Unable to do it posting via Google Groups.

>I now decease from "rant"

No problem, a good purge every now and then is
recommended for mental stability.

Thanks for caring.

Greg

Richard wrote:
>I now decease from "rant"
>

Heart attack?

Jun 25, 2010 11:27:45 AM, chris.bore@gmail.com wrote:

On Jun 23, 11:30 pm, Greg Heath wrote:
> On Jun 23, 5:04 am, Chris Bore wrote:
>
> > I've seen discussion here of the correct use of MATLAB's fftshift()
> > function
>
> Please do not send separate copies of the same post to different
> newsgroups. Just post one copy but include all group names,
> separated by commas, on the newsgroup line.
>
> 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

>Thanks to all contributors.
>
>Can I state the issue another way, and try to suggest a more generic
>idiom?

>The issue, I think, is that the fft() and ifft() in MATLAB are defined
>on an interval n = [0 : N-1] where N is the number of samples, and n
>is the 'sample number' that has n ==0 at the 'center' (time or
>frequency).

No.

You have just defined n==0 to be at the beginning of the interval.

>Whereas, I addressed the problem of a signal or spectrum
>that is defined on the interval n = [ -N/2 : (N/2 - 1) ].

Your index n is only derived correctly for N even.

You are just further complicating the issue by introducing
the discrete index variable n. Now you have to get away from
the basic message by making a distinction between indexing based
on counting from 0 and indexing based on counting from 1.

If you reread my reply, I never explicitly mention the indexing
variable. The important point is, for N even or odd,

tb = dt*[-ceil(N-1)/2) : floor((N-1)/2) ]   % "b"ipolar

and

t  = dt*[ 0 : N-1 ]

regardless of how you index.

Similarly for fb and f with df = 1/(N*dt).


>First, because the FFT assumes a periodic function, we can use a
>circular shift to translate from one interval to another.
>Second, the ifftshift() does such a circular shift for the interval I
>specified (with appropriate caveats as to correct use)

Now you are trivializing a most important point:
The correct usage of fftshift and ifftshift is key.

>Second, the issue is partly that it is easy to confuse the matrix
>index (call this 'i' where in matlab i = [1 : N ] )

Arghhh!

Please not use i (or j) as an index. It will get confused with
sqrt(-1).

> and the 'sample number' n (which are not the same).

That confusion is exactly why I explained it without getting
into the distraction of the type of indexing variable that is
used.

>So in matlab's native usage, n( 1 ) == 0 whereas in my interval n( N/2
>+ 1 ) == 0. That is, to re-state the obvious , the 'centers' are at
>index i == 1 and i == (N/2+1) respectively.

...except when N is odd.

In general (without specifying the indexing convention),

tb(ceil(N+1)/2) = 0  % exanple: tb = dt*[ -2 -1 0 1 2 ]
t(1) = 0                  % example   t = dt*[ 0 1 2 -2 -1 ]

>But if I extend this to the case of a signal or spectrum that is
>defined on any interval of length N, and where the i corresponding to
>the 'zero' (time or frequency) is i0, then this interval is n = [ (1 -
>i0 ) : (N - i0 ) ] and what is needed to translate this to matlab's
>interval with n( 1 ) == 0 is a circular shift by n0.


>Thus we need to do:

>circshift( h, [ 0, ( 1 - i0 ) ] );

>which is matlab-ese for a circular shift of h left by ( i0 - 1 ) - the
>-1 because matlab starts array indices at 1..

>Thus, by explicitly stating the matlab array index of the sample that
>corresponds to 'zero' time or frequency, any interval of length N may
>be translated correctly. And the matlab array index of 'zero' time
>could alternatively be implicitly calculated from a vector of the
>sample numbers (n) or from a time- or frequency- base.

>This is one way I would expect to see this addressed in a slightly
>object-oriented DSP architecture, where a 'signal' would also contain
>metadata such as its 'zero' time or frequency position (offset if you
>like).

In general, consider x(t) defined over the arbitrary interval

t = t0 + dt*[0:(N-1)]

If t0 ~= 0 or t0 ~= -dt*ceil((N-1)/2)

just make a change of variable

s = (t-t0)

Then

X(f) = exp(i*2*pi*f*t0)*fft(x(s)).

This makes more sense than dealing with a more complicated
shift function.

Hope this helps.

Greg