[music-dsp] fast sin(x) function

[Ross Bencina]

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, 
October 1992.
http://ccrma.stanford.edu/~jos/wgo/

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...

Cheers

Ross.




----- 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 
> other
> 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 
frequency response.

________________________________________________________________

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 
more efficient:

   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))

or

    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[0]).  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[0] and x[1]) 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 
code.

So we chuck the input to the filter, x[n], and get this recursive difference 
equation

    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)
and
    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 
links
http://music.columbia.edu/cmc/music-dsp
http://music.columbia.edu/mailman/listinfo/music-dsp