Prev: Dmitry Terez frequency estimation method.
Next: OMAP DSP starting vector
From: Chris Bore on 25 Jun 2010 11:36
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. 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
From: Greg Heath on 25 Jun 2010 20:11

Jun 25, 2010 11:27:45 AM, chris.bore(a)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
From: Nasser M. Abbasi on 25 Jun 2010 20:43
On 6/25/2010 5:11 PM, Greg Heath wrote:

>
> 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).
>

On a side point, and I do not mean to interrupt the main discussion
about fftshift, but wanted to mention the indexing issue.

I think the fact that in Matlab one can't start indexing at 0 can cause
one to easily make a programming error.

When one tries to code a mathematical equation from the textbook, which
uses 0 as an index, into matlab, and have to adjust things, one can make
the famous 1-off error.

I wish Matlab would allow one to change the indexing. In Fortran for
example, one can do that, and it can make implementing DSP algorithm
easier I think. (Strange that even though Matlab was originally
implemented in Fortran, if I remember correctly, that this was not put
into Matlab arrays? May because at the time Fortran did not allow 0
based arrays?)

As an exercise the other day, I wrote few lines of code to implement DFT
in Fortran 90 and Ada. Both of these allowed one to define an array with
an index that starts from 0, and I think the code was easier to
implement, since I did not have to worry about the 1-off problem as I
would have in Matlab or in any other language that does not have this
flexibility.

In DFT, the summation starts from n=0 to n=N-1, and the index n is also
used, inside the loop, to index into the array of samples. So, one have
to remember to add 1 to n before referencing the array. But in the
textbook equation, there was no such addition of 1.

Hence the Fortran and the Ada code matched the text book exactly, but
the Matlab code would not. This is just an example.

http://12000.org/my_notes/mma_matlab_control/KERNEL/node94.htm

I know that this point have been talked about may times before, and it
is too late now to change Matlab.

--Nasser
From: Chris Bore on 26 Jun 2010 13:07
On Jun 26, 1:11 am, Greg Heath <he...(a)alumni.brown.edu> wrote:
> Jun 25, 2010 11:27:45 AM, chris.b...(a)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

So why is this, which I agree makes complete sense, not done in matlab
normally but instead the fftshift() and ifftshift() functions are
used? That is the thing that initially puzzled me.

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