#include <stdio.h>
#define TBLSZ 64
#define PI 3.14159265358979323846
#define SCALE 32767
int my_sin(int i){
double rad = 2*PI*i/TBLSZ, term=rad, sin=rad, old_sin=0;
for (int n=3; sin != old_sin; n+=2){
old_sin = sin;
term *= (rad*rad) / (-n * (n-1));
sin += term;
}
return (sin * SCALE) + (sin >= 0 ? 0.5 : -0.5);
}
int main(){
printf("short sintbl[] = {\n");
for (int i=0; i<TBLSZ; ++i)
printf(" %6d,\n", my_sin(i));
printf("};\n");
return 1;
}
Sorry, guys, of course you are right. I should just generate the LUTs for 64/128/256/512/1024 values per cycle and #include them according to the SPC (samples per cycle) value which is #defined in the embedded code.
I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC. Last time I did that sort of thing was in the days of Borland C++ Or I would have used Basic, which came with every "IBM PC". I have found this: https://www.freebasic.net/.
Re "#error Workflow error" that is something else I wondered about: how does one generate warnings, or fatal compilation errors, in GCC. That's very useful.
What, you want to custom implement sin() only to run at the init, generating a fixed LUT, only because the size of that LUT is a compile-time constant?
I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC. Last time I did that sort of thing was in the days of Borland C++ Or I would have used Basic, which came with every "IBM PC". I have found this: https://www.freebasic.net/.
I should just generate the LUTs for 64/128/256/512/1024 values per cycle
I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC. Last time I did that sort of thing was in the days of Borland C++ Or I would have used Basic, which came with every "IBM PC".
I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC.[/url]
One problem is that the table size (64/128/256/etc) is set at compile time, otherwise I could just precompute the table, and say a 128 sized table would be just 256 bytes and that would be it. Yes of course anything is possible; I could get somebody on freelancer.com, most likely in the former USSR, to knock up a prog for £50 which reads the .h file where the table size is defined and generates it, and then it can be #included, but then I have to make sure the chain of dependencies is correctly set up (I am using ST Cube IDE) etc. A float sin(x) would be a lot neater.
import math
def calculateTable(length):
table = []
for i in range(length):
value = math.cos( 2 * math.pi * i / length );
table.append(value)
return table
n = 256
table = calculateTable(n)
f = open("table"+str(n)+".c", "w", encoding="ansi")
f.write("float table"+str(n)+"["+str(n)+"] = {\n")
for x in table:
f.write(" "+str(x)+"f,\n")
f.write("};\n")
f.close()
I do that all the time. This is why python is amazing.
dim twopi as double
dim x as double
twopi = 2*3.1415926
for x = 0 to twopi step twopi/1024
print x, sin(x)
next x
/* phase = x / (2*PI) */
int sinscale(float phase, float scale)
{
/* Optional: Handle negative phases */
if (phase < 0.0f) {
phase = -phase;
scale = -scale;
}
/* Optional: Handle phases outside a full period */
if (phase >= 1.0f) {
phase -= (float)(int)phase;
}
/* Second half of each period is the first half negated */
if (phase >= 0.5f) {
phase -= 0.5f;
scale = -scale;
}
/* Second quarter is first quarter mirrored: */
if (phase > 0.25f) {
phase = 0.50f - phase;
}
/* Convert phase to radians. */
float x = phase * 6.28318530717958647692528676655900576839433879875021164f;
const float xx = x*x;
x *= scale;
/* Power series expansion. First term is x. */
const float term1 = x;
/* term -x^3/3! */
x *= xx;
x /= 2*3;
const float term3 = x;
/* +x^5/5! */
x *= xx;
x /= 4*5; /* 4*5 */
const float term5 = x;
/* -x^7/7! */
x *= xx;
x /= 6*7;
const float term7 = x;
/* +x^9/9! */
x *= xx;
x /= 8*9;
/* Sum terms in increasing order of magnitudes, to minimize rounding errors. */
x -= term7;
x += term5;
x -= term3;
x += term1;
/* Rounding. */
return (x < 0.0f) ? (x - 0.5f) : (x + 0.5f);
}
If you are interested in why one would implement the power series expansion of sin() that way: we want to do the addition of the terms in order of increasing magnitude, to minimize domain cancellation errors.#define HALFPI 1.57079632679489661923132169163975144209858469968755291
#define PI 3.14159265358979323846264338327950288419716939937510582
#define TWOPI 6.28318530717958647692528676655900576839433879875021164
template <class T, int N, T M>
class sine_table {
public:
T value[N];
constexpr sine_table(): value() {
for (int i = 0; i < N; i++) {
double x = TWOPI * i / N;
double s = (double)M;
if (x >= PI) {
x -= PI;
s = -s;
}
if (x >= HALFPI) {
x = PI - x;
}
const double xx = x*x;
/* Power series begins with the x term */
x *= s;
const double t1 = x;
/* -x**3/3! */
x *= xx;
x /= 2*3;
const double t3 = x;
/* +x**5/5! */
x *= xx;
x /= 4*5;
const double t5 = x;
/* -x**7/7! */
x *= xx;
x /= 6*7;
const double t7 = x;
/* +x**9/9! */
x *= xx;
x /= 8*9;
const double t9 = x;
/* -x**11/11! */
x *= xx;
x /= 10*11;
const double t11 = x;
/* +x**13/13! */
x *= xx;
x /= 12*13;
const double t13 = x;
/* -x**15/15! */
x *= xx;
x /= 14*15;
const double t15 = x;
/* +x**17/17! */
x *= xx;
x /= 16*17;
const double t17 = x;
/* -x**19/19! */
x *= xx;
x /= 18*19;
const double t19 = x;
/* +x**21/21! */
x *= xx;
x /= 20*21;
x -= t19;
x += t17;
x -= t15;
x += t13;
x -= t11;
x += t9;
x -= t7;
x += t5;
x -= t3;
x += t1;
value[i] = (x < 0.0) ? (x - 0.5) : (x + 0.5);
}
}
};
#undef HALFPI
#undef PI
#undef TWOPI
gives about 50-bit approximation at compile time (±8 ULP in double precision). For example,#if N-0 == 0
#error Sine table size (N) is zero
#elif S-0 == 0
#error Sine table scale (S) is zero
#elif N == 4096 && S == 32000
#include "sin-4096-32000.h"
#elif N == 16384 && S == 32000
#include "sin-16384-32000.h"
#elif N == 65536 && S == 32000
#include "sin-65536-32000.h"
#else
#error Unknown sine table size (N) and/or scale (S).
#endif
with whatever macro names (N and S) you prefer. In the preprocessor tests, N-0 yields zero for both N=0 and when N is undefined, and generates an error if N is not an integral constant.
If you are interested in why one would implement the power series expansion of sin() that way: we want to do the addition of the terms in order of increasing magnitude, to minimize domain cancellation errors.
For N=65536 and S=32000, i.e. sinscale((float)i/(float)N, S) ≃ (int)round(S*sin(2*PI*i/N)), 99% of the values are exact, and the rest only differ by ±1. Same for N=65536 and S=32767.
#define kPI 3.14159265358979323846
static void MakeSinCosTable(double *pSinTable, double *pCosTable, unsigned int n)
{
if (n == 0)
return;
double x = kPI / (double) n, y = sin(x), Alpha = 2.0 * y * y, Beta = sin(2.0 * x),
Si = 0.0, Ci = 1.0;
for (unsigned int i = 0; i < n; i++)
{
if (pSinTable != NULL)
pSinTable[i] = Si;
if (pCosTable != NULL)
pCosTable[i] = Ci;
double Si1 = Si - (Alpha*Si - Beta*Ci);
double Ci1 = Ci - (Alpha*Ci + Beta*Si);
Si = Si1;
Ci = Ci1;
}
}