Back to main

# Audio Programming Lab 5

Exercises:
• 13. Small Angle Interpolation
• 14. A Recursive Oscillator

## Exercise 13: Small Angle Interpolation

### Background

As well as programming a large computer for audio output, it's often useful to be able to minimise the time and storage spent generating a sine-wave. Some devices contain very low power processors but might still need to implement audio-frequency oscillators; others may need to generate sine-waves at very high frequency sample rates. It is now possible to build FM receivers using software-only IF processing (at 10.7MHz) on a relatively modest dsp chip.

The size of the look-up table, even in its quarter-cycle version, would be regarded by hardware designers as a profligate waste of memory, when a simple piece of trig can cut the storage requirement dramatically.

sin(A+B) = sinAcosB + cosAsinB
sin(x) ≈ x; cos(x) ≈ 1 - x2/2 (for small x).

Here's the clever bit: let's store a smaller number of entries in the look-up table, and use the identity above to interpolate between them. So A will be a "big" angle (cos and sine obtained from the table) and B will be the difference between the value which can be obtained by the table lookup and the angle we need. B is a "small" angle, so the approximations can be used to calculate the final value.

Do the Maths Homework to find out how much of the table can be thrown away without compromising the accuracy of the result.

### Maths homework

Before investing effort in coding up an interpolated solution, we'd better find out exactly how much storage can be saved from the look-up table.

• The interpolation described begins to fail as the angle increases. If we're using 16-bit audio and A=0, find how the largest B before the error exceeds 0.5bit (1/65536) by solving
|sin(A+B) - ξ| = 1/216
where ξ = (sinA)(1-B2/2) + BcosA
• How many samples are there between 0 and this B radians in your original look-up table?
• Is the answer different if A = π/2?
• How does the situation change if we ignore the x2 term and simply take cos(x)≈1 for small x?

### Reference: Newton's Method

You can find the roots of (most) differentiable polynomials iteratively using Newton's method. (Hint: it's easy to code in Python if you get bored using a calculator)

${x}_{n+1}={x}_{n}-\frac{f\left({x}_{n}\right)}{f\text{'}\left({x}_{n}\right)}$

Modify the oscillator class you have already written so that it uses a smaller look-up table and interpolation.

## Exercise 14: A Recursive Oscillator

### Background

As described in Chapter 2 of mcm.pdf the position of poles and zeros in Laplacian space characterises a system response. Poles on the left half of the s-plane represent a damped oscillatory response to an impulse stimulus, and poles on the right indicate unstable responses.

Discrete-time systems are characterised by the positions of poles and zeros in the z-plane. Poles inside the unit circle characterise damped responses, and poles outside it unstable ones.

If we could build a discrete-time system with poles on the unit circle at the appropriate frequency, it would be possible to create an oscillator with a single pair of poles (representing a second-order system).

The z-transform is as follows:

$X\left(z\right)=\sum _{n\ge 0}x\left[n\right]{z}^{-n}$

An oscillator has a conjugate pair of poles on the unit circle. Moving the poles from z=1 to z=-1 increases the frequency of oscillation from 0Hz to the Nyquist frequency. Note that as the angle around the circle is increased further still, the poles appear to swap over and move back towards z=1, which is exactly the behaviour of the alias frequency which would be generated if an attempt was made to represent a frequency higher than half the sample rate.

The transfer function of a system with poles at an and zeros at bn in the z-plane can be written out as follows:

$\begin{array}{lll}H\left(z\right)& =& \frac{Y\left(z\right)}{X\left(z\right)}\\ & =& \frac{A\left(z-{a}_{1}\right)\left(z-{a}_{2}\right)...}{\left(z-{b}_{1}\right)\left(z-{b}_{2}\right)...}\\ & =& \frac{\sum _{n=0}^{{n}_{z}}{u}_{n}{z}^{-n}}{1-\sum _{n=1}^{{n}_{p}}{v}_{n}{z}^{-n}}\end{array}$

Considering the definition of the z-transform, it is clear that division by z is equivalent to a delay of one sample in time. The rearrangement into the quotient of two polynomials in z-1 therefore realises the oscillator directly.

### Maths homework

Perform the following preliminary design for a 440Hz oscillator using a second-order (two-delay) section as follows.

• Draw a diagram of the pole positions necessary to realise the oscillator at the sample rate you are using.
• Write down the transfer response by substituting the pole positions into bn in the transfer function.
• Multiply out the denominator and hence determine the values of its coefficients -vn. (Sanity check: the poles are conjugates so if there is any imaginary part left, you've got it wrong). You will have to divide the top and bottom by a power of z to get terms with negative exponents.
• Recalling that for an arbitrary signal s[n], z-1S(z) is the z-transform of s[n-1], determine the oscillator output given its previous two outputs.

This all sounds quite daunting, but it isn't hard after you've worked through what's happening once. Read Chapter 2 of mcm.pdf for more detail.

• Write the necessary Python to generate the oscillation and use it to replace the oscillator class you wrote for the previous exercises. You will need to convince the oscillator that it has been running forever. In `__init__`, it's OK to "cheat" by calling the sine function a couple of times to set the previous samples (actually, you only need to call it once if you take 0 as the previous output).