Once again!!
Once more I find that I’d be in deep do-do without Lisp…I had to resurrect an ancient DSP demo for ARM processor this afternoon, and it was clear to my listening that something just wasn’t quite right about the operation.
So break out the Lisp, and recode the ARM C into parallel and equivalent Lisp to allow us to look for errors. Doing this on the ARM development environment is tedious and slow and painful. But in Lisp, it is a breeze because I can introspect so easily, plot graphs of function outputs, etc.
Here is a bit of code from the ARM and the errors that Lisp allowed me to find:
This routine computes a base-2 antilog (2^x) for x in x.18 format, with results in x.18 format. That format is a fixed point fractional value with 18 bits of fraction.
inline F23_t f18_floor(F23_t x)
{
// for x >= 0
return (0xff800000 & x); ;; <— OOPS!! Lisp found this one… that should be 0xFFFC0000 for 18 bits
}
inline F23_t f18_ceiling(F23_t x)
{
// for x >= 0
return f18_floor(x + 0x007fffff); ;; <— and this one too, should be 0x0003FFFF
}
F23_t f18_pow2(F23_t x)
{
// x is x.18 format
// result is x.18 format
// good to +/-0.0025 for x in [0.0,1.0)
// (1.4139997232692892D0 0.49752394360425795D0 0.08600027673071043D0)
int32_t expon;
if(x < 0)
expon = -(f18_ceiling(-x));
else
expon = f18_floor(x);
x -= expon; // x now in [0,1) in x.18 format
expon >>= 18; // expon now integer
// x in [0,1) in 24-bit frac form. Needs to be range reduced to [-1,1)
// x' = 2*x - 1
// shift 1 bit for *2, shift 5 bits for x.18 to x.23
x = (x - F23(0.5)) << (1 + 5); ;; <<<— and this one too! Should be F18(0.5)
x = f23mul(x, f23mul(x, F22(0.08600027673071043))
+ F22(0.49752394360425795))
+ F22(1.4139997232692892);
expon -= 4; // for x.22 to x.18 format
if(expon > 0)
return (x << expon);
else if(expon < 0)
return (x >> (- expon));
else
return x;
}