A better sin(x) approximation

Introduction

Last time I showed a simple algorithm for approximating \sin(2 \cdot \pi \cdot x). The third harmonic was -28.6 dBc, which is not good enough for high-quality sound synthesis.

In this post, I will derive a more accurate approximation polynomial.

Method

The first change we will make in the algorithm is to split the domain [0..1) into four pieces, instead of two, and exploit additional symmetry properties of the sine function. This well-known trick allows us to decrease the input range over which our polynomial P(x) must be accurate.

Now, we pick the order of the polynomial to use. Our previous algorithm used a second-order polynomial so we’ll try a third-order one here. The third-order polynomial has four coefficients, so we’ll be able to enforce four constraints:

P(x) = \alpha_3 \cdot x^3 + \alpha_2 \cdot x^2 + \alpha_1 \cdot x + \alpha_0.

Setting constraints

To ensure the curve starts and ends at the correct amplitude, we should constrain the end points, i.e. P(0) and P(\frac{1}{4}). These should be equal to \sin \left( 2 \cdot \pi \cdot 0 \right ) and \sin \left( 2 \cdot \pi \cdot \frac{1}{4} \right ), respectively. This leads to the following constraints:

P(0) = 0 and P(\frac{1}{4}) = 1.

Now, we still have two addition constraints we can set. We could either choose two additional points our polynomial must intersect, or we could set two constraints on the derivative of the polynomial, P'(x):

P'(x) = 3 \cdot \alpha_3 \cdot x^2 + 2 \cdot \alpha_2 \cdot x + \alpha_1.

I chose to constrain the derivatives at the end points because it makes the piece-wise sin approximation continuous in the first derivative, which is desirable from a spectral point-of-view.

The derivative of \sin(2\cdot \pi \cdot x) is 2 \cdot \pi \cdot \cos(2 \cdot \pi \cdot x) and we equate the derivative of the polynomial at x=0 and x=\frac{1}{4} to it:

P'(0) = 2 \cdot \pi and P'(\frac{1}{4}) = 0.

Solving for the coefficients

Given our constraints, we can now solve the system of equations. The system we must solve is:

\alpha_3 \cdot 0^3 + \alpha_2 \cdot 0^2 + \alpha_1 \cdot 0 + \alpha_0 = 0.
\alpha_3 \cdot \frac{1}{4}^3 + \alpha_2 \cdot \frac{1}{4}^2 + \alpha_1 \cdot \frac{1}{4} + \alpha_0 = 1.
3 \cdot \alpha_3 \cdot 0^2 + 2 \cdot \alpha_2 \cdot 0 + \alpha_1 = 2 \cdot \pi.
3 \cdot \alpha_3 \cdot \frac{1}{4}^2 + 2 \cdot \alpha_2 \cdot \frac{1}{4} + \alpha_1 = 0,

The first and third equations give values for \alpha_0 = 0 and \alpha_1 = 2 \cdot \pi. Using these values and simplifying leads to:

\alpha_3 \cdot \frac{1}{64} + \alpha_2 \cdot \frac{1}{16} + \frac{1}{2}\pi = 1.
\alpha_3 \cdot \frac{3}{16} + \alpha_2 \cdot \frac{1}{2} + 2\pi = 0,

which, after some manipulation, gives:

\alpha_3 = 32\pi - 128 \approx -27.469 and \alpha_2 = -16\pi + 48 \approx -2.2655.

How good is it?

The polynomial P(x) follows the \sin(2\cdot \pi \cdot x) curve (for x=[0..\frac{1}{4}]) very accurately:

As can be seen from the error graph below, the error is less than +/- 1.1% (with respect to full-scale):

The following graph shows the spectrum of the approximated sinusoidal wave. It was generated using an 65536-point FFT and 2129 sine wave periods:

The spectrum graph shows that the third harmonic, which is found at -46.9 dBc, is much lower than the simple sine approximation from the previous post.

Show us the code!

/*
  A better sin(2*pi*x) approximation.
  The caller must make sure that the argument lies
  within [0..1).

  Note: this can be optimized further by using Horner's
  rule for evaluating polynomials, see:
  
  http://en.wikipedia.org/wiki/Horner%27s_method
  
  Omitted here for clarity.
*/

float sin_poly(float x)
{
  const float a3 = -27.469f;
  const float a2 = -2.2655f;
  const float a1 = 6.2832f;

  return a3*x*x*x + a2*x*x + a1*x;
}

float better_sin(float x)
{
  if (x > 0.75f)
    return -sin_poly(1.0f - x);
  else if (x  > 0.5f)
    return -sin_poly(x - 0.5f);
  else if (x > 0.25f)
    return sin_poly(0.5f - x);
  else
    return sin_poly(x);
};

Conclusion

This approximation is quite good and it has a spectral purity approximately equivalent to a 7-bit signal.

Advertisements

One thought on “A better sin(x) approximation

  1. If you want the approximation to be continuous in the first derivative, you need to set P′(1/4) = 0, but you do not need to constrain P′(0). Relaxing the unneeded constraint can lower the error to about 0.43%.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s