Digital Harmonic Oscillators

If you need to generate sinusoidal waweform or more complex harmonic waveform, you have several methods at hand. First you may want to reach for application notes of several major DSP manufacturers. These are Motorola Application Notes, Texas instruments application notes, or Analog Devices application notes.

Methods presented in these a reliable and accurate, are those of the lowest performance. Quite often there are way more important task for the DSP to calculate and the construction of a sine wave os only an elementary step to the goal. In this case more powerful methods of a complex oscillator or second order harmonic oscillator can be used. Description of these algorithms is explained in plain term on SineWaves and another on the Don Cross' DSP pages.

These articles also compare the speed of the above mentioned approaches. But there is some reason, why the slower algorithms are for some applications more suitable than than fast ones.

One of basic signal processing component is Digital oscillator. They are widely used as a replacement for memory intensive sine tables, often achieving higher speeds than the table interpolation. Application like Quadrature demodulators, DCT, DFT, FFT, computer graphics can use digital oscillator in lieu of computationally intensive Taylor series approximations, or instead of memory hungry lookup tables.

Digital oscillator can be constructed using several methods, ordered by descending requirements on CPU resources.
 
Algorithm Speed Quality 

at low freq.

 Accuracy Stability Orthogonality 

of sin/cos pair

 Memory 

requirements

 Modulation 

capability
Polynoms Slow Excellent Good Excellent Excellent Medium Excellent
Lookup Medium Excellent Fair Excellent Good Huge Good
Complex  

Oscillator

 Fast Limited Good Poor Good Low Limited
2nd order  

oscillator

 Fastest Poor Good Excellent Good Low Difficult

Polynoms

One of the best known polynominal series are Taylor series. These can be either done as a development around 0, which is:
sin(x) = x - x3/6 + x 5/120 - x7/3620 ....
Where the division is of course converted into multiplicative polynomial coefficients {0, 1, -1/6, 0 1/120, 0, -1/3720}

This is especially suitable for architectures with double accumulator and with possibility of immediate feedback of an accumulator as multiply operand, also with no pipeline latency in multiplication.

Or by development around another point x0 , which usually converges faster.

y = x - x0
sin(x) = k0 + k1*y + k2*y2 + k3*y3 + k4*y4....
Although fewer powers of the argument are necessary, all coefficients are non zero and require about the same number of multiply-accumulate instructions.  For practical implementation, the number of multiplications with calculating the power can be reduced:
sin(x) = K0 + K1*(y + K2*(y + K3*(y +K4*(y + ... ))))
K0 = k0 ; K1 = k1 ; K2 = k2/k1 ; K3 = k3/k2/k1 ...
Continuous function can be approximated with a polynom. This is used in most of the floating point libraries used for high level languages. On most CISC machines, that do not have a dedicated floating point processor this floating point solution to the sine wave approximation is the most computationally intensive alternative. This is less true for fixed point DSPs, which are optimized for sum of products. Taylor series computation can often beat the linear interpolated lookup Tables in both speed and precision. The memory requirements are moderate, since it is only necessary to store the coefficients, in some implementation few intermediate results.

Complexity:  Depending on the necessary accuracy, n-th order polynom requires n2/2 multiplications

Lookup Tables

Lookup tables might be attractive in some applications, where memory is abundant and arithmetic capabilities of the processor are limited. On majority of DSP processors, however, the other three discussed methods are better on almost on all counts, except for some special cases of a sine wave .

Complex Oscillator

Complex oscillator provides a pair of orthogonal sine and cosine waveform. This can be thought of as a point moving in a circle in the complex plane. Multiplication of complex point with amplitude by a number on unit circle ei rotates it by angle i:

 A*e(x+i) = A*ex *ei
The Complex Oscillator is repeated iterative application of this relation:
Substitute:
Zn= Xn+jYn= A*(sin(x)+j*cos(x)) =>  Xn= cos(x); Yn= sin(x)
Zn+1 = Zn *ei          ; Zn = A*e(x+i) = Xn + jYn
Well, the DSP really has to work with real numbers, so we have to expand (2) into:
Xn+1 = A*cos(x+i) = cos(i)*Xn - sin(i)*Xn
Yn+1 = A*sin(x+i) = cos(i)*Xn + sin(i)*Yn
Substitute also convenient constant names for  incremental angle i:
ei= A + jB = cos(i) + j*sin(i)  =>  A = cos(i); B = sin(i)
Then rotating this "base" by the increment can be achieved using formula
Xn+1 = A*Xn - B*Yn
Yn+1 = B*Xn + A*Yn
Complexity: 2 multiplications + 2 multiply-accumulate per one iteration

Stability: In finite arithmetic, the pair  A+jB is off the unit circle except for four four points. This causes  the magnitude of the oscillator slowly but exponentially converge or diverge. If the A+jB falls inside the unit circle, the
oscillator fades. If the A+jB falls outside the unit circle, the oscillator grows unto overflow.

Solution 1: Reset the oscillator points every n-th oscillation
Solution 2: Use pair of second order oscillators

 


Second Order Oscillator

Second  order filter is  modeling difference equation of harmonic movement:

d2Y/dt2 = -omega2*X
With discrete time differential equation changes into difference equation:
D(D(Y)) = (Yn+1-Yn)-(Yn-Yn-1) = Yn+1+2Yn+Yn-1 = -OMEGA2*Yn
Yn+1 = (2-OMEGA2)*Yn - Yn-1
The coefficient  OMEGA is approaching omega  for the discrete Time period approaching zero.  The real meaning of the K =(2-OMEGA2) can be obtained by solving equation, that assumes increment i in angle for each sample:
sin(x+i) = K*sin(x) - sin(x-i)
K*sin(x) = sin(x+i) + sin(x-i) = 2*cos(i)*sin(x)
hence
Yn+1 = K*Yn - Yn-1        ;  K = 2*cos(i)
 is the iterative formula for the second order harmonic oscillator.

Complexity: 1 multiply-accumulate per iteration

Stability: The convergence of this algorithm is determined by the accuracy of the second coefficient with the Yn-1 term, which is always -1. Since this can be modeled in all integer arithmetic accurately, this oscillator always remains stabile.
 
 


In want Home!                        If you have comments or suggestions, email me at petr@ied.com