How to Create Oscillators in Software

How to Create Oscillators in Software

by Matt Donadio


Oscillators can be created in software directly, using the "sine" function, or they can be calculated indirectly using several different iterative methods. We survey those methods here.

Notation:

First, let's assume that:

f = desired oscillator frequency
w = 2 pi f / Fs
phi = initial phase
y = ouput array containing the oscillator signal

Real Oscillator

Now, creating a real oscillator in software is equivalent to sampling a sinusoid, so

for i = 0 to N
y[n] = sin(w * i + phi)
end

is sufficient. This isn't very efficient, though, because we calculate a sine for each output point (though using a sine approximation can speed things up a little.)

Another way is to use the following recurrence relationship for sinusoids:

y[n] = a * y[n-1] - y[n-2], where a = 2 cos(w)

In pseudocode, it would look like:

a = 2 * cos(w)
y[0] = sin(phi)
y[1] = sin(w + phi)
for i = 2 to N
y[n] = a * y[n-1] - y[n-2]
end

This has a problem, though. The recurrence relationship is a 2nd-order all-pole system with:

H(z) = 1 / (1 - a*z^-1 + z^-2)

Plugging:

a = 1
b = -2cos(w)
c = 1

into the standard quadradic equation solver:

x = (-b +/- sqrt(b*b - 4*a*c)) / 2a

we get

x = (2*cos(w) +/- sqrt(4*cos(w)*cos(w) - 4)) / 2
= (2*cos(w) +/- 2*sqrt(cos(w)*cos(w) - 1)) / 2
= cos(w) +/- sqrt(cos(w)*cos(w) - 1)
= cos(w) +/- sqrt(-1 * (1 - cos(w)*cos(w)))
= cos(w) +/- j*sqrt(1 - cos(w)*cos(w))
= cos(w) +/- j*sqrt(sin(w)*sin(w))
= cos(w) +/- j*sin(w)
= e^+/-jw

Therefore

H(z) = 1 / ((1 - e^jw * z^-1) * (1 - e^-jw * z^-1))

This has a pair of conjugate poles on the unit circle, and as a result, is nearly unstable. Despite this, the above approach tends to work well in practice, depending on the application, numerical representation, and precision used. ([Fre94] has an excellent discussion of this method.)

Complex Oscillator

For a complex (quadrature) oscillator, we can use the above method using cosine for the real arm and sine for the quadrature arm (you can actually generate both at the same time; see the Frerking for details). Another method is to view the problem as a rotating vector. We can create an initial vector and rotate is using a complex multiply for each point. You can derive this method from the cos(A + B) and sin(A + B) trig relationships.

dr = cos(w) /* dr,di are used to rotate the vector */
di = sin(w)
y[0].r = cos(phi) /* initial vector at phase phi */
y[0].i = sin(phi)
for i = 1 to N
y[n].r = dr * y[n-1].r - di * y[n-1].i
y[n].i = dr * y[n-1].i + di * y[n-1].r
end

This approach still suffers from stability problems, though. The new recurrence relationship is

y[n] = e^jw * y[n-1]

with

H(z) = 1 / (1 - e^jw * z^-1)

This system has a single pole on the unit circle and is unstable. In this case, we can deal with it. We know that the length of the vector is 1, so we can normalize each point as we create it. If we define the magnitude of y[n] as

|y[n]| = sqrt(y[n].r * y[n].r + y[n].i * y[n].i)

we can add

y[n].r = y[n].r / |y[n]| 
y[n].i = y[n].i / |y[n]|

inside of the "for" loop to normalize the magnitude of y[n]. We loose efficiency, though, because we calculate a square root for each point. However, if x is close to 1, then

   1 / sqrt(x) ~= (3 - x) / 2

Since the length of the rotating vector will always be close to 1 (it will only stray from 1 due to roundoff error), we can use this relationship. We can change the AGC inside the for loop to

mag_sq = y[n].r * y[n].r + y[n].i * y[n].i
y[n].r=y[n].r * (3 - (mag_sq)) / 2
y[n].i=y[n].i * (3 - (mag_sq)) / 2

Another method for creating an oscillator is to use a lookup table. The key to understanding the lookup table approach is understanding the concept of the phase accumulator.

Let's look at a sampled sinusoid

for i=0 to N
y[n] = sin(w * i + phi)
end

Since

sin(x) = sin(x + 2pi) = sin(x - 2pi)

we can write the sampled sinusoid as

for i = 0 to N
y[n] = sin((w * i + phi) mod 2pi)
end

Since the variable "i" is monotnically increasing by 1 each iteration, we can rewrite this as

dw = phi
for i=0 to N
y[n] = sin(dw mod 2pi)
dw=dw + w
end

and again as

dw = phi
for i = 0 to N
y[n] = sin(dw)
dw =(dw + w) mod 2pi
end

Now, here is the key. A very easy way to do modulo addition is to use two's complement math. So, if we make dw a two's complement register, the mod 2pi is implicitly done, and we can also use it to index into a precalculated lookup table. If we have M bits, then the most negative number will represent -pi and the most positive number will represent pi (minus a smidge).

In practical applications, if we need N bits of phase resolution in the lookup table, we only use a table with 2 ^ N/4 entries and use the symetrical properties of sine to choose the proper entry. We also usually use a phase accumulator, delta phase, and initial phase with M bits, and M > N. When we index into the lookup table we take N top bits of the phase accululator (either round or truncate) as the index.

Reference:

     [Fre94] Digital Signal Processing in Communication Systems by Marvin E. Frerking