Chebyshev approximation

Work is coming along well on my Webassembly based synthesizer. I've drawn upon my earlier notes to create an s-expression→wasm compiler that lives in the browser. I'm on course to start making actual sound with it now that I've overcome one of Webassembly's biggest pain points: its lack of a native sin function.

Sinusoids are the foundation of all sound. No matter how harmonically complex a signal may be, it can always be decomposed into a series of sines. Any audio program without a sin function is obviously crippled. It is unfortunately very tricky to implement a general use sine generating procedure; performance has to be weighed against accuracy, and what may satisfy a composer might be wholly unacceptable for an engineer. I understand Webassembly's decision to divorce themselves from the matter and not include a primitive f32.sin instruction—however convenient that might have been.

This meant that I had to get my sines from somewhere else. My options were:

1. Import Javascript's sin function and have the wasm call that.
2. Draw a sine cycle in wasm's memory, then perform lookups on it with linear interpolation.
3. Use an approximate procedure.

Option 1 was out of the question for performance reasons. Option 2 was probably the sanest course of action, but I've done it before and find it boring. So I went with the third choice and learned something cool in the process. Whether the final product actually sounds good is another question.

I've seen Chebyshev polynomials here and there in signal processing discussions, but I've never had reason to study them up until now. I imagine the Wikipedia article is full of theoretical insight, but I can't say I understand most of it. The only thing I was concerned with were the Chebyshev polynomials of the first kind.

This series of polynomial functions, indexed T₀(u) ... T_n(u), is derived from a recursive procedure, much like with Fibonacci numbers. It's hard to imagine just looking at them, but they allow one to create a new polynomial function over an arbitrary interval that approximates any other function.

They look like this in Scheme.

(import srfi-1)

(define (t0 u) 1)
(define (t1 u) u)
(define (t2 u) (- (* 2 (expt u 2)) 1))
(define (t3 u) (- (* 4 (expt u 3)) (* 3 u)))
(define (t4 u) (+ (- (* 8 (expt u 4)) (* 8 (expt u 2))) 1))
(define (t5 u) (+ (- (* 16 (expt u 5)) (* 20 (expt u 3))) (* 5 u)))
(define (t6 u) (- (+ (- (* 32 (expt u 6)) (* 48 (expt u 4))) (* 18 (expt u 2))) 1))
(define ts (list t0 t1 t2 t3 t4 t5 t6))

Not pretty.

You don't even need to expand the series that far to get an accurate approximation. My Chebyshev sin was derived from the first 7 polynomials and looks like this.

How did I get a sine from that? How did I get that? A theoretically-sparse account of my implementation follows. I suggest you stop reading now if you genuinely appreciate mathematics.

Quick overview

Approximating a function f(x), x ∈ [a,b] with n Chebyshev polynomials involves:

1. Calculating n Chebyshev nodes.
2. Finding the value of f(x) at each of these nodes, using a scaling procedure.
3. Finding a coefficient for each node.
4. Expanding the nodes and their coefficients into a regular polynomial function g(x), using a reverse scaling procedure.

If that all makes sense to you, great. Read on if it doesn't.

Concerning intervals

Each T_n(u) is a function with a minimum/maximum value of ± 1 over the interval [-1,1]. The special properties of T_n(u), u ∈ [-1,1] can be used to approximate a function over an arbitrary interval f(x), x ∈ [a,b]. Two bridging operations facilitate this: one that scales the value of u in [-1,1] to that of x in [a,b], and one in the opposite direction.

In the case of my desired interval [0,2π], I needed the following scales.

Equivalent Scheme code follows.

(define pi 3.14159265359)

(define (fu x) (/ (- x pi) pi))
(define (fx u) (+ pi (* pi u)))    

Both operations will be used once.

Calculating nodes

Chebyshev nodes are points in [-1,1] where the roots of the polynomials are found. The accuracy of the approximate function improves with the number of nodes, but they generally provide diminishing returns. I found seven nodes to be accurate enough for my sine. If there are n nodes, a node x_i in i = 1 ... n is found with the following formula.

Scheme code.

(define (nodes n)
  (reverse (map (lambda (i) (cos (* pi (/ (- (* 2 i) 1) (* 2 n))))) (iota n 1))))

reverse is used here because polynomials are counted down from their greatest term to their lowest.

The first seven.

(nodes 7)

(-0.974927912181866 -0.781831482468131 -0.433883739117678 -1.03411553555107e-13 0.433883739117492 0.781831482468002 0.97492791218182)

f(x) at each node

In their default state, the nodes can be thought of as vector of u values across [-1,1].

(define us (nodes 7))

They need to be scaled appropriately to approximate sin(x) across [0, 2π].

(map fx us)

(6.20441902028024 5.59778869525685 4.50467862091367 3.14159265358967 1.77850668626574 0.68539661192274 0.078766286899612)

(map sin (map fx us))

(-0.0786848661403183 -0.632980105575548 -0.978505648902669 1.18250194500031e-13 0.978505648902704 0.632980105575554 0.0786848661405859)

This vector of f(x) values across the scaled nodes will be called ys.

(define ys (map sin (map fx us)))

Finding coefficients

Since there are seven nodes, there will be seven Chebyshev coefficients in the approximate function. The values in both us and ys are used to find those coefficients. Any coefficient C_k in the approximate function is found with the following formula.

In a few more words, C_k is found by mapping T_k across us, and then getting the product of this vector against ys. The resultant vector is then averaged. Any coefficient past C₀ is multiplied by two. We're then left with a series of coefficients that will weight the terms in our approximation function.

The process, broken down.

These Scheme procedures produce the list of coefficients cs.

(define (avg x) (/ (foldl + 0 x) (length x)))

(define (cn tf u y) (avg (map * (map tf u) y)))

(define cs
  (let ((fs (map (lambda (x) (c x us ys)) (list t0 t1 t2 t3 t4 t5))))
    (cons (car fs) (map (lambda (x) (* x 2)) (cdr fs)))))

cs

(6.11316552934227e-14 
 -0.569230592157212 
 1.96656183401192e-13 
 0.666910822168275 
 -6.86296257918733e-14 
 -0.104032361849377)

Expanding to an approximate function

Since many of these coefficients are close to zero, they can be ignored. Retaining only the significant ones, we're left with g(x), an approximation of sin(x).

Expanding each T_k and replacing each u for the scaling function to x yields the following.

Simplifying algebraically yields the final procedure.

Here is a graph of a similar, older, less-accurate approximation of mine. You can see how it doesn't differ too much from a true sine. The equation above is even more accurate.

As wasm

Using Horner's method to evaluate the polynomial, I came up with the following Webassembly procedure.

(module
  (func $sin (param $x f32) (result f32)
    get_local $x
    f32.const 3.141592653589793
    f32.sub
    f32.const 0.318309886183791
    f32.mul
    tee_local $x
    get_local $x
    get_local $x
    get_local $x
    get_local $x
    f32.const -1.66451778959003
    f32.mul
    f32.mul
    f32.const 4.74829052566064
    f32.add
    f32.mul
    f32.mul
    f32.const -3.09012486790893
    f32.add
    f32.mul
  )
  (export "sin" (func $sin))
)

Not a table in sight.

In the end

The equation above still doesn't sound great. It will buzz at certain frequencies, which suggests an unacceptable level of error. I would suggest using Colin Wallace's approximation for any serious sound work. But at least we learned something.

Acknowledgements

This entry is basically my specific notes on Jason Sachs' excellent article on Chebyshev approximation.