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

Question

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?

Answer 1

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

Answer 2

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

Answer 3

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.

Answer 4

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)