Sine Lookup Table with Linear Interpolation

I've been studying an approach to implementing a lookup table (LUT) to 
implement a sine function.  The two msbs of the phase define the 
quadrant.  I have decided that an 8 bit address for a single quadrant is 
sufficient with an 18 bit output.  Another 11 bits of phase will give me 
sufficient resolution to interpolate the sin() to 18 bits.

If you assume a straight line between the two endpoints the midpoint of 
each interpolated segment will have an error of
((Sin(high)-sin(low))/2)-sin(mid)

Without considering rounding, this reaches a maximum at the last segment 
before the 90 degree point.  I calculate about 4-5 ppm which is about 
the same as the quantization of an 18 bit number.

There are two issues I have not come across regarding these LUTs.  One 
is adding a bias to the index before converting to the sin() value for 
the table.  Some would say the index 0 represents the phase 0 and the 
index 2^n represents 90 degrees.  But this is 2^n+1 points which makes a 
LUT inefficient, especially in hardware.  If a bias of half the lsb is 
added to the index before converting to a sin() value the value 0 to 
2^n-1 becomes symmetrical with 2^n to 2^(n+1)-1 fitting a binary sized 
table properly.  I assume this is commonly done and I just can't find a 
mention of it.

The other issue is how to calculate the values for the table to give the 
best advantage to the linear interpolation.  Rather than using the exact 
match to the end points stored in the table, an adjustment could be done 
to minimize the deviation over each interpolated segment.  Without this, 
the errors are always in the same direction.  With an adjustment the 
errors become bipolar and so will reduce the magnitude by half (approx). 
  Is this commonly done?  It will require a bit of computation to get 
these values, but even a rough approximation should improve the max 
error by a factor of two to around 2-3 ppm.

Now if I can squeeze another 16 dB of SINAD out of my CODEC to take 
advantage of this resolution!  lol

One thing I learned while doing this is that near 0 degrees the sin() 
function is linear (we all knew that, right?) but near 90 degrees, the 
sin() function is essentially quadratic.  Who would have thunk it?

-- 

Rick

this *should* be a relatively simple issue, but i am confused

On 3/16/13 6:30 PM, rickman wrote:
> I've been studying an approach to implementing a lookup table (LUT) to
> implement a sine function. The two msbs of the phase define the
> quadrant. I have decided that an 8 bit address for a single quadrant is
> sufficient with an 18 bit output.

10 bits or 1024 points.  since you're doing linear interpolation, add 
one more, copy the zeroth point x[0] to the last x[1024] so you don't 
have to do any modulo (by ANDing with (1023-1) on the address of the 
second point.  (probably not necessary for hardware implementation.)


     x[n]  = sin( (pi/512)*n )      for  0 <= n <= 1024

> Another 11 bits of phase will give me
> sufficient resolution to interpolate the sin() to 18 bits.

so a 21 bit total index.  your frequency resolution would be 2^(-21) in 
cycles per sampling period  or 2^(-21) * Fs.  those 2 million values 
would be the only frequencies you can meaningful

>
> If you assume a straight line between the two endpoints the midpoint of
> each interpolated segment will have an error of
> ((Sin(high)-sin(low))/2)-sin(mid)
>

do you mean "+" instead of the first "-"?  to be explicit:

  ( sin((pi/512)*(n+1)) + sin((pi/512)*n) )/2  -  sin((pi/512)*(n+0.5))

that's the error in the middle.  dunno if it's the max error, but it's 
might be.

> Without considering rounding, this reaches a maximum at the last segment
> before the 90 degree point.

at both the 90 and 270 degree points. (or just before and after those 
points.)

> I calculate about 4-5 ppm which is about the
> same as the quantization of an 18 bit number.
>
> There are two issues I have not come across regarding these LUTs. One is
> adding a bias to the index before converting to the sin() value for the
> table. Some would say the index 0 represents the phase 0 and the index
> 2^n represents 90 degrees. But this is 2^n+1 points which makes a LUT
> inefficient, especially in hardware. If a bias of half the lsb is added
> to the index before converting to a sin() value the value 0 to 2^n-1
> becomes symmetrical with 2^n to 2^(n+1)-1 fitting a binary sized table
> properly. I assume this is commonly done and I just can't find a mention
> of it.

do you mean biasing by 1/2 of a point?  then your max error will be *at* 
the 90 and 270 degree points and it will be slightly more than what you 
had before.

>
> The other issue is how to calculate the values for the table to give the
> best advantage to the linear interpolation. Rather than using the exact
> match to the end points stored in the table, an adjustment could be done
> to minimize the deviation over each interpolated segment. Without this,
> the errors are always in the same direction. With an adjustment the
> errors become bipolar and so will reduce the magnitude by half (approx).
> Is this commonly done? It will require a bit of computation to get these
> values, but even a rough approximation should improve the max error by a
> factor of two to around 2-3 ppm.

if you assume an approximate quadratic behavior over that short segment, 
you can compute the straight line where the error in the middle is equal 
in magnitude (and opposite in sign) to the error at the end points. 
that's a closed form solution, i think.

dunno if that is what you actually want for a sinusoidal waveform 
generator.  i might think you want to minimize the mean square error.

> Now if I can squeeze another 16 dB of SINAD out of my CODEC to take
> advantage of this resolution! lol
>
> One thing I learned while doing this is that near 0 degrees the sin()
> function is linear (we all knew that, right?) but near 90 degrees, the
> sin() function is essentially quadratic. Who would have thunk it?

Newton?  Leibnitz?  Gauss?

   sin(t + pi/2)  =  cos(t)


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

On Mar 16, 11:30&#4294967295;pm, rickman <gnu...@gmail.com> wrote:
> I've been studying an approach to implementing a lookup table (LUT) to
> implement a sine function. &#4294967295;The two msbs of the phase define the
> quadrant. &#4294967295;I have decided that an 8 bit address for a single quadrant is
> sufficient with an 18 bit output. &#4294967295;Another 11 bits of phase will give me
> sufficient resolution to interpolate the sin() to 18 bits.
>
> If you assume a straight line between the two endpoints the midpoint of
> each interpolated segment will have an error of
> ((Sin(high)-sin(low))/2)-sin(mid)
>
> Without considering rounding, this reaches a maximum at the last segment
> before the 90 degree point. &#4294967295;I calculate about 4-5 ppm which is about
> the same as the quantization of an 18 bit number.
>
> There are two issues I have not come across regarding these LUTs. &#4294967295;One
> is adding a bias to the index before converting to the sin() value for
> the table. &#4294967295;Some would say the index 0 represents the phase 0 and the
> index 2^n represents 90 degrees. &#4294967295;But this is 2^n+1 points which makes a
> LUT inefficient, especially in hardware. &#4294967295;If a bias of half the lsb is
> added to the index before converting to a sin() value the value 0 to
> 2^n-1 becomes symmetrical with 2^n to 2^(n+1)-1 fitting a binary sized
> table properly. &#4294967295;I assume this is commonly done and I just can't find a
> mention of it.
>
> The other issue is how to calculate the values for the table to give the
> best advantage to the linear interpolation. &#4294967295;Rather than using the exact
> match to the end points stored in the table, an adjustment could be done
> to minimize the deviation over each interpolated segment. &#4294967295;Without this,
> the errors are always in the same direction. &#4294967295;With an adjustment the
> errors become bipolar and so will reduce the magnitude by half (approx).
> &#4294967295; Is this commonly done? &#4294967295;It will require a bit of computation to get
> these values, but even a rough approximation should improve the max
> error by a factor of two to around 2-3 ppm.
>
> Now if I can squeeze another 16 dB of SINAD out of my CODEC to take
> advantage of this resolution! &#4294967295;lol
>
> One thing I learned while doing this is that near 0 degrees the sin()
> function is linear (we all knew that, right?) but near 90 degrees, the
> sin() function is essentially quadratic. &#4294967295;Who would have thunk it?
>
> --
>
> Rick

why  not skip the lut and just do a full sine approximation

something like this, page 53-54

http://www.emt.uni-linz.ac.at/education/Inhalte/pr_mderf/downloads/ADSP-2181/Using%20the%20ADSP-2100%20Family%20Volume%201/Using%20the%20ADSP-2100%20Family%20Volume%201.pdf


-Lasse

In comp.dsp rickman <gnuarm@gmail.com> wrote:

> I've been studying an approach to implementing a lookup table (LUT) to 
> implement a sine function.  The two msbs of the phase define the 
> quadrant.  I have decided that an 8 bit address for a single quadrant is 
> sufficient with an 18 bit output.  Another 11 bits of phase will give me 
> sufficient resolution to interpolate the sin() to 18 bits.

In the early days of MOS ROMs, there were commercial such ROMs,
I believe from National. (And when ROMs were much smaller than today.)

The data sheet has (I can't find it right now, but it is around
somewere) the combination of ROMs and TTL adders to do the
interpolation.
 
> If you assume a straight line between the two endpoints the midpoint of 
> each interpolated segment will have an error of
> ((Sin(high)-sin(low))/2)-sin(mid)
 
> Without considering rounding, this reaches a maximum at the last segment 
> before the 90 degree point.  I calculate about 4-5 ppm which is about 
> the same as the quantization of an 18 bit number.

I presume the ROM designers had all this figured out.
 
> There are two issues I have not come across regarding these LUTs.  One 
> is adding a bias to the index before converting to the sin() value for 
> the table.  Some would say the index 0 represents the phase 0 and the 
> index 2^n represents 90 degrees.  But this is 2^n+1 points which makes a 
> LUT inefficient, especially in hardware.  If a bias of half the lsb is 
> added to the index before converting to a sin() value the value 0 to 
> 2^n-1 becomes symmetrical with 2^n to 2^(n+1)-1 fitting a binary sized 
> table properly.  I assume this is commonly done and I just can't find a 
> mention of it.

I don't remember anymore. But since you have the additional bit to
do the interpolation, it should be easy. 
 
> The other issue is how to calculate the values for the table to give the 
> best advantage to the linear interpolation.  Rather than using the exact 
> match to the end points stored in the table, an adjustment could be done 
> to minimize the deviation over each interpolated segment.  Without this, 
> the errors are always in the same direction.  With an adjustment the 
> errors become bipolar and so will reduce the magnitude by half (approx). 
>  Is this commonly done?  It will require a bit of computation to get 
> these values, but even a rough approximation should improve the max 
> error by a factor of two to around 2-3 ppm.

I am not sure how you are thinking about doing it. I believe that
some of the bits that go into the MSB ROM also go into the lower ROM
to select the interpolation slope, and then additional bits to select
the actual value. Say you have a 1024x10 ROM for the first one,
then want to interpolate that. How many bits of linear interpolation
can be done? (That is, before linear isn't close enough any more.)
Then the appropriate number of low bits and high bits, but not the
in between bits, go into the interpolation ROM, which is then added
to the other ROMs output.

> Now if I can squeeze another 16 dB of SINAD out of my CODEC to take 
> advantage of this resolution!  lol
 
> One thing I learned while doing this is that near 0 degrees the sin() 
> function is linear (we all knew that, right?) but near 90 degrees, the 
> sin() function is essentially quadratic.  Who would have thunk it?

Well, cos() is known to be quadratic near zero.

I once knew someone with a homework problem something like:

Take a string all the way around the earth at the equator, and add
(if I remember right) 3m. (Assuming the earth is a perfect sphere.)
If the string is at a uniform height around the earth, how far from
the surface is it?

Now, pull it up at one point. How far is that point above the surface?
More specifically, as I originally heard it, is it higher than the
height of a specific nine story library?

As I remember it, the usual small angle approximations to trig.
functions aren't enough to do this. The next term is needed.
The 3m added might not be right, as the answer is close enough 
to the certain building height to need the additional term.

-- glen

rickman wrote:
> I've been studying an approach to implementing a lookup table (LUT) to
> implement a sine function. The two msbs of the phase define the
> quadrant. I have decided that an 8 bit address for a single quadrant is
> sufficient with an 18 bit output. Another 11 bits of phase will give me
> sufficient resolution to interpolate the sin() to 18 bits.
>
> If you assume a straight line between the two endpoints the midpoint of
> each interpolated segment will have an error of
> ((Sin(high)-sin(low))/2)-sin(mid)
>
> Without considering rounding, this reaches a maximum at the last segment
> before the 90 degree point. I calculate about 4-5 ppm which is about the
> same as the quantization of an 18 bit number.
>
> There are two issues I have not come across regarding these LUTs. One is
> adding a bias to the index before converting to the sin() value for the
> table. Some would say the index 0 represents the phase 0 and the index
> 2^n represents 90 degrees. But this is 2^n+1 points which makes a LUT
> inefficient, especially in hardware. If a bias of half the lsb is added
> to the index before converting to a sin() value the value 0 to 2^n-1
> becomes symmetrical with 2^n to 2^(n+1)-1 fitting a binary sized table
> properly. I assume this is commonly done and I just can't find a mention
> of it.
>
> The other issue is how to calculate the values for the table to give the
> best advantage to the linear interpolation. Rather than using the exact
> match to the end points stored in the table, an adjustment could be done
> to minimize the deviation over each interpolated segment. Without this,
> the errors are always in the same direction. With an adjustment the
> errors become bipolar and so will reduce the magnitude by half (approx).
> Is this commonly done? It will require a bit of computation to get these
> values, but even a rough approximation should improve the max error by a
> factor of two to around 2-3 ppm.
>
> Now if I can squeeze another 16 dB of SINAD out of my CODEC to take
> advantage of this resolution! lol
>
> One thing I learned while doing this is that near 0 degrees the sin()
> function is linear (we all knew that, right?) but near 90 degrees, the
> sin() function is essentially quadratic. Who would have thunk it?
>
   Sounds like you are making excellent improvements on the standard 
ho-hum algorithms; the net result will be superior to anything done out 
there (commercially).
   With the proper offsets, one needs only 22.5 degrees of lookup 
("bounce" off each multiple of 45 degrees).

On 3/16/2013 5:30 PM, rickman wrote:
> I've been studying an approach to implementing a lookup table (LUT) to
> implement a sine function.   The two msbs of the phase define the
> quadrant.  I have decided that an 8 bit address for a single quadrant is
> sufficient with an 18 bit output.  Another 11 bits of phase will give me
> sufficient resolution to interpolate the sin() to 18 bits.

An optimal one quadrant LUT with 256 entries and linear interpolation 
makes for sine approximation with accuracy about 16 bits.
If there is a need for better precision, consider different approaches.

Vladimir Vassilevsky
DSP and Mixed Signal Designs
www.abvolt.com

In comp.dsp rickman <gnuarm@gmail.com> wrote:
> I've been studying an approach to implementing a lookup table (LUT) to 
> implement a sine function.  The two msbs of the phase define the 
> quadrant.  I have decided that an 8 bit address for a single quadrant is 
> sufficient with an 18 bit output.  Another 11 bits of phase will give me 
> sufficient resolution to interpolate the sin() to 18 bits.
 
> If you assume a straight line between the two endpoints the midpoint of 
> each interpolated segment will have an error of
> ((Sin(high)-sin(low))/2)-sin(mid)

See:

http://ia601506.us.archive.org/8/items/bitsavers_nationaldaMOSIntegratedCircuits_20716690/1972_National_MOS_Integrated_Circuits.pdf

The description starts on page 273.

-- glen

On Sat, 16 Mar 2013 18:30:51 -0400, rickman wrote:

> There are two issues I have not come across regarding these LUTs.  One 
> is adding a bias to the index before converting to the sin() value for 
> the table.  Some would say the index 0 represents the phase 0 and the 
> index 2^n represents 90 degrees.  But this is 2^n+1 points which makes a 
> LUT inefficient, especially in hardware.

An n-bit table has 2^n+1 entries for 2^n ranges. Range i has endpoints of
table[i] and table[i+1]. The final range has i=(1<<n)-1, so the last
entry in the table is table[1<<n], not table[(1<<n)-1].

> One thing I learned while doing this is that near 0 degrees the sin() 
> function is linear (we all knew that, right?) but near 90 degrees, the 
> sin() function is essentially quadratic.  Who would have thunk it?

sin((pi/2)+x) = sin((pi/2)-x) = cos(x), and the Maclaurin series for
cos(x) is:

	cos(x) = 1 - (x^2)/2! + (x^4)/4! - ...

On 3/16/2013 7:13 PM, robert bristow-johnson wrote:
>
> this *should* be a relatively simple issue, but i am confused
>
> On 3/16/13 6:30 PM, rickman wrote:
>> I've been studying an approach to implementing a lookup table (LUT) to
>> implement a sine function. The two msbs of the phase define the
>> quadrant. I have decided that an 8 bit address for a single quadrant is
>> sufficient with an 18 bit output.
>
> 10 bits or 1024 points. since you're doing linear interpolation, add one
> more, copy the zeroth point x[0] to the last x[1024] so you don't have
> to do any modulo (by ANDing with (1023-1) on the address of the second
> point. (probably not necessary for hardware implementation.)
>
>
> x[n] = sin( (pi/512)*n ) for 0 <= n <= 1024

So you are suggesting a table with 2^n+1 entries?  Not such a great idea 
in some apps, like hardware.  What is the advantage?  Also, why 10 bit 
address for a 1024 element table?  My calculations indicate a linear 
interpolation can be done with 4 ppm accuracy with a 256 element LUT. 
I'm not completely finished my simulation, but I'm pretty confident this 
much is corrrect.


>> Another 11 bits of phase will give me
>> sufficient resolution to interpolate the sin() to 18 bits.
>
> so a 21 bit total index. your frequency resolution would be 2^(-21) in
> cycles per sampling period or 2^(-21) * Fs. those 2 million values would
> be the only frequencies you can meaningful

No, that is the phase sent to the LUT.  The total phase accumulator can 
be larger as the need requires.


>> If you assume a straight line between the two endpoints the midpoint of
>> each interpolated segment will have an error of
>> ((Sin(high)-sin(low))/2)-sin(mid)
>>
>
> do you mean "+" instead of the first "-"? to be explicit:
>
> ( sin((pi/512)*(n+1)) + sin((pi/512)*n) )/2 - sin((pi/512)*(n+0.5))
>
> that's the error in the middle. dunno if it's the max error, but it's
> might be.

Yes, thanks for the correction.  The max error?  I'm not so worried 
about that exactly.  The error is a curve with the max magnitude near 
the middle if nothing further is done to minimize it.


>> Without considering rounding, this reaches a maximum at the last segment
>> before the 90 degree point.
>
> at both the 90 and 270 degree points. (or just before and after those
> points.)

I'm talking about the LUT.  The LUT only considers the first quadrant.


>> I calculate about 4-5 ppm which is about the
>> same as the quantization of an 18 bit number.
>>
>> There are two issues I have not come across regarding these LUTs. One is
>> adding a bias to the index before converting to the sin() value for the
>> table. Some would say the index 0 represents the phase 0 and the index
>> 2^n represents 90 degrees. But this is 2^n+1 points which makes a LUT
>> inefficient, especially in hardware. If a bias of half the lsb is added
>> to the index before converting to a sin() value the value 0 to 2^n-1
>> becomes symmetrical with 2^n to 2^(n+1)-1 fitting a binary sized table
>> properly. I assume this is commonly done and I just can't find a mention
>> of it.
>
> do you mean biasing by 1/2 of a point? then your max error will be *at*
> the 90 and 270 degree points and it will be slightly more than what you
> had before.

No, not quite right.  There is a LUT with points spaced at 90/255 
degrees apart starting at just above 0 degrees.  The values between 
points in the table are interpolated with a maximum deviation near the 
center of the interpolation.  Next to 90 degrees the interpolation is 
using the maximum interpolation factor which will result in a value as 
close as you can get to the correct value if the end points are used to 
construct the interpolation line.  90 degrees itself won't actually be 
represented, but rather points on either side, 90&#4294967295;delta where delta is 
360&#4294967295; / 2^(n+1) with n being the number of bits in the input to the sin 
function.


>> The other issue is how to calculate the values for the table to give the
>> best advantage to the linear interpolation. Rather than using the exact
>> match to the end points stored in the table, an adjustment could be done
>> to minimize the deviation over each interpolated segment. Without this,
>> the errors are always in the same direction. With an adjustment the
>> errors become bipolar and so will reduce the magnitude by half (approx).
>> Is this commonly done? It will require a bit of computation to get these
>> values, but even a rough approximation should improve the max error by a
>> factor of two to around 2-3 ppm.
>
> if you assume an approximate quadratic behavior over that short segment,
> you can compute the straight line where the error in the middle is equal
> in magnitude (and opposite in sign) to the error at the end points.
> that's a closed form solution, i think.

Yes, it is a little tricky because at this point we are working with 
integer math (or technically fixed point I suppose).  Rounding errors is 
what this is all about.  I've done some spreadsheet simulations and I 
have some pretty good results.  I updated it a bit to generalize it to 
the LUT size and I keep getting the same max error counts (adjusted to 
work with integers rather than fractions) &#4294967295;3 no matter what the size of 
the interpolation factor.  I don't expect this and I think I have 
something wrong in the calculations.  I'll need to resolve this.


> dunno if that is what you actually want for a sinusoidal waveform
> generator. i might think you want to minimize the mean square error.

We are talking about the lsbs of a 20+ bit word.  Do you think there 
will be much of a difference in result?  I need to actually be able to 
do the calculations and get this done rather than continue to work on 
the process.  Also, each end point affects two lines, so there are 
tradeoffs, make one better and the other worse?  It seems to get 
complicated very quickly.


>> Now if I can squeeze another 16 dB of SINAD out of my CODEC to take
>> advantage of this resolution! lol
>>
>> One thing I learned while doing this is that near 0 degrees the sin()
>> function is linear (we all knew that, right?) but near 90 degrees, the
>> sin() function is essentially quadratic. Who would have thunk it?
>
> Newton? Leibnitz? Gauss?
>
> sin(t + pi/2) = cos(t)

How does that imply a quadratic curve at 90 degrees?  At least I think 
like the greats!

-- 

Rick

On 3/16/2013 8:02 PM, langwadt@fonz.dk wrote:
> On Mar 16, 11:30 pm, rickman<gnu...@gmail.com>  wrote:
>> I've been studying an approach to implementing a lookup table (LUT) to
>> implement a sine function.  The two msbs of the phase define the
>> quadrant.  I have decided that an 8 bit address for a single quadrant is
>> sufficient with an 18 bit output.  Another 11 bits of phase will give me
>> sufficient resolution to interpolate the sin() to 18 bits.
>>
>> If you assume a straight line between the two endpoints the midpoint of
>> each interpolated segment will have an error of
>> ((Sin(high)-sin(low))/2)-sin(mid)
>>
>> Without considering rounding, this reaches a maximum at the last segment
>> before the 90 degree point.  I calculate about 4-5 ppm which is about
>> the same as the quantization of an 18 bit number.
>>
>> There are two issues I have not come across regarding these LUTs.  One
>> is adding a bias to the index before converting to the sin() value for
>> the table.  Some would say the index 0 represents the phase 0 and the
>> index 2^n represents 90 degrees.  But this is 2^n+1 points which makes a
>> LUT inefficient, especially in hardware.  If a bias of half the lsb is
>> added to the index before converting to a sin() value the value 0 to
>> 2^n-1 becomes symmetrical with 2^n to 2^(n+1)-1 fitting a binary sized
>> table properly.  I assume this is commonly done and I just can't find a
>> mention of it.
>>
>> The other issue is how to calculate the values for the table to give the
>> best advantage to the linear interpolation.  Rather than using the exact
>> match to the end points stored in the table, an adjustment could be done
>> to minimize the deviation over each interpolated segment.  Without this,
>> the errors are always in the same direction.  With an adjustment the
>> errors become bipolar and so will reduce the magnitude by half (approx).
>>    Is this commonly done?  It will require a bit of computation to get
>> these values, but even a rough approximation should improve the max
>> error by a factor of two to around 2-3 ppm.
>>
>> Now if I can squeeze another 16 dB of SINAD out of my CODEC to take
>> advantage of this resolution!  lol
>>
>> One thing I learned while doing this is that near 0 degrees the sin()
>> function is linear (we all knew that, right?) but near 90 degrees, the
>> sin() function is essentially quadratic.  Who would have thunk it?
>>
>> --
>>
>> Rick
>
> why  not skip the lut and just do a full sine approximation
>
> something like this, page 53-54
>
> http://www.emt.uni-linz.ac.at/education/Inhalte/pr_mderf/downloads/ADSP-2181/Using%20the%20ADSP-2100%20Family%20Volume%201/Using%20the%20ADSP-2100%20Family%20Volume%201.pdf

I'm not sure this would be easier.  The LUT an interpolation require a 
reasonable amount of logic.  This requires raising X to the 5th power 
and five constant multiplies.  My FPGA doesn't have multipliers and this 
may be too much for the poor little chip.  I suppose I could do some of 
this with a LUT and linear interpolation... lol

-- 

Rick