Prev: Dmitry Terez frequency estimation method.
Next: OMAP DSP starting vector
From: Chris Bore on 23 Jun 2010 05:04
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

From: Rune Allnor on 23 Jun 2010 07:38
On 23 Jun, 11:04, Chris Bore <chris.b...(a)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
From: robert bristow-johnson on 23 Jun 2010 15:58
On Jun 23, 5:04 am, Chris Bore <chris.b...(a)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 –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’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.
>
....
> 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():
>
....
>
>   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
From: Greg Heath on 23 Jun 2010 18:30
On Jun 23, 5:04 am, Chris Bore <chris.b...(a)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


From: Chris Bore on 25 Jun 2010 11:23
On Jun 23, 11:30 pm, Greg Heath <he...(a)alumni.brown.edu> wrote:
> On Jun 23, 5:04 am, Chris Bore <chris.b...(a)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

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




 |  Next  |  Last
Pages: 1 2
Prev: Dmitry Terez frequency estimation method.
Next: OMAP DSP starting vector