4
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Maybe this - mooooo.ooo/chebyshev-sine-approximation $\endgroup$ – Juha P Jun 30 at 13:15
  • $\begingroup$ Maybe. Maybe it is an overkill. That is the point of my question - how do I know some approximation is good enough? $\endgroup$ – Suma Jun 30 at 13:19
  • $\begingroup$ 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. $\endgroup$ – Cedron Dawg Jul 1 at 11:42
  • $\begingroup$ 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. $\endgroup$ – Suma Jul 1 at 13:27
3
$\begingroup$

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
| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Suma Jun 30 at 14:07
  • $\begingroup$ Nice one. The maximum error is going to be at the 0.5 boundaries. You should be able to solve that analytically. $\endgroup$ – Cedron Dawg Jun 30 at 14:20
  • $\begingroup$ 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). $\endgroup$ – Suma Jun 30 at 14:25
  • $\begingroup$ @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. $\endgroup$ – Cedron Dawg Jun 30 at 14:29
  • $\begingroup$ 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. $\endgroup$ – Suma Jun 30 at 14:33
1
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Suma Jun 30 at 12:57
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$
0
$\begingroup$

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)
| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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... . $\endgroup$ – Juha P Jun 30 at 16:01
  • $\begingroup$ Do you have some link for more information about this particular computation? What is the assumed x domain for the error? -pi .. pi? $\endgroup$ – Suma Jun 30 at 16:10
  • $\begingroup$ @JuhaP Yep, if memory is expensive and processing is cheap.... $\endgroup$ – Cedron Dawg Jun 30 at 16:47
  • $\begingroup$ @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. $\endgroup$ – Cedron Dawg Jun 30 at 16:55
  • $\begingroup$ I did calculations using Sollya. Here's the calculation: i.postimg.cc/dtnTss4x/sinremez.png and here's the error plot: i.postimg.cc/ncLh1vG3/sinremezerror.png $\endgroup$ – Juha P Jun 30 at 17:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.