Newsgroups: comp.dsp
Path: utzoo!utgpu!cunews!wilf
From: wilf@sce.carleton.ca (Wilf Leblanc)
Subject: Re: Sine-wave generator algorithms
Message-ID: <wilf.672514636@rigel.sce.carleton.ca>
Sender: news@ccs.carleton.ca (news)
Organization: Carleton University, Ottawa, Canada
References: <1750@tekig7.MAP.TEK.COM> <1991Apr15.081550.26195@bernina.ethz.ch> <1991Apr24.050126.8668@cs.cmu.edu>
Date: 24 Apr 91 17:37:16 GMT

dandb+@cs.cmu.edu (Dean Rubine) writes:

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

>[...]


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

Very unstable.  Unless you chose theta, such that the values of
cos(theta) and sin(theta) can be represented exactly in floating
point, the output will eventually spiral in or spiral out.

A better algorithm is:

complex a,y;

    a = exp(j*theta)
    y = 1 + j*0;
    forever
    {
        y *= a;
        output = imaginary_part(y)
        normalize(y)
    }

Normalize just makes sure y remains on the unit circle.  Can be simple
since y is very close (i.e. an iterative algorithm can be used. The algorithm
mentioned in the previous post is very similar to this one, except
this algorithm has a single pole at exp(j*theta), and not complex
conjugate pole pairs, and this algorithm has normalization to ensure
that y remains on the unit circle.

>--
>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.
--
Wilf LeBlanc, Carleton University, Systems & Comp. Eng. Ottawa, Canada, K1S 5B6
Internet: wilf@sce.carleton.ca   UUCP: ...!uunet!mitel!cunews!sce!wilf
                Oh, cruel fate! Why do you mock me so! (H. Simpson)