🚀 go-pugleaf

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?
>
>Thanks,
>Mike

I implemented a sine-wave generator using a DSP56001 using the following
algorithm:

x(t) = x(t - Ts) * exp (j * w * Ts)

where x(t) is a complex value, Ts ist the sampling interval, w ist the
desired frequency in radians/s.
The algorithm has some constraints regarding the amplitude accuracy:
using 24 bit arithmetic, the amplitude changes in a way that it is just
visible on an analogue level meter. Doing the calculation with double
precision (48 bit) should obtain a stable amplitude for a very long time.
using 32 bit on another DSP should be adequate to guarantee an implitude
accuracy of fractions of a percent over several minutes.

Norbert Hanke
Swiss Federal Institute of Technology
ETH Zurich, Switzerland

I believe most people do this with ROM look-up.  Some DSP chips
have this table as part of the default ROM, but for high speeds,
give up on DSP chips and use a commercial phase-accumulator chip.
Sorry, I don't know any part numbers.

David L. Rick
drick@hplvdo.HP.COM

Disclaimer:  My employer would rather you bought an HP product.

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.

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)

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?


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

A standard CG trick to generate circle points is the "register saving" technique
originating from a programmer's mistake.

Rotating a point around the origin can be done by multiplying by a
rotation matrix:

                /  cos(phi)  -sin(phi)  \
[x' y'] = [x y] |                       |
		\  sin(phi)   cos(phi)  /

As an approximation for small angles you can use the matrix

    / 1       -phi \
    |              |
    \ phi       1  /

A program rotating the point:

   tmp = x + y * phi
   y = y - x * phi
   x = tmp

The mistake was:

   x = x + y * phi
   y = y - x * phi  (saving one scratch register)

This corresponds to the matrix

    / 1       -phi \
    |              |
    \ phi  1-phi^2 /

which has a determinant equal to 1. This algorithm draws an ellipse very close
to a circle. For further information see a CG book like
Graphics Gems, by A. Glassner.

I don't know whether this algorithm is stable, but you could restart it
after one full circle and use the X component as output.

Ge' Weijers
--
Ge' Weijers                                    Internet/UUCP: ge@sci.kun.nl
LCN Foundation
Reguliersstraat                                tel. +3180238130 (UTC+1,
6525 ED Nijmegen, the Netherlands               UTC+2 march/september

How about this for a sine wave generator ::

x = initial val
y = initial val
epsilon = small number

loop repeat 
 { x =  x + epsilon*y
   y =  y - epsilon*x
 }