gaulandm@tekig7.MAP.TEK.COM (Mike Gauland) (04/12/91)
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 -- ------------------------------------------------------------------------------ Michael Gauland KG7OS || "No man's life, liberty, or property is safe, (503) 627-5067 || while the legislature is in session." gaulandm@tekig7.map.tek.com ||--unattributed; can't find the source.
hanke@nessie.cs.id.ethz.ch (Norbert Hanke) (04/15/91)
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
drick@hplvec.LVLD.HP.COM (David Rick) (04/16/91)
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.
dandb+@cs.cmu.edu (Dean Rubine) (04/24/91)
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.
wilf@sce.carleton.ca (Wilf Leblanc) (04/25/91)
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)
ge@dbf.kun.nl (Ge' Weijers) (04/25/91)
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
hqm@ai.mit.edu (Henry Minsky) (04/27/91)
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
}alk@hpspdra.spd.HP.COM (Al Kovalick) (05/04/91)
Several responded with the sine recursion formula method of generating a sine wave. IMHO this is the FASTEST method for generating a sine wave short of the table look up method. Actually, for say 12 bit precision, a chordal fit to a sine wave will yield a faster algorithm at the expense of amplitude resolution and quantization error. In fact with only 2 straight lines fitting from 0 to 90 degrees the max departure error is only about 2.7%. The straight line fit method still needs a look up table but with only 32 or 64 entries for the first quadrant. The other quadrants may be derived thru quadrant logic based on the first. I've always been interested in the tradeoffs made when deriving sine compute algorithms. The speed VS accuracy tradeoff is the most interesting. BTW, the previously mentioned recursion method is used in almost all FFT/IFFT transforms that are worth their salt. Using C and all variables in DOUBLE, I have generated 1 million points with the limit cycle noise less than -100 dB ref 1.0. Also the recursion method is only practical when you need to compute a succession of sine points with equally spaced time increments. For random look up, the chordal fit with small table is best, albeit the amplitude noise may be high. Finally, the recursion method is based on the trig identity SIN(N*a) = 2*SIN((N-1)a)COS(a) - SIN((N-2)a) ! simple, beautiful. Happy generating !! Al Kovalick Hewlett-Packard Stanford Park Division, Palo Alto, Ca.