How to compute sine wave with accuracy over the time

Question

Use case is to generate a sine wave for digital synthesis, so, we need to compute all values of sin(d t) where:

t is an integer number, representing the sample number. This is variable. Range is from 0 to 158,760,000 for one hour sound of CD quality.

d is double, representing the delta of the angle. This is constant. And the range is: greater than 0 , less than pi.

Goal is to achieve high accuracy with traditional int and double data types. Performance is not important.

Naive implementation is:

double next()
{
    t++;
    return sin( ((double) t) * (d) );
}

But, the problem is when t increases, accuracy gets reduced because big numbers provided to "sin" function.

An improved version is the following:

double next()
{
    d_sum += d;
    if (d_sum >= (M_PI*2)) d_sum -= (M_PI*2);

    return sin(d_sum);
}

Here, I make sure to provide numbers in range from 0 to 2*pi to the "sin" function.

But, now, the problem is when d is small, there are many small additions which decreases the accuracy every time.

The question here is how to improve the accuracy.

---

Appendix 1

"accuracy gets reduced because big numbers provided to "sin" function":

#include <stdio.h>
#include <math.h>

#define TEST      (300000006.7846112)
#define TEST_MOD  (0.0463259891528704262050786960234519968548937998410258872449766)
#define SIN_TEST  (0.0463094209176730795999323058165987662490610492247070175523420)

int main()
{
    double a = sin(TEST);
    double b = sin(TEST_MOD);

    printf("a=%0.20f \n" , a);
    printf("diff=%0.20f \n" , a - SIN_TEST);
    printf("b=%0.20f \n" , b);
    printf("diff=%0.20f \n" , b - SIN_TEST);
    return 0;
}

Output:

a=0.04630944601888796475
diff=0.00000002510121488442
b=0.04630942091767308033
diff=0.00000000000000000000

Answer 1

You can try an approach that is used is some implementations of fast Fourier transformation. Values of trigonometric function are calculated based on previous values and delta.

Sin(A + d) = Sin(A) * Cos(d) + Cos(A) * Sin(d)

Here we have to store and update cosine value too and store constant (for given delta) factors Cos(d) and Sin(d).

Now about precision: cosine(d) for small d is very close to 1, so there is risk of precision loss (there are only few significant digits in numbers like 0.99999987). To overcome this issue, we can store constant factors as

dc = Cos(d) - 1 =  - 2 * Sin(d/2)^2
ds = Sin(d) 

using another formulas to update current value
(here sa = Sin(A) for current value, ca = Cos(A) for current value)

ts = sa //remember last values
tc = ca
sa = sa * dc + ca * ds
ca = ca * dc - ts * ds
sa = sa + ts
ca = ca + tc

P.S. Some FFT implementations periodically (every K steps) renew sa and ca values through trig. functions to avoid error accumulation.

Example result. Calculations in doubles.

d=0.000125
800000000 iterations
finish angle 100000 radians

                             cos               sin
described method       -0.99936080743598  0.03574879796994 
Cos,Sin(100000)         -0.99936080743821  0.03574879797202
windows Calc           -0.9993608074382124518911354141448 
                            0.03574879797201650931647050069581           

Answer 2

sin(x) = sin(x + 2N∙π), so the problem can be boiled down to accurately finding a small number which is equal to a large number x modulo 2π.

For example, –1.61059759 ≅ 256 mod 2π, and you can calculate sin(-1.61059759) with more precision than sin(256)

So let's choose some integer number to work with, 256. First find small numbers which are equal to powers of 256, modulo 2π:

// to be calculated once for a given frequency
// approximate hard-coded numbers for d = 1 below:
double modB = -1.61059759;  // = 256  mod (2π / d)
double modC =  2.37724612;  // = 256² mod (2π / d)
double modD = -0.89396887;  // = 256³ mod (2π / d)

and then split your index as a number in base 256:

// split into a base 256 representation
int a = i         & 0xff;
int b = (i >> 8)  & 0xff;
int c = (i >> 16) & 0xff;
int d = (i >> 24) & 0xff;

You can now find a much smaller number x which is equal to i modulo 2π/d

// use our smaller constants instead of the powers of 256
double x = a + modB * b + modC * c + modD * d;
double the_answer = sin(d * x);

For different values of d you'll have to calculate different values modB, modC and modD, which are equal to those powers of 256, but modulo (2π / d). You could use a high precision library for these couple of calculations.

Answer 3

Scale up the period to 2^64, and do the multiplication using integer arithmetic:

// constants:
double uint64Max = pow(2.0, 64.0);
double sinFactor = 2 * M_PI / (uint64Max);

// scale the period of the waveform up to 2^64
uint64_t multiplier = (uint64_t) floor(0.5 + uint64Max * d / (2.0 * M_PI));

// multiplication with index (implicitly modulo 2^64)
uint64_t x = i * multiplier;

// scale 2^64 down to 2π
double value = sin((double)x * sinFactor);

As long as your period is not billions of samples, the precision of multiplier will be good enough.

Answer 4

The following code keeps the input to the sin() function within a small range, while somewhat reducing the number of small additions or subtractions due to a potentially very tiny phase increment.

double next() {
    t0 += 1.0;
    d_sum = t0 * d;
    if ( d_sum > 2.0 * M_PI ) {
        t0 -= (( 2.0 * M_PI ) / d );
    }
    return (sin(d_sum));
}

Answer 5

For hyper accuracy, OP has 2 problems:

1. multiplying d by n and maintaining more precision than double. That is answered in the first part below.

2. Performing a mod of the period. The simple solution is to use degrees and then mod 360, easy enough to do exactly. To do 2*π of large angles is tricky as it needs a value of 2*π with about 27 more bits of accuracy than (double) 2.0 * M_PI

---

Use 2 doubles to represent d.

Let us assume 32-bit int and binary64 double. So double has 53-bits of accuracy.

0 <= n <= 158,760,000 which is about 2^27.2. Since double can handle 53-bit unsigned integers continuously and exactly, 53-28 --> 25, any double with only 25 significant bits can be multiplied by n and still be exact.

Segment d into 2 doubles dmsb,dlsb, the 25-most significant digits and the 28- least.

int exp;
double dmsb = frexp(d, &exp);  // exact result
dmsb = floor(dmsb * POW2_25);  // exact result
dmsb /= POW2_25;               // exact result
dmsb *= pow(2, exp);           // exact result
double dlsb = d - dmsb;        // exact result

Then each multiplication (or successive addition) of dmsb*n will be exact. (this is the important part.) dlsb*n will only error in its least few bits.

double next()
{
    d_sum_msb += dmsb;  // exact
    d_sum_lsb += dlsb;
    double angle = fmod(d_sum_msb, M_PI*2);  // exact
    angle += fmod(d_sum_lsb, M_PI*2);
    return sin(angle);
}

Note: fmod(x,y) results are expected to be exact give exact x,y.

---

#include <stdio.h>
#include <math.h>

#define AS_n 158760000
double AS_d = 300000006.7846112 / AS_n;
double AS_d_sum_msb = 0.0;
double AS_d_sum_lsb = 0.0;
double AS_dmsb = 0.0;
double AS_dlsb = 0.0;

double next() {
  AS_d_sum_msb += AS_dmsb;  // exact
  AS_d_sum_lsb += AS_dlsb;
  double angle = fmod(AS_d_sum_msb, M_PI * 2);  // exact
  angle += fmod(AS_d_sum_lsb, M_PI * 2);
  return sin(angle);
}

#define POW2_25 (1U << 25)

int main(void) {
  int exp;
  AS_dmsb = frexp(AS_d, &exp);         // exact result
  AS_dmsb = floor(AS_dmsb * POW2_25);  // exact result
  AS_dmsb /= POW2_25;                  // exact result
  AS_dmsb *= pow(2, exp);              // exact result
  AS_dlsb = AS_d - AS_dmsb;            // exact result

  double y;
  for (long i = 0; i < AS_n; i++)
    y = next();
  printf("%.20f\n", y);
}

Output

0.04630942695385031893

---

Use degrees

Recommend using degrees as 360 degrees is the exact period and M_PI*2 radians is an approximation. C cannot represent π exactly.

If OP still wants to use radians, for further insight on performing the mod of π, see Good to the Last Bit