Author Topic: Code used for sine table  (Read 19795 times)

0 Members and 1 Guest are viewing this topic.

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Code used for sine table
« on: January 06, 2017, 10:35:25 am »
Hi Guys :)
I have an explicitly defined sine table in C that I’ve been using for years,
but now with a micro, as program memory gets short, I’d like to know how it could be generated in code if possible
....or pretty close to it.

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
};
 

Offline alexanderbrevig

  • Frequent Contributor
  • **
  • Posts: 700
  • Country: no
  • Musician, developer and EE hobbyist
    • alexanderbrevig.com
Re: Code used for sine table
« Reply #1 on: January 06, 2017, 10:38:16 am »
You could maybe move it from program memory to RAM, or more likely from RAM to program memory (guessing by the lack of qualifiers that you use RAM now).
For approximating with code I've had luck with Taylor series.
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #2 on: January 06, 2017, 10:45:49 am »
If I defined it as constant I presume that would save RAM, but RAM isn’t the issue.... I’m not short on it yet.
Because it’s explicitly defined this way, all 91 values have to first be stored in program memory somewhere before being loaded to RAM.
 

Offline alexanderbrevig

  • Frequent Contributor
  • **
  • Posts: 700
  • Country: no
  • Musician, developer and EE hobbyist
    • alexanderbrevig.com
Re: Code used for sine table
« Reply #3 on: January 06, 2017, 10:48:57 am »
Because it’s explicitly defined this way, all 91 values have to first be stored in program memory somewhere before being loaded to RAM.
:palm: of course you're right. Maybe the Taylor series was useful for you then?
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #4 on: January 06, 2017, 10:55:36 am »
I’m not familiar with it.. don’t claim to be a math guru of course.

I do have functions I’ve used to generate audio at given sample rates and frequencies which could do it,
but have no idea how to land it near the values I want.

It would look something like looping through this:

Code: [Select]
   const float PI2 = 3.14159265359f * 2;
    double x = PI2 * frequency * rate;
    buffer[index] = sin(x) * amplitude;

I’m thinking maybe run the audio on a big platform with debug,
and maybe play with it and see if I can get a table close.

 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #5 on: January 06, 2017, 10:59:58 am »

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.
« Last Edit: January 06, 2017, 11:16:25 am by BrianHG »
 

Offline tggzzz

  • Super Contributor
  • ***
  • Posts: 19493
  • Country: gb
  • Numbers, not adjectives
    • Having fun doing more, with less
Re: Code used for sine table
« Reply #6 on: January 06, 2017, 11:24:10 am »
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.
There are lies, damned lies, statistics - and ADC/DAC specs.
Glider pilot's aphorism: "there is no substitute for span". Retort: "There is a substitute: skill+imagination. But you can buy span".
Having fun doing more, with less
 
The following users thanked this post: albert22

Offline iaeen

  • Regular Contributor
  • *
  • Posts: 65
  • Country: us
Re: Code used for sine table
« Reply #7 on: January 06, 2017, 11:28:27 am »
Yeah, a half cycle should return to 0 at 180 degrees. This looks like 90 degrees which is a quarter cycle. It's still enough to get the value of sine for any angle though...

I do believe that Taylor expansions are used by modern calculators/computers to get these values, so that is what you should look into. On the other hand, calculators need to be able to return a whole range of functions (all the trig functions, exp, log, etc...), so they can't afford to waste memory on tables. If you don't have that limitation, it might make sense to use a table if your processor can retrieve the value from memory faster than it can do the calculation.
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #8 on: January 06, 2017, 11:32:54 am »
You could call it a 1/4 wave because it’s all positive, but it’s indexed from 0-90,
and then the indexing reverses direction and reads the array backward,
so in the end, if used for sine text or plasma, etc, you end up with a sine wave.

I presume to originally save memory because it’s many years old.


I’ve sort of had luck with my audio function... kinda :D I think maybe using the wrong variable types.

The micro is dsPic with software FP. All of the trig is already in code, so introducing trig calls shouldn’t consume more memory.
 

Offline GK

  • Super Contributor
  • ***
  • Posts: 2607
  • Country: au
Re: Code used for sine table
« Reply #9 on: January 06, 2017, 11:33:15 am »
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.




« Last Edit: January 06, 2017, 11:36:58 am by GK »
Bzzzzt. No longer care, over this forum shit.........ZZzzzzzzzzzzzzzzz
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #10 on: January 06, 2017, 11:36:56 am »
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.

Wow... Reminds me of doing everything in early PIC assembly in the 90s.
Today with C programming, I think just an initialize for-loop eats more than that.
Without knowing the CPU being used and which compiler, I cant say doing the math might be the right solution.  As the question was given, with a smaller 8 bit PIC, I would choose a processor with user flash eerom and store the sine table there.  It appears every bit of program rom space is being used up already.  This probably means floating point is out of the question, and maybe even divide.
-------------------------------------------------------------------
Oops, updates were given while I was typing....
-------------------------------------------------------------------
If you are using a dspic, with floats already in use, use GK's solution, If I am reading my basic right...

« Last Edit: January 06, 2017, 11:43:35 am by BrianHG »
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #11 on: January 06, 2017, 11:42:26 am »
What’s the display resolution?

You probably wouldn’t believe I wrote this from scratch by scaling the audio generator.
If I can dig up any of that I’m set.

 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #12 on: January 06, 2017, 11:47:11 am »
Looking at GK's code, changing the CONST value will scale the amplitude and definition of the sine wave.  Also, the range on the 'for x' will change as well.
 

Offline GK

  • Super Contributor
  • ***
  • Posts: 2607
  • Country: au
Re: Code used for sine table
« Reply #13 on: January 06, 2017, 11:52:13 am »
Yes, the smaller the CONST variable, the finer the effective integration step. Good luck implementing this in code that takes up less room than a pi/2 look-up table though.
 
Bzzzzt. No longer care, over this forum shit.........ZZzzzzzzzzzzzzzzz
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #14 on: January 06, 2017, 11:55:18 am »
To get the resolution of the video, with such a low res fixed sine table, you will need to interpolate/up sample the fixed table.
GKs basic code should be able to synth a really high res table, but, since it is progressively computed, if you want rapid random point access, you will need to render and store at least the first 1/4 of the sine in a ram table.
Summing up sine waves from different points in the table to synth a complex waveform with multiple frequencies and amplitudes becomes easy from that point on.
« Last Edit: January 06, 2017, 12:01:58 pm by BrianHG »
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #15 on: January 06, 2017, 11:59:30 am »
Which DSPic will you be using?  Many have user eeproms for efficiently storing your table.  In fact, you can make a huge table with the 1k or 2k devices which can be read directly as needed, instead of eating up CPU ram.

« Last Edit: January 06, 2017, 12:03:28 pm by BrianHG »
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #16 on: January 06, 2017, 12:04:24 pm »
No eeprom, but a thought occurred to me.

I have already implemented a function that rotates a point (at angle) about an origin.
If I define a point 16384 pixels below an origin, I should only have to step through a loop of 90 steps for each angle integer,
and read only the y coordinate of the point into the sine array.

Other than the loop, that would be free since the point rotation function already exists.

Finding the code for the video above, that’s how I generated the waves, not with an audio function.

« Last Edit: January 06, 2017, 12:06:04 pm by @rt »
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #17 on: January 06, 2017, 12:48:32 pm »
and done yay :)

It’s the function pixel_rotation that already existed. Normally for rotating a point on display, but I don’t need to draw this.
In the added code, the variables that are cleared are already clear, so those can probably go too.


Code: [Select]
// generate sine table
n = 0;
rotation = 0;
while (n < 91) {
ox = 0; oy = 0;
x = 16384; y = 0;
pixelrotation();
sin_table[n] = y;
rotation++;
n++;
}


void pixelrotation() {
theta = rotation; // pre calculation
thetab = 360 - theta;
thetab = d2r * theta;
xx = x; yy = y;
rx = xx - ox; ry = yy - oy;
xnew = cos(thetab)*rx - sin(thetab)*ry;
ynew = sin(thetab)*rx + cos(thetab)*ry;
xnew = xnew + ox; ynew = ynew + oy;
x = xnew; y = ynew;
}


 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #18 on: January 06, 2017, 01:04:41 pm »
Hold on a sec, you wanted a part of a sine table and you were willing to use the 'sin()' command?  In this case, you've over engineered your code.
Or, did you just add additional functionality to your code?

sin_table[n] = sin(pi*(n/90)/2) * 16384;

I hope I got my #s right...
« Last Edit: January 06, 2017, 01:50:05 pm by BrianHG »
 

Offline rs20

  • Super Contributor
  • ***
  • Posts: 2318
  • Country: au
Re: Code used for sine table
« Reply #19 on: January 06, 2017, 01:08:03 pm »
I guess we all forgot to answer the OPs question by suggesting to use sin() :-)

sin_table[n] = sin(pi*(n/91)/4) * 16384;

I suspect you need to cast to a float in there somewhere, or else 90/91 will evaluate to integer 0 (and also please, use "i" as your iteration variable, "n" almost always means the total count...)
« Last Edit: January 06, 2017, 01:09:56 pm by rs20 »
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #20 on: January 06, 2017, 01:36:41 pm »
In the math line using the Pic XC32/16 compiler, thanks to the sin() function, it should evaluate everything as a float until it is converted into an integer for the table word.  This code is so simple, it is easy to test in M-Chip's compiler and view your table's memory contents.

Here is my test code in basic:
Code: [Select]
const pi = atn(1) * 4
dim as integer n,sin_t

 for n=0 to 90
  sin_t = sin(pi*(n/90)/2) * 16384
  Print sin_t;",";
 next n

END

Here is the output:
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,

There is a slight rounding error of 1 compared to the original table, but, This is the targeted 1/4 sine wave.
« Last Edit: January 06, 2017, 02:03:05 pm by BrianHG »
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #21 on: January 06, 2017, 01:57:00 pm »
Now, for my code, if you want to quickly use the 1/4 sine wave with 4 quadrants, if it were me, I would change the 90s to 127 making the table a little wider and when reading a full sine, I would expect a range from 0 to 511, using bits 8 and 9 in the pointer to quickly flip and mirror the read table read position of 0-127 giving you a full sine wave.  It also makes ignoring higher bits in the pointer address if you are synthesizing variable frequency sine waves from this one reference table by adding a fixed number to an integer, or float, then converting it to an integer.

Remember, everything I've recommended up until now has been in effort to minimize compute time when in use and minimal memory.  If you have the computing processing time, I would just use the sin() command all the time and forget about the table all together.
« Last Edit: January 06, 2017, 02:01:16 pm by BrianHG »
 

Offline snarkysparky

  • Frequent Contributor
  • **
  • Posts: 414
  • Country: us
Re: Code used for sine table
« Reply #22 on: January 06, 2017, 02:08:16 pm »

Page 4 of this paper gives some polynomial approximations to the sin function.  You can use them to build your table.

http://krisgarrett.net/papers/l2approx.pdf

A.1 Sine on [??, ?]
+ 3.03963550927013314332e-01x
Pointwise Error Estimate: 9.54929658551371892153e-01

- 9.33876972831853619134e-02x^3 + 8.56983327795249506462e-01x
Pointwise Error Estimate: 2.03312253648039771447e-01

+ 5.64311797634681035370e-03x^5 - 1.55271410633428644799e-01x^3
+ 9.87862135574673806965e-01x
Pointwise Error Estimate: 1.59772935681895954411e-02

- 1.47740880797318521837e-04x^7 + 7.99858143743132551201e-03x^5
- 1.65838452698030873892e-01x^3 + 9.99450193893262089505e-01x
Pointwise Error Estimate: 6.64973349774765786287e-04

+ 2.17326217498596729611e-06x^9 - 1.93162796407356830500e-04x^7
+ 8.31238887417884598346e-03x^5 - 1.66632595072086745320e-01x^3
+ 9.99984594193494365437e-01x

Pointwise Error Estimate: 1.72333528683927437958e-05
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #23 on: January 06, 2017, 02:15:27 pm »
One last thing about positive and negative sine quadrants, you have 16384 as the sines peak. Since you are using words, I assume you might want to change it to a gain of 32767, making a sine from -32767 to +32767.  Don't use 32768, since that number is actually -1 as a signed word.
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #24 on: January 06, 2017, 02:15:52 pm »
Yes I didn't say I wasn't willing to use the trig.
It still has to generate the table so that can be used fast.

The function Pixel_rotation was already in my graphics code.
The top part just has to be a worthwhile trade off for the values being defined.
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #25 on: January 06, 2017, 02:45:54 pm »
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.
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #26 on: January 06, 2017, 03:15:10 pm »
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 :D

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.
 

Offline magetoo

  • Frequent Contributor
  • **
  • Posts: 284
  • Country: se
Re: Code used for sine table
« Reply #27 on: January 06, 2017, 04:47:19 pm »
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.
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #28 on: January 06, 2017, 04:59:47 pm »
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
 

Offline NivagSwerdna

  • Super Contributor
  • ***
  • Posts: 2495
  • Country: gb
Re: Code used for sine table
« Reply #29 on: January 06, 2017, 05:08:43 pm »
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 };
« Last Edit: January 06, 2017, 09:41:07 pm by NivagSwerdna »
 

Offline AndyC_772

  • Super Contributor
  • ***
  • Posts: 4228
  • Country: gb
  • Professional design engineer
    • Cawte Engineering | Reliable Electronics
Re: Code used for sine table
« Reply #30 on: January 06, 2017, 05:14:08 pm »
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  :-\

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #31 on: January 06, 2017, 05:28:54 pm »
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.




 

Offline tszaboo

  • Super Contributor
  • ***
  • Posts: 7374
  • Country: nl
  • Current job: ATEX product design
Re: Code used for sine table
« Reply #32 on: January 06, 2017, 05:36:34 pm »
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.
 

Offline orolo

  • Frequent Contributor
  • **
  • Posts: 352
  • Country: es
Re: Code used for sine table
« Reply #33 on: January 06, 2017, 07:16:03 pm »
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.
« Last Edit: January 06, 2017, 07:31:39 pm by orolo »
 

Offline hans

  • Super Contributor
  • ***
  • Posts: 1638
  • Country: nl
Re: Code used for sine table
« Reply #34 on: January 06, 2017, 07:59:30 pm »
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.
« Last Edit: January 06, 2017, 08:19:51 pm by hans »
 

Offline orolo

  • Frequent Contributor
  • **
  • Posts: 352
  • Country: es
Re: Code used for sine table
« Reply #35 on: January 06, 2017, 08:14:41 pm »
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
 

Offline snarkysparky

  • Frequent Contributor
  • **
  • Posts: 414
  • Country: us
Re: Code used for sine table
« Reply #36 on: January 06, 2017, 09:03:45 pm »
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



 
The following users thanked this post: jancumps

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #37 on: January 06, 2017, 11:38:10 pm »
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:
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #38 on: January 07, 2017, 12:06:29 am »
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.
« Last Edit: January 07, 2017, 12:25:24 am by BrianHG »
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #39 on: January 07, 2017, 01:23:06 am »
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
« Last Edit: January 07, 2017, 01:28:03 am by BrianHG »
 

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #40 on: January 07, 2017, 01:30:13 am »

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.
 

Offline orolo

  • Frequent Contributor
  • **
  • Posts: 352
  • Country: es
Re: Code used for sine table
« Reply #41 on: January 07, 2017, 01:56:43 am »
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.
 
The following users thanked this post: JPortici, BrianHG

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #42 on: January 07, 2017, 02:32:08 am »
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.
« Last Edit: January 07, 2017, 03:12:08 am by BrianHG »
 
The following users thanked this post: orolo

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #43 on: January 07, 2017, 06:05:41 am »
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.
 

Offline JPortici

  • Super Contributor
  • ***
  • Posts: 3461
  • Country: it
Re: Code used for sine table
« Reply #44 on: January 07, 2017, 08:23:34 am »
Using only 16-bit integer arithmetic, you can recreate the table with relative error less than 0.001 with this:

 :clap: :clap: :clap:
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)
 
The following users thanked this post: orolo

Offline mikeselectricstuff

  • Super Contributor
  • ***
  • Posts: 13745
  • Country: gb
    • Mike's Electric Stuff
Re: Code used for sine table
« Reply #45 on: January 07, 2017, 08:36:46 am »
Using only 16-bit integer arithmetic

Quote
    t  = ((long)......
..well not all 16-bit, but neat nonetheless
Youtube channel:Taking wierd stuff apart. Very apart.
Mike's Electric Stuff: High voltage, vintage electronics etc.
Day Job: Mostly LEDs
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #46 on: January 07, 2017, 09:07:43 am »
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.
 

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #47 on: January 07, 2017, 09:36:41 am »
Using only 16-bit integer arithmetic

Quote
    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.)
« Last Edit: January 07, 2017, 10:05:06 am by BrianHG »
 
The following users thanked this post: orolo

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #48 on: January 07, 2017, 09:57:40 am »
Using only 16-bit integer arithmetic, you can recreate the table with relative error less than 0.001 with this:

 :clap: :clap: :clap:
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?

« Last Edit: January 07, 2017, 10:44:49 am by BrianHG »
 
The following users thanked this post: orolo

Offline JPortici

  • Super Contributor
  • ***
  • Posts: 3461
  • Country: it
Re: Code used for sine table
« Reply #49 on: January 07, 2017, 11:34:41 am »
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.)
 

Offline orolo

  • Frequent Contributor
  • **
  • Posts: 352
  • Country: es
Re: Code used for sine table
« Reply #50 on: January 07, 2017, 01:10:10 pm »
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.

Quote
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  :horse: , 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).
« Last Edit: January 07, 2017, 03:18:15 pm by orolo »
 
The following users thanked this post: BrianHG

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #51 on: January 08, 2017, 09:33:55 am »
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.
 
The following users thanked this post: orolo

Online BrianHG

  • Super Contributor
  • ***
  • Posts: 7733
  • Country: ca
Re: Code used for sine table
« Reply #52 on: January 08, 2017, 11:26:58 am »
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;
--------------------------------------------------------------
« Last Edit: January 08, 2017, 11:32:26 am by BrianHG »
 

Offline jesuscf

  • Frequent Contributor
  • **
  • Posts: 499
  • Country: ca
Re: Code used for sine table
« Reply #53 on: January 09, 2017, 12:01:58 am »
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


Homer: Kids, there's three ways to do things; the right way, the wrong way and the Max Power way!
Bart: Isn't that the wrong way?
Homer: Yeah, but faster!
 
The following users thanked this post: sorin, orolo

Offline jesuscf

  • Frequent Contributor
  • **
  • Posts: 499
  • Country: ca
Re: Code used for sine table
« Reply #54 on: January 15, 2017, 04:33:38 pm »
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).



« Last Edit: January 15, 2017, 04:39:10 pm by jesuscf »
Homer: Kids, there's three ways to do things; the right way, the wrong way and the Max Power way!
Bart: Isn't that the wrong way?
Homer: Yeah, but faster!
 
The following users thanked this post: rs20

Offline CatalinaWOW

  • Super Contributor
  • ***
  • Posts: 5231
  • Country: us
Re: Code used for sine table
« Reply #55 on: January 16, 2017, 12:14:34 am »
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/
 
The following users thanked this post: orolo, BrianHG

Offline @rtTopic starter

  • Super Contributor
  • ***
  • Posts: 1059
Re: Code used for sine table
« Reply #56 on: January 16, 2017, 06:39:23 am »
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?
 

Offline jesuscf

  • Frequent Contributor
  • **
  • Posts: 499
  • Country: ca
Re: Code used for sine table
« Reply #57 on: January 16, 2017, 06:07:31 pm »
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

Homer: Kids, there's three ways to do things; the right way, the wrong way and the Max Power way!
Bart: Isn't that the wrong way?
Homer: Yeah, but faster!
 

Offline jesuscf

  • Frequent Contributor
  • **
  • Posts: 499
  • Country: ca
Re: Code used for sine table
« Reply #58 on: January 16, 2017, 07:29:20 pm »
Attached is the Excel spread sheet I used to test some of the ideas presented in this thread.
Homer: Kids, there's three ways to do things; the right way, the wrong way and the Max Power way!
Bart: Isn't that the wrong way?
Homer: Yeah, but faster!
 


Share me

Digg  Facebook  SlashDot  Delicious  Technorati  Twitter  Google  Yahoo
Smf