2011
08.09

Kronos offers extensive metaprogramming capabilities. In fact, the metacompiler is Turing complete while the run time is not!

In this section, we are going to build a very good approximation of a sinusoidal oscillator. In part 12 we will construct a playable FM-synth from it.

Computing sinusoids

Sinusoidal oscillator is a common theme in DSP. Performance is important, as many applications require oodles of them. There are also surprisingly many different tradeoffs that can be made to optimize an oscillator to a certain application.

Part 10 ends up with a representative example of a certain type of optimized sinusoidal oscillator. It is a recursive method based on filters. This oscillator is cheap to stream and expensive to modulate; it is suitable when frequency updates are rare. The sinusoid is very pure, although some amplitude error may accumulate as it runs on.

If we require an oscillator to use as a FM operator, recursive methods are rarely useful. We could use a phase accumulator and sinusoid mapping, as in part 7. However, going to the C run time for a precise trigonometric function is very expensive, motivating us to build a more efficient, less precise version. The technique of choice is a Taylor series of the cosine function, chosen over sine for its nice symmetry properties.

The first thing for the series expansion is a factorial function. Somewhat of an ‘Hello world’ of functional programming, the function is predictably simple:

Factorial(n)
{
    Factorial = #1
    Factorial = When(n > #1 n * Factorial(n - #1))
}

The distinction is that all numbers used here are prefixed with ‘#’. In Kronos, this means that they are compiler directives. Any computations involving directives are resolved at compile time, meaning that the result of any computation on pure directives appears as a  constant in the machine code. A further distinction is that you can actually branch on directives as opposed to data. This is demonstrated by the ‘When’ clause in the factorial example above.

To optimize the oscillator for our periodic ramp, let us expand ‘sin(2 pi x)’ instead of ‘sin(x)’. This gives us a function with a period of 1, just like our periodic ramp oscillators.

Looking at the Taylor expansion, we can note that every other coefficient is zero. That’s why we can skip over every other coefficient in the series.

Here is a function to compute the n:th non-zero Taylor coefficient:

Sine-Coef(n)
{
	Sine-Coef = Crt:pow(#2 * Pi n * #2 - #1) * Crt:pow(#-1 n - #1) / Factorial(n * #2 - #1)
}

It makes use of exponentiation and the ‘Factorial’ function we wrote previously. It’s probably a good idea to test it, so let’s write a small helper function to sample a function at a set of points:

Sample(func min max steps)
{
	Use Algorithm
	inc = (max - min) / (steps - #1)
	xs = Expand(steps 'arg + inc min)
	Sample = Map(func xs)
}
PS > .\k2cli.exe -e "Sample('(arg Approx:Sine-Coef(arg)) #1 #10 #10)" .\approx.k
K2CLI 0.1
(c) 2011 Vesa Norilo

Sample('(arg Approx:Sine-Coef(arg)) #1 #10 #10) => ( (#1 #6.28319) (#2 #-41.3417) (#3 #81.6052) (#4 #-76.7059) (#5 #42.0587) (#6 #-15.0946) (#7 #3.81995) (#8 #-0.718122) (#9 #0.104229) #10 #-0.0120316)

Seems reasonable.

Evaluating polynomials

One of the most efficient methods of computing a polynomial is the Horner method.

Horner-Compute(x result coefs)
{
	Horner-Compute = result * x + coefs
	(c cs) = coefs
	Horner-Compute = Recur(x result * x + c cs)
}

Horner-Polynomial(x coefficients)
{
	Use Algorithm
	Horner-Polynomial = Horner-Compute(x 0 coefficients)
}

‘Horner-Polynomial’ is the front end to this evaluation algorithm, with a recursive computation realized in ‘Horner-Compute’.

Since our coefficient series skips all the zeroes, the polynomial we want for our sinusoid approximation is actually on x squared. Finally, the entire Horner evaluation should be multiplied by x once more to get the appropriate powers.

Sin(x order)
{
	Use Algorithm
	Sin = x * Horner-Polynomial(x * x
			Reduce(Swap Map(Sine-Coef Expand(order 'arg + #1 #1))))
}

This little routine computes ‘order’ coefficients for sinusoid approximation, proceeding to evaluate the resulting taylor series with a value of ‘x’. Let’s evaluate the most accurate period of this approximation, from -0.5 to 0.5, using the ‘Sample’ function we wrote earlier:

Sample('Approx:Sin(arg #5) - Crt:sin(arg * 2 * Pi) -0.5 0.5 #10) => (-0.00692543 -0.000447333 -1.12653e-005 5.96046e-008 -0 -2.98023e-008 -0 1.12653e-005 0.000446856 0.00692495)
Sample('Approx:Sin(arg #7) - Crt:sin(arg * 2 * Pi) -0.5 0.5 #10) => (-2.13067e-005 -5.96046e-007 -0 5.96046e-008 -0 -2.98023e-008 -5.96046e-008 -0 2.98023e-007 2.08298e-005)
Sample('Approx:Sin(arg #9) - Crt:sin(arg * 2 * Pi) -0.5 0.5 #10) => (-8.74228e-008 -5.96046e-008 -0 5.96046e-008 -0 -2.98023e-008 -5.96046e-008 -0 -2.38419e-007 -1.50996e-007)

Here we see the errors at ten points from -0.5 to 0.5. Increasing the approximation order by two seems to bump our worst-case error down by several orders of magnitude. Ok! Let’s hear it, in part 12.

No Comment.

Add Your Comment


9 + nine =