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 =96N/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=92s complement style) from =96N/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 =3D fftshift( fft ( h ) );

Similarly, the inverse fft also requires this fftshift:

  h =3D 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 =3D=3D 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 =3D 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 =3D 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 =3D=3D 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 =3D 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 =3D fftshift( fft( h ) );                  % most common usage,
wrong result
   H =3D 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 =3D 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
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
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=A0am, 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 =96N/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=92s complement style) from =96N/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():
>
....
>
> =A0 H =3D 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 =3D some_big_index_into_array_x;

   X =3D 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 =3D 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=A0am, 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 =3D dt*(0:N-1) and f =3D 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 =3D dt * [-ceil((N-1)/2) : floor((N-1)/2)]
fb =3D 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  =3D ifftshift(xb)
X  =3D fft(x)
Xb =3D fftshift(X)

and to return

X  =3D ifftshift(Xb)
x  =3D ifft(X)
xb =3D 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=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 23, 5:04=A0am, 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 =3D dt*(0:N-1) and f =3D 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 =3D dt * [-ceil((N-1)/2) : floor((N-1)/2)]
> fb =3D 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 =A0=3D ifftshift(xb)
> X =A0=3D fft(x)
> Xb =3D fftshift(X)
>
> and to return
>
> X =A0=3D ifftshift(Xb)
> x =A0=3D ifft(X)
> xb =3D 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 =3D [0 : N-1] where N is the number of samples, and n
is the 'sample number' that has n =3D=3D0 at the 'center' (time or
frequency). Whereas, I addressed the problem of a signal or spectrum
that is defined on the interval n =3D [ -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 =3D [1 : N ] ) and the 'sample
number' n (which are not the same).

So in matlab's native usage, n( 1 ) =3D=3D 0 whereas in my interval n( N/2
+ 1 ) =3D=3D 0. That is, to re-state the obvious , the 'centers' are at
index i =3D=3D 1 and i =3D=3D (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 =3D [ (1 -
i0 ) : (N - i0 ) ] and what is needed to translate this to matlab's
interval with n( 1 ) =3D=3D 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
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
Chris Bore
BORES Signal Processing
www.bores.com

On Jun 23, 11:30=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 23, 5:04=A0am, 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 =3D dt*(0:N-1) and f =3D 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 =3D dt * [-ceil((N-1)/2) : floor((N-1)/2)]
> fb =3D 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 =A0=3D ifftshift(xb)
> X =A0=3D fft(x)
> Xb =3D fftshift(X)
>
> and to return
>
> X =A0=3D ifftshift(Xb)
> x =A0=3D ifft(X)
> xb =3D 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 =3D [0 : N-1] where N is the number of samples, and n
is the 'sample number' that has n =3D=3D0 at the 'center' (time or
frequency). Whereas, I addressed the problem of a signal or spectrum
that is defined on the interval n =3D [ -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 =3D [1 : N ] ) and the 'sample number'
n (which are not the same).

So in matlab's native usage, n( 1 ) =3D=3D 0 whereas in my interval n( N/2
+ 1 ) =3D=3D 0. That is, to re-state the obvious , the 'centers' are at
index i =3D=3D 1 and i =3D=3D (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 =3D [ (1 -
i0 ) : (N - i0 ) ] and what is needed to translate this to matlab's
interval with n( 1 ) =3D=3D 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
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
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

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

On Jun 26, 1:11=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> Jun 25, 2010 11:27:45 AM, chris.b...@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 =3D dt*(0:N-1) and f =3D 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 =3D dt * [-ceil((N-1)/2) : floor((N-1)/2)]
> > fb =3D 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 =A0=3D ifftshift(xb)
> > X =A0=3D fft(x)
> > Xb =3D fftshift(X)
>
> > and to return
>
> > X =A0=3D ifftshift(Xb)
> > x =A0=3D ifft(X)
> > xb =3D 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 =3D [0 : N-1] where N is the number of samples, and n
> >is the 'sample number' that has n =3D=3D0 at the 'center' (time or
> >frequency).
>
> No.
>
> You have just defined n=3D=3D0 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 =3D [ -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 =3D dt*[-ceil(N-1)/2) : floor((N-1)/2) ] =A0 % "b"ipolar
>
> and
>
> t =A0=3D dt*[ 0 : N-1 ]
>
> regardless of how you index.
>
> Similarly for fb and f with df =3D 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 =3D [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 ) =3D=3D 0 whereas in my interval n( N=
/2
> >+ 1 ) =3D=3D 0. That is, to re-state the obvious , the 'centers' are at
> >index i =3D=3D 1 and i =3D=3D (N/2+1) respectively.
>
> ...except when N is odd.
>
> In general (without specifying the indexing convention),
>
> tb(ceil(N+1)/2) =3D 0 =A0% exanple: tb =3D dt*[ -2 -1 0 1 2 ]
> t(1) =3D 0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% example =A0 t =3D 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 =3D [ (1 -
> >i0 ) : (N - i0 ) ] and what is needed to translate this to matlab's
> >interval with n( 1 ) =3D=3D 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 =3D t0 + dt*[0:(N-1)]
>
> If t0 ~=3D 0 or t0 ~=3D -dt*ceil((N-1)/2)
>
> just make a change of variable
>
> s =3D (t-t0)
>
> Then
>
> X(f) =3D 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

On Jun 25, 8:43=A0pm, "Nasser M. Abbasi" <n...@12000.org> wrote:
> 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.

no shit.

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

my bottom lip is nearly bleeding, i am biting it so hard.

> I wish Matlab would allow one to change the indexing.

now it *is* bleeding.

> In Fortran for
> example, one can do that,

in the *new* Fortran.  because they allowed it to evolve.

Cleve (Moler), did you know that languages can evolve?  even MATLAB?

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

i think so.  but it's still amazing that the overlooked this problem
then and continue to do so now.  it was more excusable then.

> 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=3D0 to n=3DN-1, and the index n is al=
so
> 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.

NO IT ISN'T!!  and it wasn't 10 years ago, even though they tried to
imply to us that it was.  they can *extend* the language by
incorporating into every MATLAB variable a vector which has length
equal to the number of dimensions that's just like the vector they
already have that defines the numbers of rows and columns and whatever
they call it for the 3rd or 4th dimension.  except this vector would
define the index number of the "first" (and i don't mean the "one_th")
element in each dimension.  if, when a MATLAB variable is first
created, if these origin values are assigned the default value of "1",
it would be perfectly backward compatible and no legacy code would
break.  then there would be functions similar to size() and reshape()
that would instead be called base() and rebase() and would allow the
user to examine and change (respectively) these origin values.  THAT
WOULD FIX THIS STUPID PROBLEM THAT CLEVE DENIES (or, at least
steadfastly denied it a decade ago) IS A PROBLEM.

r b-j

On Sat, 26 Jun 2010 18:40:49 -0700 (PDT), robert bristow-johnson
<rbj@audioimagination.com> wrote:

>On Jun 25, 8:43=A0pm, "Nasser M. Abbasi" <n...@12000.org> wrote:
>> 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.
>
>no shit.
>
>> 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.
>
>my bottom lip is nearly bleeding, i am biting it so hard.
>
>> I wish Matlab would allow one to change the indexing.
>
>now it *is* bleeding.
>
>> In Fortran for
>> example, one can do that,
>
>in the *new* Fortran.  because they allowed it to evolve.
>
>Cleve (Moler), did you know that languages can evolve?  even MATLAB?
>
>> 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?)
>
>i think so.  but it's still amazing that the overlooked this problem
>then and continue to do so now.  it was more excusable then.
>
>> 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=3D0 to n=3DN-1, and the index n is al=
>so
>> 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.
>
>NO IT ISN'T!!  and it wasn't 10 years ago, even though they tried to
>imply to us that it was.  they can *extend* the language by
>incorporating into every MATLAB variable a vector which has length
>equal to the number of dimensions that's just like the vector they
>already have that defines the numbers of rows and columns and whatever
>they call it for the 3rd or 4th dimension.  except this vector would
>define the index number of the "first" (and i don't mean the "one_th")
>element in each dimension.  if, when a MATLAB variable is first
>created, if these origin values are assigned the default value of "1",
>it would be perfectly backward compatible and no legacy code would
>break.  then there would be functions similar to size() and reshape()
>that would instead be called base() and rebase() and would allow the
>user to examine and change (respectively) these origin values.  THAT
>WOULD FIX THIS STUPID PROBLEM THAT CLEVE DENIES (or, at least
>steadfastly denied it a decade ago) IS A PROBLEM.
>
>r b-j

IMHO the most plausible explanation for why this has never been
addressed is that a conscious decision has been made within the
MathWorks that Matlab should not have flexible or zero-based indexing
capability.

Matlab has evolved a fair amount, as they've added structured features
and a few other things that probably appeal to a few people, so it's
not like they're just averse to changing anything.


Eric Jacobsen
Minister of Algorithms
Abineau Communications
http://www.abineau.com

On 27 Jun, 07:56, eric.jacob...@ieee.org (Eric Jacobsen) wrote:

> IMHO the most plausible explanation for why this has never been
> addressed is that a conscious decision has been made within the
> MathWorks that Matlab should not have flexible or zero-based indexing
> capability.

The reason why it has never been addressed is that MATLAB is
an acronym for MATrix LABoratory. In linear algebra the numerical
arrays are - like it or not - indexed base 1.

Changing this would undermine nearly 30 years worth of code base.
So the main decision to be made at TMW is whether they want to
sell a marginally predictable product on one side, or piss off
RBJ on the other.

The choise might be easier than I suspect RBJ appreciates.

Rune

On Jun 27, 8:09=A0am, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 27 Jun, 07:56, eric.jacob...@ieee.org (Eric Jacobsen) wrote:
>
> > IMHO the most plausible explanation for why this has never been
> > addressed is that a conscious decision has been made within the
> > MathWorks that Matlab should not have flexible or zero-based indexing
> > capability.
>
> The reason why it has never been addressed is that MATLAB is
> an acronym for MATrix LABoratory. In linear algebra the numerical
> arrays are - like it or not - indexed base 1.

whether it's 0 or 1 is not the issue.  whether the user can define it
for a particular array *is* the issue.

> Changing this would undermine nearly 30 years worth of code base.

no it wouldn't.  it could be perfectly backward compatible, because
newly-created arrays would have default origin of 1 for every
dimension of the array.  the user would have to call a not-yet-
existing function to change the origin.

> So the main decision to be made at TMW is whether they want to
> sell a marginally predictable product on one side,

it's still just as deterministic as before.  no one would be
randomizing the origins of arrays.

> or piss off RBJ on the other.

it's not just me, Rune.

r b-j

On 27 Jun, 20:06, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Jun 27, 8:09=A0am, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > On 27 Jun, 07:56, eric.jacob...@ieee.org (Eric Jacobsen) wrote:
>
> > > IMHO the most plausible explanation for why this has never been
> > > addressed is that a conscious decision has been made within the
> > > MathWorks that Matlab should not have flexible or zero-based indexing
> > > capability.
>
> > The reason why it has never been addressed is that MATLAB is
> > an acronym for MATrix LABoratory. In linear algebra the numerical
> > arrays are - like it or not - indexed base 1.
>
> whether it's 0 or 1 is not the issue. =A0whether the user can define it
> for a particular array *is* the issue.

In that case, why can't you just use the OO capabilities of
matlab and declare your own array class, say, RBJarray, and
overload the () operator?

Ought to be a piece of cake, provided matlab's OO is half decent.

Rune

Apologies to those who might think I'm stealing their thread.

*CAVEAT LECTOR* Though I've *NEVER* supported myself by any trade 
normally associated with DSP, my avocations regularly deal with 
topics covered by this group.

I've dealt with problems best expressed with a 0 base index.
I've dealt with problems best expressed with a 1 base index.

For problems involving a FFT (or inverse) I use Scilab 4.x
Choice NOT based on indexing - couldn't afford MATLAB and at the 
time Octave had installation problems on *my* Windows setup.

I think it is foolish for for any purported "GENERAL PURPOSE" 
package to not address the index base issue.

A before I get really sarcastic comments - If I'm not old enough 
to have changed their diapers/wipes, I'm more than old enough to 
have been their "sitter" ;}

On Jun 27, 2:16=A0pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 27 Jun, 20:06, robert bristow-johnson <r...@audioimagination.com>
> wrote:
>
> > On Jun 27, 8:09=A0am, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > > On 27 Jun, 07:56, eric.jacob...@ieee.org (Eric Jacobsen) wrote:
>
> > > > IMHO the most plausible explanation for why this has never been
> > > > addressed is that a conscious decision has been made within the
> > > > MathWorks that Matlab should not have flexible or zero-based indexi=
ng
> > > > capability.
>
> > > The reason why it has never been addressed is that MATLAB is
> > > an acronym for MATrix LABoratory. In linear algebra the numerical
> > > arrays are - like it or not - indexed base 1.
>
> > whether it's 0 or 1 is not the issue. =A0whether the user can define it
> > for a particular array *is* the issue.
>
> In that case, why can't you just use the OO capabilities of
> matlab and declare your own array class, say, RBJarray, and
> overload the () operator?

i think this is described in subsref() and subsasgn() and in this page
 http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_oop/br09eqz.h=
tml

i had thought about it before.  i think i would call the class
"rebase" or "reorigin", because a call to something that looks like a
function called that would be what creates an object in that class
(with a regular MATLAB array as an argument as well as the new origin
values).  this psuedo-function would look like a counterpart to
"reshape".  then there would be another basic function called "base"
or "origin" that would return a vector with the origins (just like
"size" returns the values set with "reshape").

the issue is "just".  never marketed myself as an OOPs person.  i have
asked people at comp.soft-sys.matlab about this, oh, maybe 5 or 8
years ago.  can't remember what came of it.

speed is an issue.  and cleanliness (much cleaner done at the root
level since MATLAB has to subtract 1 from indices (at least for arrays
having more than one dimension) for individual  to assemble a net
pointer into the array.

> Ought to be a piece of cake, provided matlab's OO is half decent.

if you see it as a piece of cake, i would invite you to collaborate
since i am shit-for-brains regarding the nasty details of setting this
up.

a decade ago, i posted to comp.soft-sys.matlab a description on how
array arithmetic would also have to be modified (by overloading "+",
"-", "*", ".*", "\", "/", "./", "^", ".^" and all of the simple math
functions that work on a straight-forward element-by-element basis).
the array created would have the origins defined in the correct manner
(not always at 1).  Google Groups doesn't work as good as it used to
but i could try to dig up those posts.

that would be a beginning.  a different fft (i would call them "dft"
and "idft") that would look at the origins and apply the proper phase
modification (or rearrange what goes in, similarly to fftshift).
different functions for ployval() (and related functions like
polyfit()) that get the order of coefficients wrong (should be the 0th-
order or constant coef first, at the 0th index).

it would be a task, but doable to start.  still should instead be done
at the fundamental level, just like "complex" should be a basic type
in C.  it might be a piece of cake for a seasoned MATLAB programmer
who's used their OO (and maybe subsref and subsasgn) before, but i
just want to be able to do DSP and translate the common equations
directly from the lit to MATLAB without having to fiddle with the
addition or subtraction of 1 for the fft or some other constant if i
want to model a non-causal impulse response.  i should be able to just
do that directly.

if this is a piece of cake, i'd be interested in collaboration.  (but
it would still be better if it were done at the fundamental level,
without having to create a class.)

r b-j

On 27 Jun, 22:03, robert bristow-johnson <r...@audioimagination.com>
wrote:
>
> > Ought to be a piece of cake, provided matlab's OO is half decent.
>
> if you see it as a piece of cake, i would invite you to collaborate
> since i am shit-for-brains regarding the nasty details of setting this
> up.

Sorry, can't help with that. I've never really progressed
from how I started using matlab some 20 years ago, as an
algorithm testbed and general prototyping tool. Once things
turn 'real' I tend to switch to C++ rather quickly. But in
C++ what you want is trivial. For instance, switching from
a base 0 vector indexing to base 1 would be something like
(allowing for the element dereference operator being [] rather
than () in C++)

template<class T>
class base1vector : public std::vector<T>{
public:
   T& operator[](size_t i){return this->operator[](i-1);};
};

You would probably need a few variations over this theme,
and some constructors, but the gist is a simple as this.

Rune

On Jun 27, 3:22=A0pm, Richard Owlett <rowl...@pcnetinc.com> wrote:
> Apologies to those who might think I'm stealing their thread.
>
> *CAVEAT LECTOR* Though I've *NEVER* supported myself by any trade
> normally associated with DSP, my avocations regularly deal with
> topics covered by this group.
>
> I've dealt with problems best expressed with a 0 base index.
> I've dealt with problems best expressed with a 1 base index.
>
> For problems involving a FFT (or inverse) I use Scilab 4.x
> Choice NOT based on indexing - couldn't afford MATLAB and at the
> time Octave had installation problems on *my* Windows setup.
>
> I think it is foolish for for any purported "GENERAL PURPOSE"
> package to not address the index base issue.
>
> A before I get really sarcastic comments - If I'm not old enough
> to have changed their diapers/wipes, I'm more than old enough to
> have been their "sitter" ;}

with me, Richard, you're preaching to the choir.  i find this deficit
in MATLAB to be astonishing.  particularly so long after multiple
people have pointed it out (including an early developer for the Sig
Proc toolbox).

r b-j

robert bristow-johnson wrote:
> On Jun 27, 8:09 am, Rune Allnor <all...@tele.ntnu.no> wrote:

>> Changing this would undermine nearly 30 years worth of code base.

> no it wouldn't.  it could be perfectly backward compatible, because
> newly-created arrays would have default origin of 1 for every
> dimension of the array.  the user would have to call a not-yet-
> existing function to change the origin.

Robert, your enthusiasm for this idea has led you to neglect thinking 
the problems through.

Suppose I have a mex file written to the current API. One or more of the 
arguments to the function are array indices. Now introduce your proposed 
new indexing scheme and have the user create such an array and pass it 
in to the existing mex file.

Your claim that the proposed change "could be perfectly backwards 
compatible" implies that the indices the user passes in must not rely on 
the new indexing scheme, because "perfectly backwards compatible" means 
that old code must work UNCHANGED. Now *without changing the API for 
existing routines*, how is Matlab going to know if a (say) 5 being 
passed in as a numeric value is already adjusted to be 1-relative or 
needs to be silently re-biased by Matlab to the appropriate basis in 
order to preserve backwards compatibility? What if the dimension to be 
indexed is itself is a parameter so that the unbiasing that needs to 
take place is not constant? What if the dimension number is not an 
_obvious_  parameter, such as if you had some encoding such as the mesh 
encoding that mixes control values and data values in the same array? 
What if the indices that have to be rebiased have been packed, such as 
two 8-bit indices numerically jammed into a 16 bit number -- how is 
Matlab going to know the jamming algorithm to know how to rebias and 
construct the appropriate new index?

As long as indices are computable data then in order to support 
different index biases you MUST break backwards compatibility, in that 
the existing code would have to be enhanced to know about and take into 
account the new indexing scheme for any parameter that is not provably 
an old-style matrix.

On Jul 4, 12:42=A0pm, Walter Roberson <rober...@hushmail.com> wrote:
> robert bristow-johnson wrote:
> > On Jun 27, 8:09 am, Rune Allnor <all...@tele.ntnu.no> wrote:
> >> Changing this would undermine nearly 30 years worth of code base.
> >
> > no it wouldn't. =A0it could be perfectly backward compatible, because
> > newly-created arrays would have default origin of 1 for every
> > dimension of the array. =A0the user would have to call a not-yet-
> > existing function to change the origin.
>
> Robert, your enthusiasm for this idea has led you to neglect thinking
> the problems through.

that might be a premature judgment.  you don't know how much i have
thought this through for more than a decade.  if you can get Google
Groups to search adequately, you might find the places where the
objections (similar to yours) were brought up and i swatted them
down.  even Cleve eventually admitted that it was, from the strict
definition of the term, backward compatible.

> Suppose I have a mex file written to the current API.

as long as no one applies your MEX function (written under the old
assumptions) to any array with any origin not 1, there would be no
problem.  still backward compatible.  no one's code breaks.

if this extension or enhancement to MATLAB were adopted and you wanted
your .mex function to work with arrays of origin different than 1, you
would have to modify the .mex function to look for the origins (using
the new API).  otherwise your function would work as if all of the
origins were 1 when they may not be.

it's no different than with any other extension or enhancement to a
language.  the issue you brought up is not a violation of backward
compatibility.  old code (with old .mex files) would still work the
same way they did before.

r b-j

On Jun 27, 9:53=A0pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 27 Jun, 22:03, robert bristow-johnson <r...@audioimagination.com>
> wrote:
>
>
>
> > > Ought to be a piece of cake, provided matlab's OO is half decent.
>
> > if you see it as a piece of cake, i would invite you to collaborate
> > since i am shit-for-brains regarding the nasty details of setting this
> > up.
>
> Sorry, can't help with that. I've never really progressed
> from how I started using matlab some 20 years ago, as an
> algorithm testbed and general prototyping tool. Once things
> turn 'real' I tend to switch to C++ rather quickly. But in
> C++ what you want is trivial. For instance, switching from
> a base 0 vector indexing to base 1 would be something like
> (allowing for the element dereference operator being [] rather
> than () in C++)
>
> template<class T>
> class base1vector : public std::vector<T>{
> public:
> =A0 =A0T& operator[](size_t i){return this->operator[](i-1);};
>
> };
>
> You would probably need a few variations over this theme,
> and some constructors, but the gist is a simple as this.
>
> Rune

Yes. Things have turned real for few months now for me. I even tried
to write a "number crunching manager" that would regularily write down
data to disk to free up matlab workspace but didn't get me too far
with my 22hours processing time. I highly recommend the "templated" C+
+ approach. You can use generics to even move down from floating-point
to fixed-point using SystemC headers. What more could  you possibly
ask for?

-Momo

"robert bristow-johnson" <rbj@audioimagination.com> wrote in message 
news:2f984955-7c2c-4414-8b5b-5843e314a8d7@b35g2000yqi.googlegroups.com...
> On Jul 4, 12:42 pm, Walter Roberson <rober...@hushmail.com> wrote:
> > robert bristow-johnson wrote:
> > > On Jun 27, 8:09 am, Rune Allnor <all...@tele.ntnu.no> wrote:
> > >> Changing this would undermine nearly 30 years worth of code base.
> > >
> > > no it wouldn't. it could be perfectly backward compatible, because
> > > newly-created arrays would have default origin of 1 for every
> > > dimension of the array. the user would have to call a not-yet-
> > > existing function to change the origin.
> >
> > Robert, your enthusiasm for this idea has led you to neglect thinking
> > the problems through.
>
> that might be a premature judgment.  you don't know how much i have
> thought this through for more than a decade.  if you can get Google
> Groups to search adequately, you might find the places where the
> objections (similar to yours) were brought up and i swatted them
> down.  even Cleve eventually admitted that it was, from the strict
> definition of the term, backward compatible.
>
> > Suppose I have a mex file written to the current API.
>
> as long as no one applies your MEX function (written under the old
> assumptions) to any array with any origin not 1, there would be no
> problem.  still backward compatible.  no one's code breaks.

Except for the person who DOES apply an older MEX function to a non-1 based 
array.  People often need/want to reuse their old code, and I've found 
people tend to get a wee bit upset if you tell them "Sorry, don't do that." 
unless there's a major benefit they can immediately see (and sometimes even 
then.)  While some (including yourself) may see 0-based indexing as such a 
major benefit, there are also many who would not.

> if this extension or enhancement to MATLAB were adopted and you wanted
> your .mex function to work with arrays of origin different than 1, you
> would have to modify the .mex function to look for the origins (using
> the new API).  otherwise your function would work as if all of the
> origins were 1 when they may not be.

And that is an incompatibility.  Let's take a look at a few others.

1)
Let's say that a user had a function that needs to loop over the columns of 
their matrix.  This is not a contrived example; I've seen plenty of code 
that does something similar.

for whichcolumn = 1:size(A, 2)
    % process column A(:, whichcolumn)
end

If A is a 1-based array, then everything works as it did before.  If A was a 
0-based array, then this would error with an error message like "Index 
exceeds matrix dimension" or something similar.  Code that used to work no 
longer works.  That's an incompatibility.

Even worse, let's say the user's code, for whatever reason, only needed to 
process all but the final column of A.

for whichcolumn = 1:size(A, 2)-1
    % process column A(:, whichcolumn)
end

Now this code can SILENTLY GIVE THE WRONG ANSWER if A is a 0-based array, as 
it will process all but the _first_ column of A.

2)
What happens if I SAVE a 0-based array in a MAT-file and LOAD it in a 
version prior to the introduction of this change?  Would you expect that to 
work?

3)
Suppose A is a 0-based array and B is a 1-based array.  What happens if I 
ask for A+B?  A.*B?  [Assume all the sizes match; just consider the 
based-ness for this scenario.]

I'm guessing you're going to say that A+B and A.*B both error; operating 
under that assumption, let's take a look at one more example.  What should 
the command pi+A do?  Does the pi built-in function return a 0-based array 
or a 1-based array?  Seems to me that users would expect pi+A to add pi to 
each element of A -- but for backwards compatibility, and to agree with your 
own proposal from earlier in this thread, PI would _have_ to return a 
1-based array and so pi+A would ERROR.  If you say that PI should be "smart 
enough" to know with what it's going to be combined and return the 
appropriate-indexed scalar:

x = pi;
y = A+x;
z = B+x;

_One_ of the latter two operations must error unless scalars are treated as 
a special case.  If they are treated specially, x = 1:10; will encounter the 
same problem assuming A and B are both 10-element vectors, since x must have 
a specified base when it is created.

> it's no different than with any other extension or enhancement to a
> language.  the issue you brought up is not a violation of backward
> compatibility.  old code (with old .mex files) would still work the
> same way they did before.

Robert, when you've brought this type of system up in the past I've 
seriously thought about how it could work, but based on the potential 
problems I called out above, speaking personally I do not think it would be 
a good idea to modify MATLAB indexing to use your proposed system.

I think your best option would be to create your own 0-based (or 
variable-based) object and overload those functions with which you want to 
work as well as subscripted indexing and assignment -- that way you can 
control the behavior of indexing for your object.  You can even make use of 
the built-in functions (rather than having to reimplement them) by using the 
BUILTIN function to call them on the 1-based data that's stored inside the 
array as a private data member, and adjust the indices afterward (if 
necessary.)  Indeed, Nabeel posted an object to do just this to the File 
Exchange several years ago:

http://www.mathworks.com/matlabcentral/fileexchange/1168-varbase

If you wanted to write a classdef-based version of this object then the 
Object-Oriented Programming section of the documentation will contain the 
information you'll need to become familiar with this style of object.  In 
particular, since indexing will be a major component of this object, the 
following page will be of special interest:

http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_oop/br09eqz.html

-- 
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on 
http://www.mathworks.com 

On Jul 6, 3:38=A0pm, "Steven Lord" <sl...@mathworks.com> wrote:
> "robert bristow-johnson" <r...@audioimagination.com> wrote in message
> news:2f984955-7c2c-4414-8b5b-5843e314a8d7@b35g2000yqi.googlegroups.com...
>
>
>
> > On Jul 4, 12:42 pm, Walter Roberson <rober...@hushmail.com> wrote:
> > > robert bristow-johnson wrote:
> > > > On Jun 27, 8:09 am, Rune Allnor <all...@tele.ntnu.no> wrote:
> > > >> Changing this would undermine nearly 30 years worth of code base.
>
> > > > no it wouldn't. it could be perfectly backward compatible, because
> > > > newly-created arrays would have default origin of 1 for every
> > > > dimension of the array. the user would have to call a not-yet-
> > > > existing function to change the origin.
>
> > > Robert, your enthusiasm for this idea has led you to neglect thinking
> > > the problems through.
>
> > that might be a premature judgment. =A0you don't know how much i have
> > thought this through for more than a decade. =A0if you can get Google
> > Groups to search adequately, you might find the places where the
> > objections (similar to yours) were brought up and i swatted them
> > down. =A0even Cleve eventually admitted that it was, from the strict
> > definition of the term, backward compatible.
>
> > > Suppose I have a mex file written to the current API.
>
> > as long as no one applies your MEX function (written under the old
> > assumptions) to any array with any origin not 1, there would be no
> > problem. =A0still backward compatible. =A0no one's code breaks.
>
> Except for the person who DOES apply an older MEX function to a non-1 bas=
ed
> array.

the definition of "backward compatible" that i am using is (from
Wikipedia): "a product or a technology is said to be backward
compatible when it is able to fully take the place of an older
product... Backward compatibility is a relationship between two
components, rather than being an attribute of just one of them. More
generally, a new component is said to be backward compatible if it
provides all of the functionality of the old component."

what i am and had been proposing for a decade is backward compatible
in that meaning.  if people try to use the "new feature" (non-1 based
arrays) and misuse them and get error messages, that does not mean
that it fails backward compatibility.

if a MEX function was unaware of this (and old ones *would* be
unaware), it would treat any array argument as if it had origin 1 for
every dimension.  it would be wrong, but that would be a misuse (to
use a function never written to deal with other origins on a non-1
based array).  and nothing would blow up.

> Robert, when you've brought this type of system up in the past I've
> seriously thought about how it could work, but based on the potential
> problems I called out above, speaking personally I do not think it would =
be
> a good idea to modify MATLAB indexing to use your proposed system.
>
> I think your best option would be to create your own 0-based (or
> variable-based) object and overload those functions with which you want t=
o
> work as well as subscripted indexing and assignment -- that way you can
> control the behavior of indexing for your object.

here's the deal:  i am a DSP person, not an OOP person.  MATLAB has
marketed itself (falsely) making claims (in the v4 and v5 user
manuals) as:

"MATLAB integrates numerical analysis, matrix computation, signal
processing, and graphics in an easy-to-use environment where problems
and
solutions are expressed just as they are written mathematically - ...
"

"MATLAB is a high-performance language for technical computing.  It
integrates computation, visualization, and programming in an easy-to-
use
environment where problems and solutions are expressed in familiar
mathematical notation."

note the phrases in claims: "familiar mathematical notation" and "just
as
they are written mathematically".  i submit that those claims are
false in
a sense that is salient particularly for those of use who use MATLAB
for
(digital) signal processing.  and i suspect the claims are false for
some
other users that also deal with data that are naturally sequenced with
subscripts or indices that are non-positive.


now, what i want is for MATLAB to live up to that.  MATLAB is not C++
where, if i complain that C doesn't give me a complex variable type,
the first response from someone is that i should use C++ and a library
with a complex class.  MATLAB *does* have an OOP capability and with
those two subscript calls (i think they're called subsref() and
subsasgn()) to make it possible to insert a shim that intercepts the
indices.  so, what you say Steven is true, a class can be done, but
what i need to use MATLAB (or Octave) for is to do DSP.  i should not
have to become an OOPs programmer just to do basic DSP with equations
that are recognizable in the DSP lit.

i remember that Thomas Krauss once told me (a decade ago) that he
brought this up early before the Sig Proc Toolbox (the very first
MATLAB toolbox) was completed and released.  this is really when you
guys should have fixed the problem.  you have forced people to adopt
non-standard non-natural indexing for the problems they work on.
besides the FFT (where it should be obvious) i think TMW blew it with
the definitions of polyval() and related functions like polyfit().
they should have used 0-based arrays of coefficients with a(n) being
the coefficient for x^n.  similarly TMW also screwed up the
definitions of conv() and related functions.  when polynomials are
multiplied to each other, we know it's like their coefficient
sequences are convolved and the FFT can be used to convolve large
coefficient sequences.  but the order of coefficients (and the index
values associated) are just plain wrong in MATLAB.  besides missing
the proper 0-based counting, TMW put the order wrong.

> =A0You can even make use of
> the built-in functions (rather than having to reimplement them) by using =
the
> BUILTIN function to call them on the 1-based data that's stored inside th=
e
> array as a private data member, and adjust the indices afterward (if
> necessary.)

if we were to do it that way, what the object converter (i would call
it "rebase" or "reorigin") should do, is check to see if the origins
for all dimensions in the output array are all 1, then it should
return a regular MATLAB matrix variable and not one of this "rebase"
type) so that all of the existing MATLAB functions can work on this 1-
based array.

> =A0Indeed, Nabeel posted an object to do just this to the File
> Exchange several years ago:

i remember.  and i remember i was unable to use it (because much more
needed to be done).

> http://www.mathworks.com/matlabcentral/fileexchange/1168-varbase
>
> If you wanted to write a classdef-based version of this object then the
> Object-Oriented Programming section of the documentation will contain the
> information you'll need to become familiar with this style of object. =A0=
In
> particular, since indexing will be a major component of this object, the
> following page will be of special interest:
>
> http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_oop/br09...

again, you're requiring me to be an OOPs programmer when what i need
MATLAB for is to do signal processing.  this is something *you*guys*
should fix.  i would be happy to help in specification (and i had and
i'll dig up the old text) but this is a deficit in MATLAB and it's
sorta non-responsible for TMW to require users to fix such deficits in
their product.

that varbase class should be renamed "reorigin" and, to begin with,
the following operators should be overloaded according to the
following spec, and if a resulting "reorigin" object happens to have
only 1 for the origins of every dimension, then a regular 1-based
MATLAB variable should be returned.  i know it's a little bit to read,
but i would ask that you read it.  to make that class work, this is
the minimum that is needed to make it work usefully and correctly (and
backward compatible).  i am, of course, guessing at what the C-like
date structure is for a MATLAB variable, but whatever it is, the
pointer to the origin[] vector should be appended to the end or use an
existing *unused* field.  then it would be backward compatible with
MEX files.

r b-j

________________________________________________


enum MATLAB_class {text, real, complex};
// I don't wanna cloud the issue considering other classes.

typedef struct
   {
   void* data;     // pointer to actual array data
   char* name;     // pointer to the variable's name
   enum MATLAB_class type; // class of MATLAB variable (real,
complex,...)
   int num_dimensions;     // number of array dimensions >=3D 2
   long* size;    // points to a vector with the number of rows,
columns, etc.
   } MATLAB_variable;

char name[32];       // suppose MATLAB names are unique to 31 chars
long size[num_dimensions];

if (type =3D=3D text)
   {
   char data[size[0]*size[1]*...*size[num_dimensions-1]];
   }
  else if (type =3D=3D real)
   {
   double data[size[0]*size[1]*...*size[num_dimensions-1];
   }
  else if (type =3D=3D complex)
   {
   double data[size[0]*size[1]*...*size[num_dimensions-1][2];
   }

The above is sorta C-like pseudocode.  I'm writing it as if the
declarations
can allocated like malloc() does.


Currently, when an element, A(n,k), of a 2 dimensional MATLAB array A
is
accessed, first n and k are confirmed to be integer value (not a
problem
in C), then confirmed to be at least 1 and less than or equal to
size[0]
and size[1], respectively.  It those constraints are satisfied, the
value
of that element is accessed as:

data[(n-1)*size[0] + (k-1)];

For a 3 dimensional array, A(m,n,k), it would be the same but now:

data[((m-1)*size[1] + (n-1))*size[0] + (k-1)];

What is proposed is to first add a new member to the MATLAB variable
structure called "origin" which is a vector of the very same length
(num_dimensions) as the "size" vector.  The default value for all
elements
of the origin[] vector would be 1 with only the exceptions outlined
below.
This is what makes this backwards compatible, in the strictest sense
of the
term.

typedef struct
 {
 void* data;
 char* name;
 enum MATLAB_class type;
 int num_dimensions;
 long* size;
 long* origin;  // points to a vector with index origin for each
dimension
 } MATLAB_variable;

char name[32];
long size[num_dimensions];
long origin[num_dimensions];

Now before each index is used, it is checked against the bounds for
that
dimension ( origin[dim] <=3D index < size[dim]+origin[dim]    where
0 <=3D dim < num_dimensions),  Since the default for origin[dim] is 1,
this
will have no effect, save for the teeny amount of processing time need
to
look up the origin, on existing MATLAB legacy code.

So to access a single element of an array A, this C array index would
look like:

data[(n-origin[1])*size[0] + (k-origin[0])];

For a 3 dimensional array, A(m,n,k), it would look like:

data[((m-origin[2])*size[1] + (n-origin[1])*size[0] + (k-origin[0])];

Okay, how someone like myself would use this to do something different
is
that there would be at least two new MATLAB facilities similar to
size() and
reshape() that I might call "origin()" and "reorigin()",
respectively.  Just
like MATLAB size() function returns the contents of the size[] vector,
origin() would return, in MATLAB format, the contents of the origin[]
vector.
And just like reshape() changes (under proper conditions) the contents
of
the size[] vector, reorigin() would change the contents of the
origin[] vector.
Since reorigin() does not exist in legacy MATLAB code (oh, I suppose
someone
could have created a function named that, but that's a naming problem
that
need not be considered here), then there is no way for existing MATLAB
programs to change the origins from their default values of 1 making
this
fix perfectly backward compatible.

Now, just as there are dimension compatibility rules that exist now
for
MATLAB operations, there would be a few natural rules that would be
added so
that "reorigined" MATLAB arrays could have operations applied to them
in a
sensible way.




ARRAY ADDITION ("+") and SUBTRACTION ("-") and element-by-element
ARRAY
MULTIPLICATION (".*"), DIVISION ("./"), POWER (".^"), and ELEMENTARY
FUNCTIONS:

Currently MATLAB insists that the number of dimensions are equal and
the
size of each dimension are equal (that is the same "shape") before
adding or
subtracting matrices or arrays.  The one exception to that is adding a
scaler to an array, in which a hypothetical array of equal size and
shape
with all elements equal to the scaler, is added to the array.  The
resulting
array has the same size and shape as the argument arrays.

The proposed system would, of course, continue this constraint and add
a new
constraint in that index origins for each dimension (the origin[]
vector)
would have to be equal for two arrays to be added or subtracted or
element-
by-element multiplied, etc.  The resulting array would have the same
shape
and origin[] vector as the input arrays.




"MATRIX" MULTIPLICATION ("*"):

A =3D B*C;

Currently MATLAB appropriately insists that the number of columns of B
are
equal to the number of rows of C (we shall call that number K).  The
resulting array has the number of rows of B and the number of columns
of C.
The value of a particular element of A would be:

          K
A(m,n) =3D SUM{ B(m,k) * C(k,n) }
         k=3D1

The proposed system would, of course, continue this constraint and add
a new
constraint in that index origins must be equal for each dimension
where the
lengths must be equal.  That is the number of columns of B are equal
to the
number of rows of C and the origin index of the columns of B are equal
to the
origin index of the rows of C.  The resulting array has the number of
rows of
B and the number of columns of C and  the origin index of the rows of
B and
the origin index of the columns of C.  The value of a particular
element of A
would be:

        org-1+K
A(m,n) =3D SUM{ B(m,k) * C(k,n) }
        k=3Dorg

where org =3D origin[0] for the B array and origin[1] for the C array
which must be the same number.  In the same manner as the present
MATLAB
K would be size[0] for the B array and size[1] for the C array and
would
have to be the same number, otherwize a diagnostic error would be
returned.

Both of these definitions are degenerations of the more general case
where:

         +inf
A(m,n) =3D SUM{ B(m,k) * C(k,n) }
        k=3D-inf

where here you consider B and C to be zero-extended to infinity in all
four
directions (up, down, left, and right).  It's just that the zero
element
pairs do not have to be multiplied and summed.

Matrix powers and exponentials (on square matrices) can be defined
to be consistent with this extension of the matrix multiply.


MATRIX DIVISION ("/" and "\"):

Like the current matrix division, it would invert the operation of

    A =3D B*C

that is

    C =3D B\A

and

    B =3D A/C

The same size (or shape) requirements and index origin requirements of
the "*" operator would apply to "\" and "/".   Given A and B, whatever
shape and origin requirement for B and C and whatever shape and origin
returned in A in the statement:

    A =3D B*C;

would be the same requirements for B and A and define the resulting
shape and origin for C in the statement:

    C =3D B\A



CONCATINATION:

This would also be a simple and straight-forward extension of how
MATLAB
presently concatinates arrays.  When we say:

A =3D [B C];

The number of rows of B and C must be equal, but the number of columns
of B
and C can be anything.  The first columns of A are identical with the
columns of B and then also must the indices of those columns.  And
independent of what the column indices of C are, they just pick up
where the
column index of B left off.  This rule extension defaults to what
MATLAB
presently does if B and C are both arrays with origin 1.  A similar
rule
extension can be made for A =3D [B ; C];  In all cases the upper left
corner
of A is identical to the upper left corner of B, both in value but
also in
subscripts (so A(1,1) becomes B(1,1) just like it does now in MATLAB).




FUNCTIONS THAT RETURN INDICES (min(), max(), find(), sort(),
ind2sub(), and
any others that I don't know about):

It must be internally consistent (and certainly can be made to be).
The
indices returned would be exactly like the 1-based indices returned
presently in MATLAB except that the origin for the corresponding
dimension
(that defaults to 1) would be added to each C-like index.  That is,
just
like now in MATLAB:

[max_value, max_index] =3D max(A);

This must mean that A(max_index) is equal to max_value.

I think that this is easy enough to define.  The only hard part is to
identify all MATLAB functions that search through an array and modify
them
to start and end at indices that might be different than 1 and
size[dim] as
are the search bounds today.  It would instead search from origin[dim]
to
size[dim]+origin[dim]-1 which would default to the current operation
if the
origin equals 1.




FOR ALL OTHER MATLAB OPERATIONS, until a reasonable extended
definition for
arrays with origins not 1 is thought up, MATLAB could either bomb out
with
an illegal operation error if the base is not 1 or could, perhaps,
ignore
the origin.  Either way, it's still backwards compatible.

because i copied some of this stuff from a decade old post, and
evidently, Google Groups has changed the number of characters per line
and we cannot control word wrapping in our posts in Google Groups, i
have edited the line endings (and we'll let the lines wrap
naturally).   this should be much more readable than my previous post.

Steven Lord and others at TMW, can you *please* read this carefully,
with enough of an open mind, and *then* respond with questions, legit
criticisms, suggestions and fixes.  and please be straight with us
about the definition of "backward compatible".  if *old* code works
with this, (without making use of any calls to "reorigin()" or
"rebase()" or whatever it would be named), it *is* "backward
compatible", no?

________________________________________________________________

(here are the quotes from old MATLAB manuals that *should* be the
uncompromized vision of the use of the product.  just because TMW no
longer prints these exact statements in their current manuals does not
remove that vision for what MATLAB should be.)

"MATLAB integrates numerical analysis, matrix computation, signal
processing, and graphics in an easy-to-use environment where problems
and solutions are expressed just as they are written
mathematically..."

"MATLAB is a high-performance language for technical computing.  It
integrates computation, visualization, and programming in an easy-to-
use environment where problems and solutions are expressed in familiar
mathematical notation."

note the phrases in claims: "familiar mathematical notation" and "just
as they are written mathematically".  i submit that those claims are
false in a sense that is salient particularly for those of use who use
MATLAB for (digital) signal processing.  and i suspect the claims are
false for some other users that also deal with data that are naturally
sequenced with subscripts or indices that are non-positive.


________________________________________________________________

(this is the foundation of the solution to this indexing problem that
i am asking TMW people to read carefully and not reject due to red
herrings.)




enum MATLAB_class {text, real, complex};
// I don't wanna cloud the issue considering other classes.

typedef struct
   {
   void* data;     // pointer to actual array data
   char* name;     // pointer to the variable's name
   enum MATLAB_class type; // class of MATLAB variable (real, etc.
   int num_dimensions;  // number of array dimensions >= 2
   long* size;    // vector with the number of rows, columns, etc.
   } MATLAB_variable;

char name[32];    // suppose MATLAB names are unique to 31 chars
long size[num_dimensions];

if (type == text)
   {
   char data[size[0]*size[1]*...*size[num_dimensions-1]];
   }
  else if (type == real)
   {
   double data[size[0]*size[1]*...*size[num_dimensions-1];
   }
  else if (type == complex)
   {
   double data[size[0]*size[1]*...*size[num_dimensions-1][2];
   }

The above is sorta C-like pseudocode.  I'm writing it as if the
declarations can allocated like malloc() does.

Currently, when an element, A(n,k), of a 2-dimensional MATLAB array A
is accessed, first n and k are confirmed to be integer value (not a
problem in C), then confirmed to be at least 1 and less than or equal
to size[0] and size[1], respectively.  If those constraints are
satisfied, the value of that element is accessed as:

  data[(n-1)*size[0] + (k-1)];

For a 3 dimensional array, A(m,n,k), it would be the same but now:

  data[((m-1)*size[1] + (n-1))*size[0] + (k-1)];



What is proposed is to first add a new member to the MATLAB variable
structure called "origin" which is a vector of the very same length
(num_dimensions) as the "size" vector.  The default value for all
elements of the origin[] vector would be 1 with only the exceptions
outlined below.  This is what makes this backwards compatible, in the
strictest sense of the term.

typedef struct
   {
   void* data;
   char* name;
   enum MATLAB_class type;
   int num_dimensions;
   long* size;
   long* origin;  // points to a vector with origin for each dimension
   } MATLAB_variable;

char name[32];
long size[num_dimensions];
long origin[num_dimensions];

Now before each index is used, it is checked against the bounds for
that dimension ( origin[dim] <= index < size[dim]+origin[dim] where 0
<= dim < num_dimensions).  Since the default for origin[dim] is 1,
this will have no effect, save for the teeny amount of processing time
need to look up the origin, on existing MATLAB legacy code.

So to access a single element of an array, A(n,k), this C array index
would look like:

  data[(n-origin[1])*size[0] + (k-origin[0])];

For a 3 dimensional array, A(m,n,k), it would look like:

  data[((m-origin[2])*size[1] + (n-origin[1])*size[0] + (k-
origin[0])];



Okay, how someone like myself would use this to do something different
is that there would be at least two new MATLAB facilities similar to
size() and reshape() that I might call "origin()" and "reorigin()",
respectively.  Just like the MATLAB size() function returns the
contents of the size[] vector, origin() would return, in MATLAB
format, the contents of the origin[] vector.

And just like reshape() changes (under proper conditions) the contents
of the size[] vector, reorigin() would change the contents of the
origin[] vector.  Since reorigin() does not exist in legacy MATLAB
code (oh, I suppose someone could have created a function named that,
but that's a naming problem that need not be considered here), then
there is no way for existing MATLAB programs to change the origins
from their default values of 1 making this fix perfectly backward
compatible.



Now, just as there are dimension compatibility rules that exist now
for MATLAB operations, there would be a few natural rules that would
be added so that "reorigined" MATLAB arrays could have operations
applied to them in a sensible way.




ARRAY ADDITION ("+") and SUBTRACTION ("-") and element-by-element
ARRAY MULTIPLICATION (".*"), DIVISION ("./"), POWER (".^"), and
ELEMENTARY FUNCTIONS:

Currently MATLAB insists that the number of dimensions are equal and
the size of each dimension are equal (that is the same "shape") before
adding or subtracting matrices or arrays.  The one exception to that
is adding a scaler to an array, in which a hypothetical array of equal
size and shape with all elements equal to the scaler, is added to the
array.  The resulting array has the same size and shape as the
argument arrays.

The proposed system would, of course, continue this constraint and add
a new constraint in that index origins for each dimension (the
origin[] vector) would have to be equal for two arrays to be added or
subtracted or element-by-element multiplied, etc.  The resulting array
would have the same size[] vector and origin[] vector as the input
arrays.




"MATRIX" MULTIPLICATION ("*"):

    A = B*C;

Currently MATLAB appropriately insists that the number of columns of B
are equal to the number of rows of C (we shall call that number K).
The resulting array has the number of rows of B and the number of
columns of C.  The value of a particular element of A would be:

              K
    A(m,n) = SUM{ B(m,k) * C(k,n) }
             k=1

The proposed system would, of course, continue this constraint and add
a new constraint in that index origins must be equal for each
dimension where the lengths must be equal.  That is the number of
columns of B are equal to the number of rows of C and the origin index
of the columns of B is equal to the origin index of the rows of C.
The resulting array, A, has the number of rows of B and the number of
columns of C and the origin index of the rows of B and the origin
index of the columns of C would also be the same as for A.  The value
of a particular element of A would be:

             org-1+K
     A(m,n) = SUM{ B(m,k) * C(k,n) }
            k=org

where org = origin[0] for the B array and origin[1] for the C array
which must be the same number.  In the same manner as the present
MATLAB K would be size[0] for the B array and size[1] for the C array
and would have to be the same number, otherwise a diagnostic error
would be returned.

Both of these definitions are degenerations of the more general case
where:

             +inf
    A(m,n) = SUM{ B(m,k) * C(k,n) }
            k=-inf

where here you consider B and C to be zero-extended to infinity in all
four directions (up, down, left, and right).  It's just that the zero
element pairs do not have to be multiplied and summed.

Matrix powers and exponentials (on square matrices) can be defined to
be consistent with this extension of the matrix multiply.




MATRIX DIVISION ("/" and "\"):

Like the current matrix division, it would invert the operation of

    A = B*C

that is

    C = B\A

and

    B = A/C

The same size (or shape) requirements and index origin requirements of
the "*" operator would apply to "\" and "/".   Given A and B, whatever
shape and origin requirement for B and C and whatever shape and origin
returned in A in the statement:

    A = B*C;

would be the same requirements for B and A and define the resulting
shape and origin for C in the statement:

    C = B\A




CONCATINATION:

This would also be a simple and straight-forward extension of how
MATLAB presently concatinates arrays.  When we say:

    A = [B C];

The number of rows of B and C must be equal, but the number of columns
of B and C can be anything.  The first columns of A are identical with
the columns of B and then also must the indices of those columns.  And
independent of what the column indices of C are, they just pick up
where the column index of B left off.  This rule extension defaults to
what MATLAB presently does if B and C are both arrays with origin 1.
A similar rule extension can be made for A = [B ; C];  In all cases
the upper left corner of A is identical to the upper left corner of B,
both in value but also in subscripts (so A(1,1) becomes B(1,1) just
like it does now in MATLAB).



FUNCTIONS THAT RETURN INDICES (min(), max(), find(), sort(),
ind2sub(), and any others that I don't know about):

It must be internally consistent (and certainly can be made to be).
The indices returned would be exactly like the 1-based indices
returned presently in MATLAB except that the origin for the
corresponding dimension (that defaults to 1) would be added to each C-
like (0-based) index.  That is, just like now in MATLAB:

   [max_value, max_index] = max(A);

This must mean that A(max_index) is equal to max_value.

I think that this is easy enough to define.  The only hard part is to
identify all MATLAB functions that search through an array and modify
them to start and end at indices that might be different than 1 and
size[dim] as are the search bounds today.  It would instead search
from origin[dim] to size[dim]+origin[dim]-1 which would default to the
current operation if the origin equals 1.



FOR ALL OTHER MATLAB OPERATIONS, until a reasonable extended
definition for arrays with origins not 1 is thought up, MATLAB could
either bomb out with an illegal operation error if the base is not 1
or could, perhaps, ignore the origin.  Either way, it's still
backwards compatible.

robert bristow-johnson <rbj@audioimagination.com> wrote in message <f3e6e4ef-8e06-4d8b-a4e6-392ed00594a4@d37g2000yqm.googlegroups.com>...
> 
    (snip)
> 
> Steven Lord and others at TMW, can you *please* read this carefully ...
>
    (snip)

To TMW:  If you guys ever *do* take the time to read through this and seriously consider doing something like this, I hope you put your ideas out ahead of time for some type of pier review to ferret out potential usage problems before implementing. e.g., Fortran had to do something similar when they went from strictly 1-based indexing to any-based indexing.

James Tursa