Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!ub!uhura.cc.rochester.edu!rochester!pt.cs.cmu.edu!o.gp.cs.cmu.edu!dandb
From: dandb+@cs.cmu.edu (Dean Rubine)
Newsgroups: comp.dsp
Subject: Re: Sine-wave generator algorithms
Summary: feedback
Message-ID: <1991Apr24.050126.8668@cs.cmu.edu>
Date: 24 Apr 91 05:01:26 GMT
References: <1750@tekig7.MAP.TEK.COM> <1991Apr15.081550.26195@bernina.ethz.ch>
Sender: netnews@cs.cmu.edu (USENET News Group Software)
Organization: School of Computer Science, Carnegie Mellon
Lines: 67

In article <1750@tekig7.MAP.TEK.COM> gaulandm@tekig7.MAP.TEK.COM (Mike Gauland) writes:
>I'm trying to implement a sine-wave generator using a DSP.  Anyone have
>an especially quick algorithm they'd like to share?

    People usually do this sort of thing by table lookup, and so should
you if memory's not a problem.  I recall some chapters in "Foundations of
Computer Music", Roads and Strawn, Eds will tell you in excruciating detail
how big the tables have to be to acheive a given S/N ratio if you are
truncating or interpolating.  Oh, BTW, there's never any benefit in
rounding the phase accumulater, contrary to popular belief.  The difference
between truncating the phase and rounding the phase amounts to half a sample
phase shift in the result; therefore, truncation will give the same results
as rounding if you compensate for the phase shift either by initializing the
phase accumulator to half a sample or by phase shifting the table accordingly.

    Anyway, what I really wanted to share was this method of generating
sine waves without a table; great if you lack memory.  Just implement
a marginally stable two-pole digital filter.  By putting the poles on the
unit circle, the filter will naturally oscillate at the frequency of the
poles.

     If we put the poles at exp(j * theta) and exp(-j * theta)  (j being
sqrt(-1) of course) we get this difference equation:

	y[n] = x[n] + 2*cos(theta)*y[n-1] - y[n-2]

    The impulse response of this filter will be a sinusoid. (For sine phase,
I think we want the response of an impulse at n=-1 rather that at n=0,
but I'll ignore that here.)   The unit impulse response produces a sinusoid
whose amplitude is 1/sin(theta), so you should start with an impulse whose
non-zero value is sin(theta).  The simple algorithm is thus:


sinewave(theta)
double theta;
{
	double c2 = 2 * cos(theta);
	double y = sin(theta);
	double y1 = 0;
	double y2;

	for(;;) {
		printf("%g\n", y);
		y2 = y1;
		y1 = y;
		y = c2 * y1 - y2;
	}
}

   Pretty slick: one multiply, one substract and some assignments for
each sample.  The drawback is that you need to calulate the sine and
cosine of theta before you start the loop, but since it's just once
you can use a library routine for that.

   You can probably vary the frequency smoothly by playing with c2 while
the oscillator is running, though I've never tried this.

   I seem to recall from a standard graphics text that this method is
slightly unstable: when used to draw circles it produces spirals.
I'm not sure about this, but I seem to remember there was a simple fix,
which unfortunately involved another multiply.

--
ARPA:       Dean.Rubine@CS.CMU.EDU	
PHONE:	    412-268-2613		[ Free if you call from work ]
US MAIL:    Computer Science Dept / Carnegie Mellon U / Pittsburgh PA 15213
DISCLAIMER: My employer wishes I would stop posting and do some work.