Was thinking hard about fract math on blackfin yesterday afternoon/evening & this morning. Kind of had an epiphany about how to make sense of the different fractional multiplication primitives with different fract base inputs/output. Promise I will condense these aleph-centric ‘how-to’ posts into a series of blog entries or a pdf at some point. Bit heavy for a forum & not very ‘discussional’, but some aleph users at least professed an interested in aleph DSP. Apologies, re-reading this post it’s not very ‘playful’. In summary I condense a general method for multiplication in any fractional base on the bfin.
I banged my head against this a lot! After muchos random hacking, it still wouldn’t stick to my brain, so I would like to describe the simple(st) way of thinking about it! IMO pretty much essential to get one’s head round this in order to author new DSP blocks. @test2, @zebra I wonder if you guys agree with this analysis!
Anyway to explain the problem & define terminology:
when you call mult_fr1x32x32in bfin codebase on 2 signed 32 bit numbers (aka fract32), the fixed-point integer FR32_MIN (aka 0x80000000) represents -1.00000000 and FR32_MAX (aka 0x7FFFFFFF) represents 0.9999999995343387. The C type is called fract32, but let’s pick a more instructive notation. Since it’s a signed type, the sign already gobbles up one bit of accuracy. Therefore the ‘raw’ fract32 type (as interpreted by the mult primitives) has 0 bits before the decimal place & 31 bits after the decimal place. So we’ll call it sFract0.31. Or simply fract0.31 because there are no unsigned fract math instructions on bfin, so we’d never have much use for uFract0.32.
All well and good to use fract0.31 for applying a fader/volume coefficient from 0 -> 0.9999999995343387. But what if you want the range of fract32 to represent, for example, -2.0 -> 1.9999999990686774, enabling a 6dB volume boost (this should be called fract1.30, because 1 bit before the decimal place, 30 bits after). Well the easy way to do this is:
fract32 a; //(an audio sample from -1.0 -> 1.0 in fract0.31)
fract32 b; //(a fader value from 0 -> 2.0 in fract1.30)
fract32 c; //(an audio sample in fract0.31)
c = shl_fr1x32(mult_fr1x32x32(a, b), 1); // (an audio sample from -1.0 -> 1.0 in fract0.31)
and, if we wanted a maximum of 24dB gain (which means multiply by 16, or shift-left 4), then
fract32 a; //(an audio sample from -1.0 -> 1.0 in fract0.31)
fract32 b; //(a fader value from 0 -> 2.0 in fract4.27)
fract32 c; //(an audio sample in fract0.31)
c = shl_fr1x32(mult_fr1x32x32(a, b), 4); // (an audio sample from -1.0 -> 1.0 in fract0.31)
So introduce the notion of a ‘fractional type’, taking inspiration from the fract math language feature in Haskell-ish HDL ClaSH. The three cases discussed so far have type signature:
fract0.31 * fract0.31 -> fract 0.31 (requires no shifting because fract primitives are in the same base)
fract1.30 * fract0.31 -> fract 0.31 (requires shift-left 1 because there’s 1 bit left of the decimal point on inputs)
fract4.27 * fract0.31 -> fract 0.31 (requires shift-left 4 because there are 4 bits left of the decimal point on inputs)
There’s a pattern here! Convince yourself the following is correct:
fract4.27 * fract6.25 -> fract 2.29 (requires shift-left 8)
fract4.27 * fract6.25 -> fract 10.21 (requires shift-left -2 or shift-right 2)
Same rule applies for the multiplication primitive mult_fr1x32, which multiplies two fract16 values to give a fract32. It has the fundamental type-signature:
fract0.15 * fract0.15 -> fract0.31 (without shift)
so, for example:
fract0.15 * fract5.10 -> fract0.31 (requires shift-left 5)
So for all these cases (apart from the one example with negative shift-left) there is a quantisation error the number of post-shift bits. In order to do mult_fr1x32x32 with maximum possible precision:
fract32 mult_fr1x32x32_autoscaling (fract32 a, fract32 b, short radix_imbalance) {
a_radix = norm_fr1x32(a); // returns the number of bits shift_left to ‘normalise’ a to fill all available bits
b_radix = norm_fr1x32(a); // ditto for b
return shl_fr1x32(shl_fr1x32(a, a_radix), shl_fr1x32(b, b_radix), radix_imbalance - a_radix - b_radix);
}
fract32 mult_fr1x32_autoscaling (fract16 a, fract16 b, short radix_imbalance) {
a_radix = norm_fr1x16(a); // returns the number of bits shift_left to ‘normalise’ a to fill all available bits
b_radix = norm_fr1x16(a); // ditto for b
return shl_fr1x32(shl_fr1x16(a, a_radix), shl_fr1x16(b, b_radix), radix_imbalance - a_radix - b_radix);
}
