[sdiy] OT - Linear interpolation in C

[Eric Brombaugh]

Tom,

You've run into one of the main problems with using C on DSPs - there is no data type to explicitly represent fixed point arithmetic. Although you didn't say how you did it, I imagine you've found the best answer in the first section of you code by using shifts and logical operations to do bit selection. When it comes time to do the multiplication things start to get nasty - you can break it down into simpler operations, but then you lose efficiency because while the CPU has fantastic math capabilities, the compiler doesn't know about them and won't use them.

Some mfgs dance around this by using compiler-specific features called intrinsics - essentially wrapping the special capabilities of the processor with small inline functions that allow the compiler to use them - this is crude and effective, but not standardized and hence completely unportable. I don't believe that MCHP has done this for the C30 compiler though, so as handy as it would be it's not an option in this case. You might be able to mimic this by using inline assembly, but that gets ugly and hard to maintain. As you note, floating point would handle the representation, but again won't use the DSP capabilities and pulls in a lot of inefficient support code.

My approach so far has been to bundle the DSP code that I need into small assembly subroutines that are in separate source files and link those in. If you study the way the MCHP DSP library is implemented you can see a number of good examples of this. If you keep the granularity at a coarse enough level you won't lose much to the overhead of function calls and register saves.

Short answer: You can't do it in C on the dsPIC, so keep to assembly for this stuff.

Eric

On Jan 30, 2011, at 7:10 AM, Tom Wiltshire wrote:

> Sorry to post off-topic (sort of - this is an LFO I'm working on). Perhaps any replies would be better to me privately if you think C programming isn't sufficiently synthy. There *does* seem to have been a lot of it about recently.
> 
> The problem I'm having is in converting the following assembly into C.
> I've got a 32-bit phase accumulator PHASE_HI and PHASE_LO. I've got 256-entry waveform tables, and the linear interpolation comes in to make straight line sections between those points.
> 
> Now, what's below is not running debugged code, just something I've knocked up to demonstrate how I normally go about this kind of thing. The plan is to use the top 8 bits as an index into the waveform table, whilst the bottom 8 bits indicate the 'sub sample' position between two table entries.
> 
> The first section (doing the table lookups) is easy, and I've got that working fine in C. It's the second section I'm struggling with. Specifically, I can't see how the C variable types I've got relate to the 1.15 format accumulator usage that I'm used to. I've got a 1.15 integer type in C (signed int), but it is treated as a full range integer (eg 32767 to -32768) rather than as a fractional number from -1 to nearly 1, which is how I use it with the accumlator.
> I could "upgrade" my variables to "float", but that's unnecessary since I already know it can be done efficiently with 16-bit variables, and that I can avoid division by using shifts.
> 
> So how do I get C to treat a signed int as if it runs from -1 to 1 rather than -32768 to 32767? Or is there a better way?
> 
> Thanks,
> Tom
> 
> 
> LinearInterp:
> ; Work out the linear interpolated value
> ; Value = (THIS * 1-index) + (NEXT * index)
> 
> 	; Deal with waveform look-up
> 	mov PHASE_HI, W2			; Get high word of phase
> 	lsr	W2, #8, W2				; Shift down so I only have the top 8 bits
> 	sl	W2, W2					; Shift up to make an index (2 bytes per entry)
> 	mov	#tblpage(LFO_SineTable), W1		; Set up TBLPAG
> 	mov	W1, TBLPAG
> 	mov	#tbloffset(LFO_SineTable), W1		; Get base position
> 	add	W1, W2, W1						; Add offset
> 	; Get the first lookup result
> 	tblrdl	[W1++], W6					; Get the data(THIS into W6)
> 	; Get the second lookup result
> 	tblrdl	[W1], W7					; Get the data (NEXT into W7)
> 
> 	; Use the 8 lower bits of the high word for the linear interp
> 	mov PHASE_HI, W4		; Index 
> 	and	W4, #255, W4		; Mask off the top word
> 	sl	W4, #7, W4			; Shift up to give range 0-1
> 	com	W4, W5				; 1-Index
> 	bclr	W5, #15				; Clear the sign bit again
> 	mpy	W5*W6, A			; (THIS * 1-index)
> 	mac	W4*W7, A			; (NEXT * index)
> 	sac.r	A, W0			; Get rounded result into W0
> 
> (code is dsPIC ASM)
> 
> 
> 
> 
> _______________________________________________
> Synth-diy mailing list
> Synth-diy at dropmix.xs4all.nl
> http://dropmix.xs4all.nl/mailman/listinfo/synth-diy