Here is what I'd use, assuming single-precision floating-point hardware is available, noting that
sinscale(p, s) ≃ (int)round(s * sin(2 * PI * p))/* 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.
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.
For
N=1048576 and
S=262144, 91% of the values are exact, and the rest only differ by ±1.
In general, for
S ≤ 262144, I do believe the maximum difference is ±1.
This works well for non-power of two
N, including primes.
If you are using C++, then consider expanding what
emece67 suggested, into something like the following:
#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,
constexpr sine_table<int16_t, 1048573, 32767> table1; constexpr sine_table<int16_t, 65536, 32000> table2;take annoyingly long to compile, but the generated tables have
table1.value[i] == round(32767*sin(2*TWOPI*i/1048573)); table2.value[j] == round(32000*sin(2*TWOPI*j/65536));with exact values including correct rounding (towards nearest integer, halfway away from zero).
If you use precalculated tables, I recommend you define preprocessor macros for the size (
N) and scale (
S), and use e.g.
#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.