waves: tuning quantization (or: better 16.16 math intrinsics)

[catfact]

'hz' and 'ratio' parameters are both 16.16 values.

the phase increment for each oscillator is calculated by multiplying base hz and tuning ratio to get final frequency:
https://github.com/monome/aleph/blob/dev/modules/waves/osc_waves.c#L62

and then applies additional operations to convert from frequency to phase increment:
https://github.com/monome/aleph/blob/dev/modules/waves/osc_waves.c#L62

the first multiply naively uses the non-blackfin-specific fixed-point arithmetic lib (this one, i think):
https://github.com/monome/aleph/blob/dev/bfin_lib/src/libfixmath/fix16.c#L122

all of this is slow, and ultimately results in the operands being truncated somewhere along the way, effectively quantizing the values of the 'tuning' parameter...

in general, more optimized 16.16 operations would be useful across the board for blackfin modules

[ghost]

Ha - yeah I copy-pasted this same code for new fmsynth, monosynth & possibly also voder. Let's sort it out!

Some possible tricks that can help here:

• phase increment == frequency (in units where nyquist frequency == 1) so that's one pointless conversion out the way.
• In order for this to work, the param scalar code needs to know the sample rate.

according to what I learnt figuring out how to do super-slo-mo log slews, 16.16 multiply with maximum possible accuracy can be done like this:

fract32 mult_16_16(fract32 a, fract32 b) {
  char a_radix = norm_fr1x32(a);
  char b_radix = norm_fr1x32(b);
  return shl_fr1x32(mult_fr1x32x32(shl_fr1x32(a, a_radix), shl_fr1x32(b, b_radix)),
                               15 - a_radix - b_radix);
}

I've forgotten your recipe to look at the asm for a function from the elf or whatever, so don't really have any way to tell whether that is more or less efficient on bfin than the generic 16.16 mult!

It's kind of brute force. Obv if you already know a & b don't vary much, this can be rewritten as static inline, specifiying a_radix & b_radix at compile time. Saves a couple cycles checking a_radix & b_radix at runtime...

[ghost]

EDIT:
Thought harder about this, then re-wrote the code in this post without the pre/post-shift (see above fract 16.16 for an example of that). 27 bits number is plenty of resolution for any sane phase_increment variables! Also I just revised the notation of previous edits to make things less confusing. We are always talking about signed numbers with fract arithmetic on bfin. So the primitive 32 bit type should really be called sFract0.31 (because it's signed, 0 bits in front of the decimal place, 31 bits after). We can shorten to fract0.31, because there is no uFract math.

pretty sure '16.16' (aka fract15.16) is just not quite the right thing here...

Look at the 'tuning' parameter of grains' pitchshift. That is fract23.8 - maybe it's a little too coarse? integer tuning from -128 -> 127 is a bit over the top now I come to think about it.

Much better would be a fract base that gives a scaling range from -32 -> 32. that divides the first musical octave into 1024 distinct pitches. ~1% semitone accuracy & 5 octaves upwards sounds like enough to me. That is fract 5.10.

This requires a recipe for fract multiply with type signature:
fract0.15 * fract5.10 -> fract0.31.

So organise the BEES tune param to send bfin a value from 0 -> S16_MAX & displays -32 -> 32 on the LCD, where S16_MAX represents a factor of 32 in fract5.10.

Then organise the frequency param scalar so that it writes Hz to screen (in linear increments), whilst sending bfin a value from 0 -> S16_MAX, where S16_MAX represents nyquist frequency in fract0.15.

Now the calculation of phase increment goes:

fract32 calc_phase_increment(fract16 freq, fract16 tune) {
 return shl_fr1x32(mult_fr1x32(freq, tune), 5);
}

Don't think it's worth worry about 5 lost bits of precision. The quantisation error due to the left shift corresponds to 3.5e-4 Hz...

If Hz/tune have to be in log increments then the above is not quite right, you'd have to send bfin values from 0 -> FR32_MAX, then (again with quantisation error of 5 bits):

fract32 calc_phase_increment(fract32 freq, fract32 tune) {
  return shl_fr1x32(mult_fr1x32x32(freq, tune), 5);
}

the type signature of this is:
fract5.26 * fract0.31 -> fract0.31

[catfact]

thanks! you are totally right, norm_fr1x32 to get appropriate radix is the key here. i think your mult_16_16 example looks like a more effecient drop-in replacement for the version in fixmath (which was always only meant to be a placeholder.)

interesting point about reducing the depth for frequency multipliers. i have to think about it a bit more.

one note: in the beginning, i indeed deliberately structured things so that bees param scaling didn't have to know the sample rate. but i think that has gone out the window anyways.

[ghost]

IIRC, blackfin sample rates can only be set to (2^n) * 48kHz?

So how about this hack? Have BEES believe every module runs at 48kHz (set with a single pre-processor #define obviously to make the code portable in future to a platform supporting different SR). Then the DSP module is responsible for right-shifting that frequency by 1/2 bits if it's samplerate is in fact 96/192 kHz ...

IIRC read somewhere that the mult_fr1x32 instruction actually runs in less cycles than mult_fr1x32x32. Again, dunno how to test that!

[catfact]

I'm actually arguing norm_fr1x32 is unnecessary for this calculation.

ah yes, i know, i mean for the general 16.16 x 16.16 case, it's the trick that libfixmath doesn't know about.

[ghost]

So re-reading & re-thinking again about this I do still think reducing depth for frequency multipliers as described above is the 'right thing'. But requires me to properly grok param scalar code, which I don't feel like doing right now.

On the other hand, am I right in thinking there's nothing very badly wrong with usage of general 16.16 (when improved by fract arithmetic discussed above)?

I've run with that for now, replacing all fix16_mul with the new fix16_mul_fract in all blackfin modules & dsp code (untested but pushed to my dev branch).