How to calculate a sine wave with accuracy in time

Question

Using the case is to generate a sine wave for digital synthesis, so we need to calculate all the values ​​of sin (dt), where:

t is an integer representing the sample number. This is a variable. The range is from 0 to 158,760,000 in one hour of CD quality.

d is double, representing an angle delta. It is permanent. And range: greater than 0, less than pi.

The goal is to achieve high precision using traditional int and double types. Performance is not important.

Naive implementation:

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

But the problem is that when t increases, accuracy decreases because large numbers are provided for the sin function.

An improved version is as follows:

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

Here I do not forget to specify numbers in the range from 0 to 2 * pi to the "sin" function.

But now the problem is that d is small, there are many small additions that each time reduce accuracy.

The question is how to increase accuracy.

---

Attachment 1

"accuracy decreases because large numbers provided by the 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 the approach that is used, these are some implementations of the fast Fourier transform. Trigonometric function values ​​are calculated based on previous values ​​and delta.

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

Here we must save and update the cosine value and save the constant (for given deltas) of the factors Cos (d) and Sin (d).

Now about accuracy: the cosine (d) for small d is very close to 1, so there is a risk of loss of accuracy (there are only a few significant digits in the figures, for example 0.99999987). To overcome this problem, we can store constant factors like

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

using other formulas to update the current value
(here sa = Sin(A) for the current value, ca = Cos(A) for the 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 

PS Some FFT implementations periodically (each step K) update the values ​​of sa and ca through a trigger. to avoid the accumulation of errors.

Result. Doubles calculations.

 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 reduced to the exact finding of 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 greater accuracy than sin(256)

So, let's choose an integer to work with, 256. First, find small numbers 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 divide your index as a base number in 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; 

Now you can 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 will have to calculate different values ​​of modB , modC and modD , which are equal to those powers of 256, but modulo (2π / d). You can use the high precision library for these two calculations.

Answer 3

Increase 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 a billion samples, multiplier accuracy will be good enough.

Answer 4

The following code stores the input of the sin () function within a small range, while slightly 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-precision, the OP has 2 problems:

• multiplying d by n and preserving greater precision than double . This answers in the first part below.

• Execution of the mod period. A simple solution is to use degrees, and then mod 360 , it is easy enough to make accurate. Making 2*π large angles is difficult because it requires a 2*π value with approximately 27 accuracy points than (double) 2.0 * M_PI

---

Use 2 double to represent d .

Suppose 32-bit int and binary64 double . So double has 53-bit precision.

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

Segment d in 2 double dmsb,dlsb , 25 most significant digits and 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 sequential addition) of dmsb*n will be exact. (This is an important part.) dlsb*n will be an error only in the smallest quantities.

 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 accurate, giving 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); } 

Exit

 0.04630942695385031893 

---

Use degrees

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

If the OP still wants to use radians, for a further understanding of the execution of the π mode see Good for the last bit