# How to determine precision needed for sin approximation used for sound synthesis?

I am doing some procedural sound synthesis in Java. I want to have a sine wave as one of the possible basic sounds. When experimenting with that, I have found the default Java Math.sin implementation causes significant performance issues to me. I would like to replace it with either a table lookup or a simple polynomial (Taylor or Chebyshev) series.

My question is: how to decide what table size or what polynomial degree to use so that the resulting artifacts are inaudible?

I can (and will) perform some basic experiments, but I do not trust my ears and speakers that much.

The sound produced will be a tone with a frequency in range 100 Hz - 10 kHz. I want the result to be indistinguishable from a real sine wave by a human. My sound output will be 16 bit at 44.1 kHz or 48 kHz, if that is important.

1. how many entries in the lookup table do I need, assuming the table will be used for -pi/2 to pi/2 range?

2. If I want to use a polynomial series (Taylor or Chebyshev), what degree does the polynomial have to be?

• Maybe this - mooooo.ooo/chebyshev-sine-approximation – Juha P Jun 30 '20 at 13:15
• Maybe. Maybe it is an overkill. That is the point of my question - how do I know some approximation is good enough? – Suma Jun 30 '20 at 13:19
• On a tangent. Java is a byte code interpreter, which means it is going to be inherently slower than native binary (assuming no JIT). You may want to consider coding in C or C++ for computation intensive real time applications. A bit less of a safety net. – Cedron Dawg Jul 1 '20 at 11:42
• I can do this later if necessary, but one can be surprised how fast JVM is nowadays (assuming JIT - as this is how it is used). Currently it seems I will be able to perform 32 voice synthesis at 48 kHz using less than 10% of one CPU core, which is not negligible, but seems is acceptable at the moment for me. I did my share of C/C++/DSP asm programming, but I do not hesitate to enjoy the productivity and safety of JVM. – Suma Jul 1 '20 at 13:27

This is very good question.

Fortunately for you, I have a very good answer.

Assuming standard CD quality, your sound level has just over 4 significant digits accuracy, thus this is the level you need to be truly indistiguishable.

Is this level necessary for auditory purposes? Let's assume so.

Your fastest solution, by far, is going to be linear interpolation on a segmented domain. You need two tables. The first hold the sine value of the center of each interval. The second holds the first derivative, aka cosine.

To calculate your sine value, multiply your angle by the spacing you are using (so a unit of one corresponds a segment). The integer portion will give you the index into your table and use the fractional portion (-0.5 to 0.5) for the interpolation.

Now the question becomes "How many entries do I need in my table?"

I'll let you take a crack at that.

I have a better answer (I think, you need to test it).

$$\sin( x + d ) = \sin( x ) \cos( d ) + \cos( x ) \sin( d )$$

Back to two tables. One for the broad range ($$x$$) , and one for the fine range ($$d$$).

If the range of $$d$$ is small enough, then you can use (Taylor or find the Remez instead):

$$\cos(d) \approx 1 - x^2/2 + x^4/24$$

$$\sin(d) \approx x - x^3/6 + x^5/120$$

But that will require more computation.

For uber accuracy, which you don't need, you could do the interpolation thing on the fine table value.

I was curious, so here you go:

import numpy as np

#==========================================================
def main():

N_coarse = 128
N_fine   = 128

#---- The coarse table

sc = np.zeros( N_coarse )   # Sine Coarse
cc = np.zeros( N_coarse )   # Cosine Coarse

theSlice = np.pi * 0.5 / N_coarse
theAngle = 0.0

for n in range( N_coarse ):
sc[n] = np.sin( theAngle )
cc[n] = np.cos( theAngle )
theAngle += theSlice

#---- The fine table

sf = np.zeros( N_fine )   # Sine Fine
cf = np.zeros( N_fine )   # Cosine Fine

theSlice /= N_fine
theAngle = 0.0

for n in range( N_fine ):
sf[n] = np.sin( theAngle )
cf[n] = np.cos( theAngle )
theAngle += theSlice

#---- The test

theFactor = N_coarse * 2.0 / np.pi

for a in range( 157 ):
theAngle = a * 0.01
n = theAngle * theFactor

nc = np.floor( n )
nf = np.floor( ( n - nc ) * N_fine )

sine = sc[nc] * cf[nf] + cc[nc] * sf[nf]

print sine, np.sin( theAngle ), sine - np.sin( theAngle )

#==========================================================
main()


This old fart got it right on the first try!

Here are the first few rows of the output, the rest is comparable:

0.0 0.0 0.0
0.00997070990742 0.00999983333417 -2.91234267487e-05
0.0199404285515 0.0199986666933 -5.82381418187e-05
0.0299081647675 0.0299955002025 -8.73354349791e-05
0.0399687249608 0.0399893341866 -2.06092258625e-05
0.0499294807897 0.0499791692707 -4.96884810173e-05
0.059885272753 0.0599640064794 -7.87337263996e-05
0.0699307504776 0.0699428473375 -1.20968599355e-05
0.0798735744039 0.0799146939692 -4.11195652507e-05
0.089808457497 0.089878549198 -7.00917010057e-05
0.0998298073783 0.0998334166468 -3.60926852212e-06
0.109745746461 0.109778300837 -3.25543758203e-05
0.119650774894 0.119712207289 -6.14323952304e-05
0.129543907942 0.12963414262 -9.02346778173e-05

• Sounds interesting. A small optimization: "You need two tables. " I do not. They are both identical, with only a different offset to add when obtaining cos. "I'll let you take a crack at that." .... Hm.That means finding a biggest difference between sin and the piecewise linear approximation of the sin I will use. That will need some effort, I am not sure from what end to tackle this, for the start I will try a brute force numerical test with the resolution about 10x that of the table. Once (if) I am done, I will post my findings here. – Suma Jun 30 '20 at 14:07
• Nice one. The maximum error is going to be at the 0.5 boundaries. You should be able to solve that analytically. – Cedron Dawg Jun 30 '20 at 14:20
• Perhaps, but being a programmer a numerical solution is much easier for me. It has also an added value that I can use the same solution to quickly evaluate possible different implementations (only determining the necessary resolution for them may be a bit more difficult). – Suma Jun 30 '20 at 14:25
• @Suma Agreed, for your purposes, implementing a solution, then measuring its accuracy numerically is perfectly appropriate. The only parameter you need to vary is the segment count. – Cedron Dawg Jun 30 '20 at 14:29
• Max errors (a table for whole 2π range): 128 entries: 3e-4 256 entries: 7e-5 512 entries: 1e-5. For a comparison, a plain table with no derivatives (no interpolation) with 1024 entries has max error of 2e-3. – Suma Jun 30 '20 at 14:33

The best way to implement a real time oscillator is based on complex phasor rotation. It's much faster and much more accurate than look up tables or polynomial approximations and you can you even adjust the frequency on the fly without artifacts.

For details on how it actually works see for example: https://dsp.stackexchange.com/a/1087/3997 with a code example at https://dsp.stackexchange.com/a/9868/3997

• Thanks, this sounds real nice. There is one limitation, though - if I understood it correctly, this has a requirement of a constant frequency (constant Δt). This is true for most sounds I will use, but sometimes their frequency is variable, therefore I am still interested in a general sin approximation. – Suma Jun 30 '20 at 12:57

Terms like ‘inaudible’ are a bit tricky, because they are taking in to account the human perception of audio. I suppose the only legitimate metric would be to run ample double blind experiments with subjects and see if they can successfully determine the difference in precision.

That being said, THD may be a useful metric for you. THD below a sufficiently low level would be inaudible, although what that level is, is debatable. From some books I’ve read it’s usually $$0.1$$% to $$0.01$$%. This is about $$-60$$, to $$-80$$dB. If you go lower than that, you’d start to approach the dynamic range of commercial DACs for audio, at which point the distortion wouldn’t be audible, because it would get swamped by the noise, distortion in the DAC.

As pointed by Cedron Dawg in his answer, to achieve a CD-like quality one needs to have a precision over 4 significant digits (1 bit of 16 bits signed is a relative value 3e-5).

I have written a program which numerically tests various implementations over the 0 .. 2π range, testing 2560K equidistant samples.

With a plain lookup table one needs 64-256K to achieve the goal:

• 256K entries: 1E-5
• 128K entries: 2E-5
• 64K entries 5E-5

Using Cedron's solution of using a first derivation as well (which is roughly equivalent to performing a linear interpolation between the table values) one needs 512 entries, max. errors are:

• 128 entries: 3e-4
• 256 entries: 7e-5
• 512 entries: 1e-5

When using a Chebyshev polynomial a polynomial of degree 9 should be enough, as error are:

• degree 7 (6 multiplications): 4.8E-4
• degree 9 (7 multiplications): 2.4E-5
• degree 11 (8 multiplications): 5.3E-8

Using a Remez-based approximation from a comment (5 multiplications) has a precision 7E-5 in range -π/2 .. π/2:

-2.77615229858732e-12 + x * (0.999696773141 + x * (2.307291219071e-11 + x * (-0.165673079310532 + x * (-1.79764345888365e-11 + x * 7.514377168088889654902e-3))))


Using another Chebyshev based approximation (5 multiplications) from a comment has the precision 6e-7 in range -π/2 .. π/2:

x2 = x*x
x*(x2*(x2*(0.00830629-0.00018363*x2)-0.16664824)+0.9999966)

• Here's Remez based approximation with error of around ±7e-5 : y = -2.77615229858732e-12 + x * (0.999696773141 + x * (2.307291219071e-11 + x * (-0.165673079310532 + x * (-1.79764345888365e-11 + x * 7.514377168088889654902e-3)))); This is 2x faster than std::sin() (which is auto vectorized in GCC). If you want to go faster then SSE/AVX is the way to go... . – Juha P Jun 30 '20 at 16:01
• @Suma I would guess by the assymetry $\sin(\theta) \ne -\sin(-\theta)$ that the range is from 0 to $\pi/4$. That is the minimal span needed for the full set and the coefficients were derived using a least squares best fit. And you need to keep $x$ small. – Cedron Dawg Jun 30 '20 at 16:55