Code used for sine table

[BrianHG]

Well if you decide to forego the trig, just use the polynomial approximation.  It just requires the floating point multiply & divide.  But, once you already have the sin() elsewhere in your code, my 1 liner is the shortest solution.

[@rt]

I was out for a bit.
I have the computing time at startup since this only has to happen once to populate the table.

Program code size is most important,
and I don’t think it’s going to be beaten considering all I added was setting variables and a call to a function in a loop

I know what a sine is, how it relates to a circle, and that a computer generates a wave in the same manner it does anything else.
It has dawned on me I only need to learn what sin as a function means.

[magetoo]

IIRC the OPL (FM synthesis) chips from Yamaha used multi-pass delta encoding to reduce the storage needed to just one bit per point in the table.  There was an article where someone decapped a chip and reverse engineered the exact values in the sin and exp tables.

For a longer table that might be another option, not sure it makes sense here for just 91 values.

[@rt]

Is every variable in the BBC Micro a float like C64 BASIC?
I had to type it out like from a book!

Code: [Select]

// BBC MICRO
// GK DEC 2014
float DS = 0.0;
float DC = 0.0;
float mode = 1.0;
const float cfixed = 0.01;
// INITIAL CONDITIONS
float S = 0.0;
float C = 200.0;
while (x < 1250) { // COMPUTE SIN, COS
DS = C; DC = -S;
S = S + (DS*CONST);
C = C + (DC*CONST);
// PLOT CURVES
// X IS ITERATING
y = S + 500; setpixel();
y = C + 500; setpixel();
} // NEXT WHILE


[NivagSwerdna]

but now with a micro, as program memory gets short, I’d like to know how it could be generated in code if possible

It would help us out if you could tell us about the MCU environment you are working in.
Another thing which may or may not be significant is the nature of the data... your table is constant.  It therefore needs to end up in Flash not get allocated and initialised in RAM.  If you are lucky adding const will allow the compiler to work out it is immutable.

Here's an old page that describes the idea and even has sine as an example... http://hades.mech.northwestern.edu/index.php/Storing_constant_data_in_program_memory

This is mine and has served me well...

#define TABLE_SIZE 64

const int16 sineTable[TABLE_SIZE+1]=
{  0,
  13, 25, 38, 50, 63, 75, 87,100,112,124,136,148,160,172,184,196,
 207,218,230,241,252,263,273,284,294,304,314,324,334,343,352,361,
 370,379,387,395,403,410,418,425,432,438,445,451,456,462,467,472,
 477,481,485,489,492,496,499,501,503,505,507,509,510,510,511,511 };

[AndyC_772]

Is every variable in the BBC Micro a float like C64 BASIC?

In BBC BASIC, an integer is denoted by a % symbol at the end of the variable name, and a string by a $ symbol.

The "resident integer" variables are always stored in the same place in memory and so are quick to access; these are A% through Z%, and they all always exist. Any other variable is created automatically the first time it's assigned a value.

It's times like now I really miss my Archimedes 

[@rt]

Ok, some of those have to be made integers then, thanks!

With the lookup table in the first post, ideally I suppose it should be constant, but that’s a lot of memory.
Now though it’s potentially not constant in nature. It would be easy to multiply the frequency by at least some multiple, for example.

It’s a dsPic, MPLAB, and outdated C30. I think the C unsigned short can be declared constant (for a 16 bit datatype) if not WORD.

[tszaboo]

I would say any MCU will have more flash memory than RAM. So did you actually run out of the memory, or this is just a "just in case" nonsense.

[orolo]

Just for fun, some months ago I was toying with trig functions and devised a recursive method for the cosine function. The idea is that, for very small angles, cos(x) ~ 1 - x^2/2. So, to get the cosine of an angle, you divide that angle for a big power of two, until it's very small. Then you approximate the cosine using the upper formula. And now, you calculate the cosine of the double angle several times: cos(2x) = cos^2x - sin^2x = 2cos^2x - 1.

In python:

Code: [Select]

x = input("Angle in radians: ")
x0 = x / 2.0**10               # Divide by 2^10.
x0 = 1 - (x0**2 / 2)            # Initial approx.
for n in range(0, 10):         # Ten times.
    x0 = 2*x0**2 - 1           # Double the angle.
print x0

The result:

Code: [Select]

Angle in radians: 3.1416
-0.999999999963

To get the sine, shift the angle by -pi/2. It is a time-memory trade-off, and you can increase precision just increasing the number of iterations, depending on your float sizes. Well, maybe this is just a curiosity, but it was fun when I thought about it.

[hans]

91 words, or about 182 bytes, for a sine table is not a lot of memory.

Consider the cost of implementing a LUT generator or calling the function sin(). And all related soft-floating point libraries on the PIC. If this is the only location where you will use these functions, it ends up way more costly.
For instance, in an empty project calling up sin() on XC16 v1.30 increased compiled size from 216 bytes to 1388 bytes of FLASH usage.

Of which __sincosf added 210 bytes alone, the rest is supporting functions called by __sincosf and the floating point operators.

edit:

Additionally; I implemented a taylor sine series function using floating point:

Code: [Select]

float taylorSineSeries(float x)
{
    int8_t i;
    int8_t k;
    float val = 0;
    for (i = 1; i < 9; i++)
    {
        float v = x;
       
        // perform x^(2n) and (2n+1)!^-1
        for (k = 1; k < i; k++)
        {
            v = v*x*x/ (2*k) / (2*k+1);;
        }
       
        if ((i % 2) == 1)
            val += v;
        else
            val -= v;
    }
   
    return val;
}
XC16 v1.30 -O1 compiles this function down to 238 bytes, excluding floating point framework.
Then, to iterate the table:

Code: [Select]

    int8_t k;
    for (k = 0; k < 91; k++)
    {
        table[k] = 16384 * taylorSineSeries(k * 3.14159265359f/2/90);
    }

Total program size is now 1206 bytes instead of 1388.

[orolo]

As for the original question, how to generate the table, you can compute sin and cos in increments. Just precompute sin(pi/180) and cos(pi/180), and use sin(x + pi/180) and cos(x + pi/180) several times.

Code: [Select]

from math import *

# Step in 90 parts of pi/2:
x = 3.1415926535/(2*90)
sx = sin(x)
cx = cos(x)

# Initialize:
sa = 0
ca = 1

print sa,

for n in range(0, 90):
    t = sa*cx + ca*sx
    ca = ca*cx - sa*sx
    sa = t
    print int(round(sa * 16384)),

Result:

Code: [Select]

0 286 572 857 1143 1428 1713 1997 2280 2563
2845 3126 3406 3686 3964 4240 4516 4790 5063 5334
5604 5872 6138 6402 6664 6924 7182 7438 7692 7943
8192 8438 8682 8923 9162 9397 9630 9860 10087 10311
10531 10749 10963 11174 11381 11585 11786 11982 12176 12365
12551 12733 12911 13085 13255 13421 13583 13741 13894 14044
14189 14330 14466 14598 14726 14849 14968 15082 15191 15296
15396 15491 15582 15668 15749 15826 15897 15964 16026 16083
16135 16182 16225 16262 16294 16322 16344 16362 16374 16382
16384


[snarkysparky]

The following polynomial generates your table to within +/- 2 using double precision values.

With D in degrees between 0 and 90

result = 1.4585e-007*D^5  +  1.077529e-005*D^4  -   0.015136*D^3 + 0.015909*D^2 + 285.789*D

[BrianHG]

Ok, if you only need a sine table for synthesizing sounds, you do not need such accuracy.  Here is a comparison of my floating point sin(x) function compared to an all 16 bit integer generator using nothing but a multiply and a subtract  The integer generator should only take a few instructions in program space......

Code: [Select]

' http://www.freebasic.net/
' for dos, win & linux
' sintable comparison generator by BrianHG.

const pi = 3.141592654  : Rem This is the same thing >>> atn(1) * 4

const xscale  = 10
const xoffset = 50
const yscale  = 30
const yoffset = 700

Declare Sub plot_sin_table()

dim Shared as integer x,y,n,cr,cg,cb,tmp
dim Shared sin_table(0 to 90) as integer

ScreenRes 1024,768,32,0,0


' Generate sine_table using the function sin(x).
   for n=0 to 90
    sin_table(n) = sin(pi*(n/90)/2) * 16384
   next n

' Plot sin_table(n) in green
cr=0:cg=255:cb=0:plot_sin_table()


' approximate sine_table using all 16 bit integer polynomile.
   for n=0 to 90
    tmp = (90-n) * 2
    sin_table(n) = (32580 - tmp^2) / 2 : REM change the divide by 2 to a right shift by 1
   next n

' Plot sin_table(n) in red
cr=255:cg=0:cb=0:plot_sin_table()


' Wait until a key if pressed before ending.
  while (inkey$="") : sleep 1 : wend
END

Sub plot_sin_table()
  for n=0 to 89
   line (n*xscale+xoffset, yoffset-(sin_table(n)/yscale)) - ((n+1)*xscale+xoffset, yoffset-(sin_table(n+1)/yscale) ), RGB(cr,cg,cb)
  next n
End Sub
See attached image for output:

[BrianHG]

I've added a better approximation in yellow, but, to compute it, it requires the math in 32bit and it has a divide:
Note that no matter what you do, since the dsPIC has a hardware multiplier, the red trace will take up the absolute least amount of code memory and will also execute in the least amount of time by a long shot.  In fact, the red trace should be calculated fast enough to be done on the fly...

Code: [Select]

' http://www.freebasic.net/
' for dos, win & linux
' sintable comparison generator by BrianHG.

const pi = 3.141592654  : Rem This is the same thing >>> atn(1) * 4

const xscale  = 10
const xoffset = 50
const yscale  = 30
const yoffset = 700

Declare Sub plot_sin_table()

dim Shared as integer x,y,n,cr,cg,cb,tmp
dim Shared sin_table(0 to 90) as integer

ScreenRes 1024,768,32,0,0


' Generate sine_table using the function sin(x).
   for n=0 to 90
    sin_table(n) = sin(pi*(n/90)/2) * 16384
   next n

' Plot sin_table(n) in green
cr=0:cg=255:cb=0:plot_sin_table()


' approximate sine_table using all 16 bit integer polynomile.
   for n=0 to 90
    tmp = (90-n) * 2
    sin_table(n) = (32580 - tmp^2) / 2 : REM change the divide by 2 to a right shift by 1
   next n

' Plot sin_table(n) in red
cr=255:cg=0:cb=0:plot_sin_table()


' slightly better approximation of sine_table using 32 bit integer polynomile with divide.
   for n=0 to 90
    tmp = (95-n) * 6
    sin_table(n) = (326619 - tmp^2) / 20 : REM change the divide by 2 to a right shift by 1
   next n

' Plot sin_table(n) in yellow
cr=255:cg=255:cb=0:plot_sin_table()


' Wait until a key if pressed before ending.
  while (inkey$="") : sleep 1 : wend
END

Sub plot_sin_table()
  for n=0 to 89
   line (n*xscale+xoffset, yoffset-(sin_table(n)/yscale)) - ((n+1)*xscale+xoffset, yoffset-(sin_table(n+1)/yscale) ), RGB(cr,cg,cb)
  next n
End Sub
Now, it is possible to perform direct full sine table approximations instead of 1/4 using with positive and negative coordinates, but, it will require a little more thought and it exceeds the construction of the 1/4 table you asked for.  For one thing, getting rid of that divide by 2 amplitude adjustment preserves the polarity MSB bit in the 16 bit result.

[BrianHG]

And, for the final one.  This uses nothing but floating point ^1.73 and a multiply by 6.82.  Everything else is the same.  This approximates a sine wave really, really close...  See attached code and plot:

Code: [Select]

' http://www.freebasic.net/
' for dos, win & linux
' sintable comparison generator by BrianHG.

const pi = 3.141592654  : Rem This is the same thing >>> atn(1) * 4

const xscale  = 10
const xoffset = 50
const yscale  = 30
const yoffset = 720

Declare Sub plot_sin_table()

dim Shared as integer x,y,n,cr,cg,cb,tmp
dim Shared sin_table(0 to 90) as integer

ScreenRes 1152,768,32,0,0


' Generate sine_table using the function sin(x).
   for n=0 to 90
    sin_table(n) = sin(pi*(n/90)/2) * 16384
   next n

' Plot sin_table(n) in green
cr=0:cg=255:cb=0:plot_sin_table()


' approximate sine_table using floating point ^1.73 and a multiply by 6.82
   for n=0 to 90
    tmp = (90-n)
    sin_table(n) = (2403.475 - tmp^1.73) * 6.82
   next n

' Plot sin_table(n) in red
cr=255:cg=0:cb=0:plot_sin_table()


' Wait until a key if pressed before ending.
  while (inkey$="") : sleep 1 : wend
END

Sub plot_sin_table()
  for n=0 to 90
   ? sin_table(n),
  next n
   ?:?

  for n=0 to 89
   line (n*xscale+xoffset, yoffset-(sin_table(n)/yscale)) - ((n+1)*xscale+xoffset, yoffset-(sin_table(n+1)/yscale) ), RGB(cr,cg,cb)
  next n
End Sub


[@rt]

But, once you already have the sin() elsewhere in your code, my 1 liner is the shortest solution.

If you mean this:

Code: [Select]

n = 0;
while (n < 91) {
sin_table[n] = sin(pi*(n/90)/2) * 16384;
n++;
}

is smaller than this:

Code: [Select]

n = 0;
while (n < 91) {
x = 16384; y = 0;
pixelrotation();
sin_table[n] = y;
rotation++;
n++;
}

I have doubt, but will find out for sure tonight.



Yes, it’s really full. This was before I started trimming it down at all.

[orolo]

Using only 16-bit integer arithmetic, you can recreate the table with relative error less than 0.001 with this:

Code: [Select]

#include <stdio.h>
int main()
{
   unsigned short sx, cx, sa, ca, t;
   char c;

   sx = 0b0000010001111000;
   cx = 0b1111111111110110;
   sa = 0;
   ca = 0b1111111111111111;
   
   for (c=0 ; c < 90 ; c++) {
    t  = ((long) (sa*cx + ca*sx))>>16;
    ca = ((long) (ca*cx - sa*sx))>>16;
    sa = t;
    t = (sa >> 2) + ((sa&2)>>1);
    printf("%d ", t);
   }
}
The precision could be improved further than 16 bits working with assembler and taking care of the multiplication overflows, I guess.

The output is:

Code: [Select]

286 572 857 1142 1427 1712 1996 2279 2562 2844
3125 3405 3684 3962 4239 4514 4788 5061 5332 5601
5869 6135 6399 6661 6921 7179 7435 7688 7940 8188
8435 8678 8919 9158 9393 9626 9856 10082 10306 10527
10744 10958 11169 11376 11580 11780 11977 12170 12359 12544
12726 12904 13078 13248 13414 13575 13733 13887 14036 14181
14321 14458 14589 14717 14840 14958 15072 15181 15286 15386
15481 15571 15657 15738 15814 15886 15952 16014 16071 16123
16170 16212 16249 16281 16308 16331 16348 16360 16367 16370

However, if you already have a function pixelrotation which computes the sine, this code will probably be inferior.

[BrianHG]

Comparing with Orolo's code, his sine wave is almost perfect, remembering to insert a 0 at the beginning of the table.  Here is the basic conversion with a plot:

Code: [Select]

' http://www.freebasic.net/
' for dos, win & linux
' sintable comparison generator by BrianHG.

const pi = 3.141592654  : Rem This is the same thing >>> atn(1) * 4

const xscale  = 10
const xoffset = 50
const yscale  = 30
const yoffset = 720

Declare Sub plot_sin_table()

dim Shared as unsigned integer x,y,n,cr,cg,cb,tmp,sx,cx,sa,ca,t,c
dim Shared sin_table(0 to 90) as integer

ScreenRes 1152,768,32,0,0


' Generate sine_table using the function sin(x).
   for n=0 to 90
    sin_table(n) = sin(pi*(n/90)/2) * 16384
   next n

' Plot sin_table(n) in green
cr=0:cg=255:cb=0:plot_sin_table()


' approximate sine_table 'orolo's code.
' please excuse basic's limitation of shifting bits, I had to replace them with divide...
' It is expected that in the CPU, you would use the bit shifting for speed and code size.

   sx = 1144  ' 0b0000010001111000;
   cx = 65526 ' 0b1111111111110110;
   sa = 0     ' 0;
   ca = 65535 ' 0b1111111111111111;
   
    sin_table(0) = 0
   for c=1 to 90
    t  = (sa*cx + ca*sx)/65536
    ca = (ca*cx - sa*sx)/65536
    sa = t
    t = (sa / 4) + ((sa and 2)/2)
    sin_table(c) = t
   next c

' Plot sin_table(n) in cyan
cr=0:cg=255:cb=255:plot_sin_table()


' Wait until a key if pressed before ending.
  while (inkey$="") : sleep 1 : wend
END

Sub plot_sin_table()
  for n=0 to 90
   ? sin_table(n),
  next n
   ?:?

  for n=0 to 89
   line (n*xscale+xoffset, yoffset-(sin_table(n)/yscale)) - ((n+1)*xscale+xoffset, yoffset-(sin_table(n+1)/yscale) ), RGB(cr,cg,cb)
  next n
End Sub
As you can see in the attached image, Orolo's code's cyan color output sits right on top of the original green sine wave with only a pixel deviation or 2 here and there.  This is the best I've seen yet requiring very little processing power to generate a table.
If you don't want any trig functions with floats, use Orolo's code.  No one will know the difference.
If you dont want to store the table and want to compute on the fly, quick random access the first ' (8100-(90-n)^2)>>1 ' trick I made will give you a quick 1 multiply, 2 subtracts + bit shift direct access to a curve, but nothing as accurate as Orolo's waveform.

[@rt]

Yeah I don’t mean to imply anything is inferior or not, some of the content I’m reading here would be better to learn something,
just not when the goal at the moment is size.
But yes I need to be trying to understand more.
A project I have in mind does involve frequency measurements perhaps with FFT to look at superposition (reflected waves) as well if possible.
So something similar to how I suppose radar works, looking at both he transmitted and reflected frequencies to break them down to determine distance or location.

This is just for the likes of sine wave text and plasma graphics though.

[JPortici]

Using only 16-bit integer arithmetic, you can recreate the table with relative error less than 0.001 with this:

 
can you point to the math behind this?
beautiful!
i sort of had this problem (i wanted to shrink down a program that had a sine wave table in prog mem so i could use a smaller dspic)

[mikeselectricstuff]

Using only 16-bit integer arithmetic

    t  = ((long)......

..well not all 16-bit, but neat nonetheless

[BrianHG]

Yeah I don’t mean to imply anything is inferior or not, some of the content I’m reading here would be better to learn something,
just not when the goal at the moment is size.
But yes I need to be trying to understand more.
A project I have in mind does involve frequency measurements perhaps with FFT to look at superposition (reflected waves) as well if possible.
So something similar to how I suppose radar works, looking at both he transmitted and reflected frequencies to break them down to determine distance or location.

This is just for the likes of sine wave text and plasma graphics though.

For synthesis, and FFT analysis, of complex waveforms, this mechanical computer does help visualize how to do so:
Part #1 of 6, introduction: (the additional parts are linked at the end of the video)

Enjoy.

[BrianHG]

Using only 16-bit integer arithmetic

    t  = ((long)......

..well not all 16-bit, but neat nonetheless

For those first 2, 1 addition, 1 subtraction, he only requires the upper 16bits of the multiplies summed together before shifting down the results.  If you are using a 16bit dsPICs & it's hardware multiplier's upper 16bit results only within your calculations, basically coded like this:  (Taking advantage of the dsPic's op-code like this might means some lines in assembly, but you will generate the sine ridiculously fast with only around 20 lines of OP-code memory for everything...)

Code: [Select]

    sin_table(0) = 0
   for c=1 to 90
    t  = (int(sa*cx/65536) + int(ca*sx/65536))
    ca = (int(ca*cx/65536) - int(sa*sx/65536))
    sa = t
    t = (sa / 4) + ((sa and 2)/2)
    sin_table(c) = t
   next c
More error is introduced into Orolo's generator.  As you can see in the attached image, his Cyan trace no longer sits right on top, masking out the original green one, but, it is still damn amazingly close for only 20 instructions in rom to generate:  (You can create Orolo's higher quality version I plotted above with a 32bit PIC in C since the assembly will translate almost identically if you write everything with longs.)

[BrianHG]

Using only 16-bit integer arithmetic, you can recreate the table with relative error less than 0.001 with this:

 
can you point to the math behind this?
beautiful!
i sort of had this problem (i wanted to shrink down a program that had a sine wave table in prog mem so i could use a smaller dspic)

Orolo is manually drawing a sine and cosine just like the first basic program post by GK -> https://www.eevblog.com/forum/microcontrollers/code-used-for-sine-table/msg1107610/#msg1107610

The difference is instead of using floating point, he's using 16 bit whole numbers, when multiplying, using the upper 16 bits as the whole number part, and the lower 16 bits as the floating number part, then he's shifting over by 16 bits holding only the fractional part.

The starting parameters are just so that the we get the beginning of a sine and sets it's amplitude.  The maximum numbers were chosen for the best possible quality sine.
The final ' t = (sa >> 2) + ((sa&2)>>1); ' just scales down the results so that the table goes from 0-16384 instead of 0-65535.

Orolo, did I get this right?

[JPortici]

oh yeah.. i missed the images in that reply (saw before he added them. and of course the code was posted after but didn't pay too much attention to it, shame on me.)