Let's Design and Build a (mostly) Digital Theremin!

"Well... I should probably apply it to my calculation of "log(x)/log(2 to the power of (1/12))" in Open Theremin with MIDI."  - Mr_Dham

With log(x) / log(2^(1/12)) the denominator is a constant, so it is probably evaluated at compile time rather than run time.  And division by a constant is multiplication by its inverse, so the otherwise time consuming division is likely just a multiplication.  In the Hive math libraries I've only implemented log2 and exp2 (i.e. nothing based on the constant e) as they are the easiest and fastest to calculate in binary, and are the most useful.  On a 32 bit machine, log2 can actually be accomplished using only shifts and subtracts to produce a result within 0.3%, and on Hive this takes 9 cycles max.  Using a simple change of variable and a second order polynomial gets the error down to +/-0.025% (for larger integer inputs) and uses 13 cycles max.  These are all integer input, producing a 5.27 fixed point output.

That said, I don't see that log2, of ANY precision, is used anywhere in the D-Lev code base.  The field linearization maths also make the pitch and volume numbers linear, so there is just exp2 at the end of it all for our conveniently doubly logarithmic (i.e. pitch and volume) ears.  For the exp2 the code uses low quality (+/-0.0005% error, 4th order poly, 18 cycles) and very low quality (+/-0.33% error, 2nd order poly, 12 cycles).  The desired quality is selected by calling the specific subroutine, I wonder if there is any way to explicitly control this in the OpenTheremin code?  Higher level languages tend to get in the way of one micromanaging the maths and value representations themselves.

There is a rather trivial non-polynomial version of exp2 as well that only requires 7 cycles, but the error is 6.15% for larger inputs.  Clearly, log2 is somewhat more amenable to calculation than exp2.

Floats add a ton of overhead too, so I've very deliberately and laboriously avoided them whenever possible, sticking to ints and fracs (<1) though that can get pretty tricky.  Put enough elbow grease into it up front and you can put a man on the moon with a pocket calculator! :-)

For a good example of elbow grease mathematics and ham radio sdr have a look at (not my project but i use it):

https://github.com/threeme3/usdx

73 de Chris DL6SEZ

"https://github.com/threeme3/usdx"  - chrisbei

Thanks, that's an impressive project Chris!  I often wonder how SDR techniques might be applied to the digital Theremin.  The Theremin is sort of both a transmitter and receiver, but not at a fixed frequency, and with no modulation (beyond possible dither / spreading).  "Upside down" Theremin designs utilize a fixed frequency and then measure the phase offset introduced by the hand, which conveniently maximizes the far field response, but higher Q coils might not be entirely amenable to this scheme as the phase change and antenna voltage may poop out prematurely.

In the Hive math libraries I've only implemented log2 and exp2 (i.e. nothing based on the constant e) as they are the easiest and fastest to calculate in binary, and are the most useful. - Dewster

Yes, I agree, I would have initially searched for a log2. Any log can be achieved with a constant multiplication then I imagine that natural log from math.h is just log2(x)/log2(e).

I don't see that log2, of ANY precision, is used anywhere in the D-Lev code base - Dewster
Yes, if D-lev's working value is directly in note number (linear) then you don't need it, that's all the beauty of it.

I need some precision in the log because I must convert a frequency (exponential) into a note number (linear) with some resolution between two notes. And as I advertise that you can play theremin's sound along an external synth through MIDI, any approximation could push me out of tune. 
That's why I rather speak about a log 2^(1/12) (which is a "log note" if log2 is a "log octave").

Open Theremin V4 has a log2, somewhere in the code (to calculate a pitch CV), with linear interpolation between 1 and 2. But the error makes that a G is taken for a F#. It not a problem as long as you don't want to play a external synth along open theremin sound.

Shifts and subtracts will bring me the range [1; 2] (an octave), 12 more multiplications and subtracts will bring me in [1; 2^1/12], then linear interpolation can do between notes. 24 steps + 1 linear interpolation would give me the precision I need. Notes (C, C#, ...) would have almost no error which is nice. But I imagine it still takes some CPU load.

3rd degree Polynomial approximation of log2 in [1; 2] doesn't distribute precision as expected around 2. I probably need to consider a variable change or find a trick. I need to reopen my math lessons... (27 years after) but I didn't say my last word.

My MIDI data shows that the open theremin sends a message every 4 ms as expected by code construction. Then I have no problem. It ain't broke, don't fix it they say... But the uglyness of calculating in float... and that redundancy with pitch CV inaccurate calculation...

Vincent, I posted my Hive algorithms spreadsheet if you want to take a look at it.  It explains the change of variable too:

  https://d-lev.com/research/HIVE_algorithms_2018-05-15.xls

Also, for fitting data to polynomials, this site comes in pretty handy:

  http://polynomialregression.drque.net/online.php

Wow!

"for fitting data to polynomials, this site comes in pretty handy" - dewster
I won't count the number of sheet of paper I wrote down yesterday, all this ending with tons of error accumulation. One problem fixed, Thanks !

Hive algorithms spreadsheet - dewster
Very nice work and very clear explanations. 
Change of variable: got it, I need to bring my calculation where the polynomial function is precise and efficient (it talks to me and my mathematics memories). 

 
Shifts and subtracts will bring me the range [1; 2] (an octave), 12 more multiplications and subtracts will bring me in [1; 2^1/12], then linear interpolation can do between notes. 24 steps + 1 linear interpolation would give me the precision I need. Notes (C, C#, ...) would have almost no error which is nice. But I imagine it still takes some CPU load. - Me...
More "dichotomic" approach requires only 4 steps to separate exponent and mantissa over 11 otaves for log2, and 4 [color=#000000]more steps to separate exponent and mantissa over 12 notes for log2^1/12 ("log note"). More optimistic on cpu load,... but a risk of corner case. [/color]

I have two options now (polynomial and log note), let's try both and compare...

Thanks for the hints.

"Change of variable: got it, I need to bring my calculation where the polynomial function is precise and efficient (it talks to me and my mathematics memories)."  - Mr_Dham

Many functions don't converge (give decreasing error with increasing polynomial order) without a change of variable.  With it, log2 converges pretty well because it's fairly linear (over the poly interval) in the first place - but 32 bits of precision require 11th order or thereabouts.  Sine and Cosine converge like crazy because they are even / odd (poly factors are applied to x^2n rather than x^n) and they are basically very close to a squared function, though not spectrally pure enough to use just one term - the 16 or so bits of precision as used in the D-Lev oscillators require a 3rd order poly.  I "wasted a lot of time" (i.e. had a lot of fun!) coming up with 32 bit precision everything only to end up not really using any of it, though it helps to know the rules of the game before one goes off breaking them.  If I had to do it over again I'd try to use golang and bessel functions rather than doing things manually and laboriously in a spreadsheet.  Writing the Hive math library was a fascinating intellectual experience, it provided a lot of DSP groundwork and was an excellent proving ground for the Hive opcode set.

For your log2, a possible error target might be 1 cent, which is 2^(1/1200) = 1.0005778 or 0.06%.  A second order polynomial here should give you +/-0.025% (for large inputs).

I think I am on the good way, thanks !

For your log2, a possible error target might be 1 cent, which is 2^(1/1200) = 1.0005778 or 0.06%.  A second order polynomial here should give you +/-0.025% (for large inputs). Dewster

It sounds more than reasonable, the 16 bit resolution of what stands for frequency makes that theres is 55 steps between A4 and A4#. and it gets even worse as we reduce the frequency... 
Log float is definitely an overkill !