[music-dsp] fast sin(x) function
rossb-lists at audiomulch.com
Sat Jun 21 04:38:37 EDT 2008
Looks good Robert...
To "the Editor": might be worth adding a "see also" for:
``The Second-Order Digital Waveguide Oscillator'', by Julius O. Smith III
and Perry R. Cook, Proceedings of the International Computer Music
Conference (ICMC-92, San Jose), pp. 150-153, Computer Music Association,
Which compares three methods of recursive sinusoid generation (and I quote):
1. the coupled form which is identical to a two-dimensional vector rotation,
2. the modified coupled form, or ``magic circle'' algorithm, which is
similar to (1) but has ideal numerical behavior, and
3. the direct-form, second-order, digital resonator with its poles set to
the unit circle.
I'm not whether sure (2) above corresponds to Robert's discussion below.
As for lookup tables, I'd love to see some discussion of efficient
high-resolution modulo phase increment computations. Often 32 bits of phase
doesn't provide enough frequency precision for audio applications.. I've
been using a double precision float, limiting the range to the LUT length
with if() and cast to int using FISTP (on IA32) -- seems like there may be
some other (faster?) techniques though, perhaps 64 bit integer phase
variables are a practical alternative these days...
----- Original Message -----
From: "robert bristow-johnson" <rbj at audioimagination.com>
To: "A discussion list for music-related DSP" <music-dsp at music.columbia.edu>
Sent: Saturday, June 21, 2008 4:33 PM
Subject: Re: [music-dsp] fast sin(x) function
> ----- Original Message -----
> From: "douglas repetto" <douglas at music.columbia.edu>
> To: "A discussion list for music-related DSP"
> <music-dsp at music.columbia.edu>
> Subject: Re: [music-dsp] fast sin(x) function
> Date: Sun, 15 Jun 2008 18:23:30 -0400
> bastian.schnuerle wrote:
> > how ?
> > Am 15.06.2008 um 06:40 schrieb douglas repetto:
> >> I don't suppose anyone wants to take a crack at summarizing this
> >> thread?
> >> It's a topic that has come up many times over the years and would be
> >> good
> >> to have in the FAQ. There's a short entry in the FAQ now (#4 What's a
> >> fast, stable, accurate method for sinewave generation?), but this
> >> thread
> >> has had much more useful info.
> >> Anyone?
> If you're interested in contributing to the FAQ or book reviews or any
> part of the site, simply send me some text in an email and I'll add it. If
> you have a little extra time you can send it as html in the same format as
> the existing entries.
> Andy Farnell and r b-j are tackling the sine generation FAQ!
my apologies for the delay in getting back to this. i've had a sorta
life-changing event happen in the meantime and have been incommunicado. if
this stuff gets into a FAQ, i would recommend or ask that someone copyedit
it (merge with what Andy said, make it clearer or more consistent, etc). i
hope and think it's accurate enough.
Four sorta different techniques for computing sin(x) will be summarized
here. As far as I can tell, there are about two immediate purposes or
applications. One is to generate sinusoidal signals; where the argument of
the sin(x) function continues to increase without bound and the algorithm
doesn't break. The other purpose for evaluating sin(x) (or cos(x)) is
simply for that: to evaluate a trancendental function (like exp(x) or log(x)
or similar), usually the argument is of limited domain and since sin(x) is
periodic, that limited domain can be one period. A music-dsp application of
the latter is in coefficient generation for filters or in evaluation of
1. Look Up Table (LUT), often with interpolation.
The LUT would contain essentially one period of the sin(x) function. This
table could either be created as constant data computed at compile time, or
sometimes could be created by an initialization routine:
for (n=0; n<N; n++)
sin_table[n] = sin( (TWOPI/(double)N) * (double)n );
To make the sin function periodic, modulo N arithmetic is used to look up.
Here, x, is the angle in radians.
sin_val = sin_table[ (int)floor(((double)N/TWOPI)*x + 0.5) % N ];
N could be any integer, but if it's not a power of 2, the modulo operator %
is as costly as division, hardly fast. if N is a power of 2, then this is
sin_val = sin_table[ (int)floor(((double)N/TWOPI)*x + 0.5) & (N-1) ];
We are assuming (hoping) that the compiler is smart enough to evaluate
((double)N/TWOPI) and (N-1) at compile time, otherwise this won't be as fast
as it could be. also, if x is known to be non-negative, then the floor()
function is unnecessary.
Above is LUT with no interpolation. If memory is cheap and the table is
very large (maybe N=2^16 or more), interpolation may not be necessary and
without interpolation the speed will be at least twice as fast as with
interpoation. If linear interpolation is used (I don't want to bother with
higher order polynomial interpolation such as cubic Lagrange or Hermite
interpolation), then one extra (redundant) point in the table is needed and
the generation of the table would look like this:
for (n=0; n<=N; n++)
sin_table[n] = sin( (TWOPI/(double)N) * (double)n );
And the table lookup would look like (assuming N is a power of 2):
x = ((double)N/TWOPI)*x;
n = (int)floor(x);
fract_x = x - (double)n;
n = n & (N-1);
sin_val = sin_table[n];
sin_val = sin_val + fract_x*(sin_table[n+1] - sin_val);
I think the concepts behind this are pretty elementary (actually more like
high school) but I thought it should be spelled out formally. I suppose
programming tricks can be used to make this faster, but I wanted to spell
out the main concept.
2. 2nd order harmonic oscillator.
This is effectively a 2-pole linear filter with complex conjugate poles
lying directly on the unit circle. I wouldn't call those poles "stable",
but some textbooks call them "marginally stable" or "critically stable". If
Fs is the sampling frequency and f0 is the target frequency of the sinusoid,
the conjugate poles are at
p1 = exp(+j*2*pi*f0/Fs) = exp(+j*w0)
p2 = exp(-j*2*pi*f0/Fs) = exp(-j*w0)
The transfer function of this filter is
H(z) = Y(z)/X(z) = z^2 / ((z-p1)*(z-p2))
= z^2 / ((z-exp(+j*w0))*(z-exp(-j*w0)))
= z^2 / (z^2 - 2*cos(w0)*z + 1)
= 1 / (1 - 2*cos(w0)*z^(-1) + z^(-2))
Y(z)*(1 - 2*cos(w0)*z^(-1) + z^(-2)) = X(z)
Y(z) = X(z) + 2*cos(w0)*z^(-1)*Y(z) - z^(-2)*Y(z)
which leads to a difference equation (in the time domain) of
y[n] = x[n] + 2*cos(w0)*y[n-1] - y[n-2]
Now a very good question would be "What the hell is the input, x[n]? Ain't
this supposed to be an oscillator? Other than control parameters, what
input should exist for an oscillator?" There are two ways to answer the
question. Let's say that the whole thing starts with n=0 (the first output
sample is y). Either the initial states, y[-1] and y[-2] are zero and
then the input x[n] must be some function (possibly with two non-zero values
at x and x) that is necessary to start this thing up with the target
amplitude and phase *or* x[n] is always zero and the initial states y[-1]
and y[-2] are set to whatever they need to be to start this up with the
desired amplitude and phase (the only other attribute of the sinusoid is the
frequency and that is determined by the coefficient, 2*cos(w0). I'm pretty
sure that the latter way to look at it is both simpler conceptually and to
So we chuck the input to the filter, x[n], and get this recursive difference
y[n] = 2*cos(w0)*y[n-1] - y[n-2]
This defines our present sinusoid value in terms of the previous two. The
initial states, y[-1] and y[-2] are simply what the output samples would
have been if the oscillator was running for those two samples. So we
extrapolate backwards. We know the normalized frequency is w0. Let the
amplitude, A, and phase, theta, be given. Then the output is expected to be
y[n] = A*sin(w0*n + theta)
Then define simply
y[-1] = A*sin(-w0 + theta)
y[-2] = A*sin(-2*w0 + theta)
Use those two initial conditions to start up your oscillator and the above
recursion formula, and you're golden. You have a sinusoidal oscillator.
One thing that must be mentioned is that I must credit James McCartney (Mr.
SuperCollider) with correcting a misconception I had regarding this simple
alg. These two complex-conjugate poles lie precisely on the unit circle
because the distance the poles are from the origin is the square root of
whatever multiplies the z^(-2) term or the y[n-2] term, and there are no
numerical problems implementing multiplication by the number 1. You can do
it exactly without any precision error. This means that the poles are
neither a teeny milli-smidgen outside the unit circle (which would cause the
sinusoid amplitude to grow slowly, but exponentially, until some numerical
limitation, the "rails", stops it) nor a teeny milli-smidgen inside the unit
circle (which would cause the sinusoid to decay slowly and exponentially
until it hits the noise floor). I originally thought that other numerical
roundoff issues (hey, nothing is perfect) would cause the amplitude to grow
or decay and that some kinda amplitude control (just as is done in analog
oscillators) would be needed to keep the amplitude pinned to a target value.
I took James on, here at music-dsp, regarding this issue and the fact is he
was right and I was wrong (I think I admitted that on the list later when I
checked it and studied it a bit). This thing will run forever without the
amplitude growing or decaying. I think that is true for any selected
frequency of oscillation. Roundoff errors in the coefficient 2*cos(w0)
result only in the oscillator frequency being slightly off of the target
frequency. Dunno if such is true for fixed-point and finite precision (if
the roundoff is always toward zero, that would result always in errors that
would reduce amplitude with no hope of errors that would increase it), but
this works fine for even single-precision floating point.
i'm outa time for now, but the next installment will be the quadrature
oscillators (useful for frequency shifting if you also toss in a Hilbert
Transformer), and finally the finite-order polynomial approximation to
sin(x) which we have already seen in this thread.
r b-j rbj at audioimagination.com
"Imagination is more important than knowledge."
dupswapdrop -- the music-dsp mailing list and website:
subscription info, FAQ, source code archive, list archive, book reviews, dsp
More information about the music-dsp