Linear interpolation in direct digital synthesis

Question

I'm working on a microcontroller DDS project in C and am having some trouble figuring out how to compute a linear interpolation to smooth the output values. The program as it stands now
uses the top 8 bits of a 24 bit accumulator as an index to an array of 8 bit output values. I need to come up with a function that will take the middle and lower byte of the accumulator and generate a value in between the "previous" and "next" value in the array. This would be straightforward enough on fast hardware, but since I'm using a microcontroller I really need to avoid doing any floating point operations or divisions!

With those restrictions, I'm not certain of a way to go about getting an 8 bit interpolated value from my two 8 bit input numbers and the lower 2 byes of the accumulator, which represents the "distance" between the two input values. Thanks in advance for any advice!

CLARIFICATION

DDS = Direct Digital Synthesis

in DDS a waveform is generated from a lookup table using a phase accumulator. The phase accumulator usually contains an integer component and a fractional component. The integer component is used as an index into the lookup table. In simple DDS implementations the fractional part is ignored, but for higher quality output the fractional component is used to interpolate (usually just linear interpolation) between adjacent lookup table values. For the above question we are looking at how to efficiently perform this linear interpolation between two lookup table values for a given fraction, f, where 0 <= f < 1.

Answer 1

Linear interpolation between two values, a and b is (a+b)/2.

This is easy in simple hardware and involves no division or floating point.

Divide-by-2 == shift right one bit.

Answer 2

If you really need the speed, check out AVR assembler.

Answer 3

If you want better precision, i'd suggest checking the lower bits of the accumulator first. For example, if we want 4 output values instead of 1:

Acc += 0x2000;
uint lower_bits = Acc & 0xffff;
int rl = LUT[ Acc >> 16];
int rh = LUT[(Acc >> 16) + 1];
if (lower_bits < 0x4000)
    return rl;
if (lower_bits < 0x8000)
    return (rl * 3 + rh) >> 2;
if (lower_bits < 0xC000)
    return (rl + rh) >> 1;
return (rl + rh * 3) >> 2;

Answer 4

Assuming you have a table of waveform values (either one quadrant or four quadrants, it doesn't matter) then one possible optimisation is to store the delta values betwen successive table values. I.e. if you have e.g. N = 256 and a waveform table LUT[N] then you also have a table of delta values, LUT_delta[N]. The relationship between the two precomputed tables is LUT_delta[i] = LUT[i+1] - LUT[i]. So instead of looking up two consecutive table values, LUT[i] and LUT[i+1], subtracting these to get the delta, then doing the interpolation, you just look up the first table value, LUT[i], and the delta, LUT_delta[i] and then calculate the interpolated value. This requires the same number of table lookups, but fewer math operations. You should be able to do the interpolation with a single multiply-accumulate instruction if you're using a DSP, otherwise it's a multiply + scale + add on a general purpose CPU. Also if you interleave the LUT and LUT_delta values you may be able to look up LUT[i] and LUT_delta[i] with a single read and then unpack the two values.

Pseudo-code:

extract integer LUT index, i, from accumulator // just need a shift for this
extract fractional part of accumulator, f // mask or subtract to get f
get p = LUT[i] // lookup waveform value
get delta = LUT_delta[i] // lookup delta
calculate p_interp = p + p_delta * f // single multiply-accumulate instruction on most DSPs - need scaling on general purpose CPUs

Answer 5

To make linear interpolation without doing division, you should make sure that your denominator is a power of 2.

value(x) = previous,
value(x+1) = next value(x + dx) = previous + (next - previous) * dx

Your question is, how do I compute dx ? The trick is to have your interpolation index (the 16 low bit of your accumulator) computed so that the maximum value (dx = 1) is a power of two :

value(x + dx) = previous + ((next - previous) * index) / 1024

Here, you have computed your step value, so that the maximum step is 1024 and corresônds to dx=1. Index = 512 is for dx=0.5 etc...