Code used for sine table

[orolo]

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 differential equation is a great idea, but I wasn't doing that. To put it briefly, it took the complex number z = cos(1) + I*sin(1), and computed its powers, 1, z, z^2, z^3,... each time multiplying, in the complex plane, by z. Since z^n = cos(n) + I*sin(n), you get the sine (and cosine) table.

In binary arithmetic, sin(1) = 0,0000010001111000... and cos(1)=0,1111111111110110... . I used the fact that 1 = 0,11111111... in binary, also.

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.

Exactly. The ((sa&2)>>1) part adds rounding to the scaling down. Just using t = sa>>2 is more compact, but loses a bit of precision.

Thank you for the interest and the graphical analysis, that was great. 

Edit:

Just to keep  , if a little increase in complexity is allowed, the table can be computed with more precision and half the number of iterations. The idea is using the cosines to compute half the table of sines:

Code: [Select]

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

   sx = 0b0000010001111000;
   cx = 0b1111111111110110;
   sa = 0;
   ca = 0b1111111111111111;
   
   for (c=0 ; c < 46 ; c++) {
    t = (sa >> 2) + ((sa&2)>>1);
    t2 = (ca >> 2) + ((ca&2)>>1);
    printf("%d %d\n", t, t2);
   
    t  = ((long) (sa*cx + ca*sx))>>16;
    ca = ((long) (ca*cx - sa*sx))>>16;
    sa = t;
   }
}
The idea is to fill the lower half of the table with t, and the upper part with t2, in descending order. With this trick, the precision of the upper part of the table increases remarkably. Both parts of the table meet in the middle with the same value. And this time the table starts at 0! (thanks, Brian).

[BrianHG]

Here you go Orolo,  Green is the original sin(x) command.  The Cyan is your new code with the (long) addition/subtraction between the 2 multiplies.  Red is your new code adding/subtracting only the upper 16 bits from the 2 multiplies as this coding for 16 bit CPU with a hardware multiplier with the upper 16 bit results would be the most compact means generating this sine table.

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 = 800

Declare Sub plot_sin_table()

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

ScreenRes 1152,816,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 NEW code using Long.
' 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;
   
   for c=0 to 46

    t  = int(sa / 4) + ((sa and 2)/2)
    t2 = int(ca / 4) + ((ca and 2)/2)

    sin_table( c    ) = t
    sin_table( 90-c ) = t2

    t  = ((sa*cx) + (ca*sx))/65536
    ca = ((ca*cx) - (sa*sx))/65536
    sa = t

   next c

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


' approximate sine_table 'orolo's NEW code using Integer multiply.
' 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;
   
   for c=0 to 46

    t  = int(sa / 4) + ((sa and 2)/2)
    t2 = int(ca / 4) + ((ca and 2)/2)

    sin_table( c    ) = t
    sin_table( 90-c ) = t2

    t  = int((sa*cx)/65536) + int((ca*sx)/65536)
    ca = int((ca*cx)/65536) - int((sa*sx)/65536)
    sa = t

   next c

' 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()
color RGB(cr,cg,cb)
  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 resulting image, you will now need to download and zoom into the image, the cheap integer version, now red, has drastically improved over my previous integer version which had green segments showing through the cyan throughout the upper half of the sine wave.

[BrianHG]

If you mean this:

Code: [Select]

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


I believe that the C compiler should simplify this one out, but, just in case it doesn't:
--------------------------------------------------------------
sin_table[n] = sin(pi*(n/90)/2) * 16384;
===is equivalent to===
sin_table[n] = sin( n*0.01745329252 ) * 16384;
--------------------------------------------------------------

[jesuscf]

One possible solution is to store the second differences.   After some quick spread sheet calculations, all the second differences fit in only three bits.  If they are stored in 32-bit constants, 10 3-bit differences can be stored per 32-bit constant.  That reduces the size of the original table from 182 bytes to just 36 bytes:

Code: [Select]

unsigned long packed_differences[] ={
   0x12451288, 0x142da492, 0x1551c4db,
   0x2c764724, 0x26965af3, 0x2dd6db6d,
   0x2ebb5bb5, 0x3dbf5db6, 0x05f7df76
}
Here is the program I used to test the whole idea.  It both generates the compressed table of second differences and un-compresses it to produce the original table.

Code: [Select]

#define WORD short int
#define TSIZE 91

WORD sin_table[TSIZE] ={
    0,   285,   571,   857,  1142,  1427,  1712,  1996,  2280,  2563,
    2845,  3126,  3406,  3685,  3963,  4240,  4516,  4790,  5062,  5334,
    5603,  5871,  6137,  6401,  6663,  6924,  7182,  7438,  7691,  7943,
    8191,  8438,  8682,  8923,  9161,  9397,  9630,  9860, 10086, 10310,
    10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
    12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
    14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
    15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
    16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
    16384
};

WORD diff1[91];
WORD diff2[91];

long pack_diff2[10];
int y=0, d1=285, d2;

void main (void)
{
int j, k;

// Compute first differences
for(j=0; j<(TSIZE-1); j++)
{
diff1[j]=sin_table[j+1]-sin_table[j];
}

// Compute second differences
for(j=0; j<(TSIZE-2); j++)
{
diff2[j]=diff1[j]-diff1[j+1]+1; // Change the order and add one so differences are all positive
}

// Pack the differences: each second difference uses 3 bits.
printf("unsigned long packed_differences[] ={\n");
for(j=0; j<9; j++)
{
pack_diff2[j]=0;
for(k=0; k<10; k++)
{
pack_diff2[j]<<=3;
pack_diff2[j]|=(diff2[(j*10)+(9-k)]&0x07);
}
printf("   0x%08x,\n", pack_diff2[j]);
}
printf("};\n");

// With the packed differences, rebuilt the original table:
printf("Un-Packed table:\n");
for(j=0; j<9; j++)
{
for(k=0; k<10; k++)
{
printf("%5d ", y);
y+=d1;
d2=pack_diff2[j]&0x7;
pack_diff2[j]>>=3;
d1=d1-d2+1;
}
printf("\n");
}
printf("%d\n", y);
}
This is the rebuilt sine table from the packed table of second differences:

Code: [Select]

    0   285   571   857  1142  1427  1712  1996  2280  2563
 2845  3126  3406  3685  3963  4240  4516  4790  5062  5334
 5603  5871  6137  6401  6663  6924  7182  7438  7691  7943
 8191  8438  8682  8923  9161  9397  9630  9860 10086 10310
10531 10748 10963 11173 11381 11585 11785 11982 12175 12365
12550 12732 12910 13084 13254 13420 13582 13740 13894 14043
14188 14329 14466 14598 14725 14848 14967 15081 15190 15295
15395 15491 15582 15668 15749 15825 15897 15964 16025 16082
16135 16182 16224 16261 16294 16321 16344 16361 16374 16381
16384
This is a complete test program using the compressed table of second differences.  It generates the original table and prints it.  I am sure it can be optimized a bit more...

Code: [Select]

#define WORD short int
#define TSIZE 91

unsigned long pack_diff2[9] = {
   0x12451288, 0x142da492, 0x1551c4db,
   0x2c764724, 0x26965af3, 0x2dd6db6d,
   0x2ebb5bb5, 0x3dbf5db6, 0x05f7df76
};

WORD sin_table[TSIZE];

void main (void)
{
WORD y=0, d1=285;
unsigned long d2;
char i=0, j, k;

// Unpack the table
for(j=0; j<9; j++)
{
d2=pack_diff2[j];
for(k=0; k<10; k++)
{
sin_table[i++]=y;
y+=d1;
d1=d1-(d2&0x7)+1;
d2>>=3;
}
}
sin_table[i]=y;

// Just to make sure, print the table
for(i=0; i<TSIZE; i++)
{
printf("%5d ", sin_table[i]);
if((i%10)==0) printf("\n");
}
printf("\n");
}
Jesus

[jesuscf]

Solving these two coupled differential equations will do it.

dsine/dt = cos
dcos/dt = -sine

I'm not a mathematician either and I'm not necessarily saying that this is an efficient solution, and the code example is in BASIC rather than C. This is how I solved the problem on one lazy day just out of idle curiosity on my Beeb. I can't figure an integration algorithm stripped any barer.

The way the solution is implemented, using Forward Euler integration, produces results that increase in magnitude over time (FE.jpg).  This is because Forward Euler integration introduces "negative" damping.  One can use the Backward Euler integration, but then the opposite problem arises as Backward Euler introduces positive damping and produces results that decrease in magnitude over time (BE.jpg) .  The solution to this problem is to use trapezoidal integration which does not introduce any damping (trap.jpg).  The resulting equations using trapezoidal integration are

f1(x)=2*h*f2(x-h)+f1(x-2h)
f2(x)=-2*h*f1(x-h)+f2(x-2h)

where h is the time step, f1(x) is the sine at time x, f1(x-h) is the sine at the previous time step, and f1(x-2h) is the sine at the time step before that.  Similarly f2(x) is the cosine at time x, f2(x-h) is the cosine at the previous time step, and f2(x-2h) is the cosine before that.  table_trap.jpg shows the numeric results obtained with the equations above and a time step of one degree (pi()/180=0.01745).

[CatalinaWOW]

Lots of good answers here.

For this and all other problems relating to digital implementations of math, one of my first stops is the Handbook of Mathematical Functions by Abramowicz and Stegun.  It has formal definitions, series approximations, rational approximations and tabular data for most functions anyone will encounter.  It is available online from several sources, one of which is here:

http://people.math.sfu.ca/~cbm/aands/

[@rt]

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.

Something I thought about more. Is this a formal Computer Science assertion?
Thinking about it more, if you can’t trust the content of RAM, you can’t trust any software.
You can’t even trust the program counter, so what’s the difference?

[jesuscf]

I have some results.  Using a table of packed second differences to generate the sine table resulted in some code size reduction.  Compiled for an 8051, I wrote the unpacking routine in assembly and changed the table of second differences from long to int.  It takes 108 bytes of code memory to initialize the table in expanded RAM.  Here is the program:

Code: [Select]

#define TSIZE 91

const int pack_diff2[18] = {
0x1288, 0x248a, 0x2492, 0x285b, 0x44db, 0x2aa3,
0x4724, 0x58ec, 0x5af3, 0x4d2c, 0x5b6d, 0x5bad,
0x5bb5, 0x5d76, 0x5db6, 0x7b7e, 0x5f76, 0x0bef
};

xdata unsigned int sin_table[TSIZE];

void expand_table (void) _naked
{
_asm
mov r4, #18 ; loop counter 1
mov dptr, #_pack_diff2
MOV P2, #(_sin_table/0x100) ;  Use _EMI0CN in C8051F38x instead
MOV R0,#(_sin_table%0x100)
; First entry in table is zero thanks to crt0 (skip setting it to zero)
;clr a
;movx @R0, a
;inc R0
;movx @R0, a
;dec R0
; Initialize first D1
mov r2, #(285%0x100)
mov r3, #(285/0x100)
L1:
; Get packed 3-bit D2 x 5 from compressed table.  Store in [r7,r6]
clr a
movc a, @a+dptr
mov r6, a
inc dptr
clr a
movc a, @a+dptr
mov r7, a
inc dptr

mov r5, #5 ;  loop counter 2
L2:
; add d1 to previous value and save in xram table
mov a, r0
add a,#2
mov r1, a

movx a, @R0
add a, r2
movx @r1, a
inc R0
inc r1
movx a, @R0
addc a, r3
movx @r1, a
inc r0

; d1=d1+1-(d2&0x7);
mov a, r6
anl a, #7
mov r1, a

clr c
mov a, r2
subb a, r1
mov r2, a
mov a, r3
subb a, #0
mov r3, a

inc r2
cjne r2, #0, L3
inc r3
L3:
; Shift compressed D2 x 5 3-bits right:
mov r1, #3
L4:
mov a, r7
rrc a
mov r7, a
mov a, r6
rrc a
mov r6, a
djnz r1, L4

djnz r5, L2
djnz r4, L1
ret
_endasm;
}

void main (void)
{
expand_table();
}

The program with the whole table uses 183 bytes of code memory (there is an extra 'ret' instruction that takes one byte added to the main program):

Code: [Select]

#define WORD short int
#define TSIZE 91

const WORD sin_table[TSIZE] ={
    0,   285,   571,   857,  1142,  1427,  1712,  1996,  2280,  2563,
    2845,  3126,  3406,  3685,  3963,  4240,  4516,  4790,  5062,  5334,
    5603,  5871,  6137,  6401,  6663,  6924,  7182,  7438,  7691,  7943,
    8191,  8438,  8682,  8923,  9161,  9397,  9630,  9860, 10086, 10310,
    10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
    12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
    14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
    15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
    16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
    16384
};

void main (void)
{
}

Using a table of second differences to generate the sine table saves 75 bytes of code memory in the 8051.

Jesus

[jesuscf]

Attached is the Excel spread sheet I used to test some of the ideas presented in this thread.

[Krischu]

Great post! I so much liked this 4-liner that I immediately tried it out in a Jupyter (python) notebook:

Code: [Select]

%pylab inline

from numpy import *
import numpy as np

pi=math.pi
arsiz=1260

x = np.linspace(0, 2*pi, arsiz)

asin=np.arange(arsiz,dtype = float)
acos=np.arange(arsiz,dtype = float)
const=0.005
s=0
c=255
for i in range(arsiz):
    ds = c
    dc = -s
    s=s+ds*const
    c=c+dc*const
    asin[i]=s
    acos[i]=c
pylab.plot(x, asin,x,acos)
for i in range(arsiz):
    print(asin[i])       

Now, what I would like to understand:

What values influence the

a) frequency
b) amplitude
c) duration

of the generated function values?

[RAPo]

Yes, that would be my suggestion also: use the cordic algoritm
See

In 1974 a 4-function plus trig functions calculator was implemented in 320 instructions. Total. Including i/o operations.

You might find inspiration by studying how it was done: http://files.righto.com/calculator/sinclair_scientific_simulator.html

Or there's always CORDIC, of course.

[bson]

What kind of microcontroller is being discussed?  Usually the exact opposite is desirable, as flash is vastly cheaper than SRAM and there's a whole lot more of it.  As in a megabyte, vs 128kB or so of SRAM.

[peter-h]

When I did a waveform generator with the 32F417, I used this

Code: [Select]



// Compact sin(x) function which doesn't need math.h
// Uses Horner's method
// i is from 0 to SPC-1

#define PI 3.14159265358979323846
#define SCALE 13107

int sin2(int i)
{
  double rad = 2*PI*i/SPC, term=rad, sin=rad, old_sin=0;
  //for (int n=3; sin != old_sin; n+=2)
  for (int n=3; (sin != old_sin) && (n < 30); n+=2)
  {
    old_sin = sin;
    term *= (rad*rad) / (-n * (n-1));
    sin += term;
  }
  return (sin * SCALE) + (sin >= 0 ? 0.5 : -0.5);
}

As noted, the reason for not just using trig functions is to avoid the need for math.h.

Then I use this to fill the table. I am using the two 12-bit DACs.

Code: [Select]


#define SPC 128 // number of samples per cycle - 32/64/128/256 (see comments in .c)
#define SHIFT 4 // # of samples to shift sine table to the left (# is related to SPC)

/*
* Precompute a sin(x) table, for a whole cycle (0 to 2*PI) with signed 16 bit int
* values from -13107 to +13107. Uses int16 rather than int32 to save RAM.
* 13107 is 16384 * 2 / 2.5, to get DAC swing from +0.25V to +2.25V (2V P-P)

* The table is shifted to the left by SHIFT samples, to compensate for the 2-sample
* delay in the wave gen, plus the phase delay in the DAC output filter (probably
* arranged to be equal to 1 sample period)
*
*/

/* This version needs #include "math.h" which adds about 5k to the FLASH usage
double xinc=2*3.1415926/SPC;
double x=xinc*SHIFT;

for (int i=0;i<SPC;i++)
{
sintab[i] = sin(x)*13107.0;
x+=xinc;
if (x>(2*3.1415926)) x=0;
}
*/

// This uses the simple sin2(x) algorithm

for (int i=0;i<SPC;i++)
{
sintab[i] = sin2(i+SHIFT);
}

// Set DAC tables to midpoint (2048) in case they contain duff data

for (int i=0;i<SPC;i++)
{
wavebuf0[i]=2048|(2048<<16);
wavebuf1[i]=2048|(2048<<16);
}


[zapta]

I assume if I begin with just the empty 91 word array, it would be possible, with the right formula, to generate the
half cycle at run time just once, and continue to use it.
Any ideas appreciated!

Code: [Select]

WORD sin_table[91] ={
    0,   285,   571,   857,  1142,  1427,  1712,  1996,  2280,  2563,
    2845,  3126,  3406,  3685,  3963,  4240,  4516,  4790,  5062,  5334,
    5603,  5871,  6137,  6401,  6663,  6924,  7182,  7438,  7691,  7943,
    8191,  8438,  8682,  8923,  9161,  9397,  9630,  9860, 10086, 10310,
    10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
    12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
    14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
    15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
    16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
    16384
};


Am I looking at these numbers wrong, this table appears to contain the first 1/4 cycle of a sine wave, not half cycle.

How much lower in resolution is acceptable for this table?
Have you considered storing this table in the processor's flash mem, or external I2C prom, then loading it into ram for use?
Could you get away with an inverted '((x/xscale)^2)/amplitudefactor' 1/4 sine wave approximation?  Note that with 8 bit micros, just the 32 bit integer multiply & divide algorithms introduced by your C compiler may be longer than the 91 words of rom.

Also depending on compiler, defining a constant may eat twice or quadruple as much rom instruction words compared to storing the table as raw data table in free rom space, as long as your MCU can read op-code as data for those MCUs which have the feature.  You will need to create a table-read loop to read the flash rom into ram yourself, but, you will cut the instruction mem size of the table by as much as 4 fold in some cases, + the instructions of the table-read loop to copy that rom space into usable ram.