Puzzle: extended precision division.

Lawson

I've built some excessively precise temperature sensors into a controller I'm building. They use 24-bit ADCs to read the resistance of thermistors arranged in a voltage divider using an ultra-stable 5Kohm pull-up resistor to a common reference. In my low level code I then sum up 128 samples from this ADC and left justify the result as a 32-bit signed integer. (though only positive values are valid.) The equation to convert ADC values to resistance is shown below.

resistance[ohms] / 5000[ohms] = adc / ( full_scale - adc )
full_scale = $7FFF_FF00
adc = 678348544 (at ~25C)

I'm also expressing resistance in 16.16 fixed point format, so there's a factor of 2^16 thrown in too. (I use an interpolated table lookup to convert to temperature) In my data I've seen at least 16 stable bits (ends up as ~+-0.001 [C] precision) and I should theoretically be able to get ~20 stable bits. (need to track down some circuit noise for that) The precision of the calculation should extend out to at least 24-bits because those extra unstable bits still contain useful information I can get at by averaging further.

I'm currently using FloatMath.spin and floating point. The 23-bit precision of floating point should be enough. I'd like to keep this all in spin, but eventually the control code will end up doing ~100 multiplications or divisions at 32Hz so speed is of some importance.

Lawson

P.S. It just crossed my mind to search the OBEX. Found some useful looking stuff, so suggestions for favorite objects are also welcome.

Dave Hein

Assuming you continue to use the adc(25) number that you gave in your post, just need to multiply times 957_961_678, and keep the high 32 bits. In Spin you would do

temp := adc ** 957_961_678

This will give you the temperature as a 16.16 number.

According to your formula, the tempurature is given by 5000/(full - adc(25)). (2^32)/(full - adc(25)) equals 2.923. Multiply this time 5000 to get 14617.335178 Scaling this by 2^16 gives 957_961_678.

Phil Pilgrim (PhiPi)

Um, Dave, I don't think that will work. His formula for resistance vs. adc is nonlinear: R = 5000 adc / (FS - adc)

-Phil

Dave Hein

OK, I misunderstood. The (FS - adc) can be done in fixed-point. Maybe there's a quick way to do 1/x without resorting to floating point. Anyhow, using floating point would be the easiest. More speed and precision could be obtained by writing a 64-bit/32-bit divide routine in PASM.

Phil Pilgrim (PhiPi)

Come to think of it, my umath object in the OBEX has a 32*32 (=64) / 32 = 32 combo mul/div method that might work.

-Phil

Tracy Allen

I think you can do that something like the following, using binary division:

PUB solve4R(adc, K, R0) | x,z  ' enter with adc reading, K=full scale $7FFF_FF00, R0=5000 ohms
  '  solve for result = R0 * adc / (K - adc)
  x  := K - adc    ' the denominator, will be less than adc when adc is greater than 1/2 scale
  z := adc / x    ' integer part, will be non-zero when adc is greater than 1/2 scale
  result := binNormal(adc // x, x,  31, 1)   ' normalize remainder ((adc//x) / x) to a denominator of 3^31
  result := 5000 * z + (2 * 5000 ** result)   ' combine integer and fractional parts.   
                                                           

PRI binNormal(y, x, b, r) : f                  ' calculate f = y/x * 2^b, round off if r=1, truncate if r=0
' b is number of bits
' enter with y,x: {x > y, x < 2^31, y <= 2^31}
' exit with f: f/(2^b) =<  y/x =< (f+1) / (2^b)
' that is, f / 2^b is the closest appoximation to the original fraction for that b.
' r=0 truncates low, r=1 rounds up or down.
  f~
  repeat b
    y <<= 1
    f <<= 1
    if y => x    '
      y -= x
      f++
  if r == 1 and (y << 1) => x
    f++

I think Phil's Umath object would work too, the double precision calculation of R=5000*adc/(K-adc):

R := multdiv(5000, adc, K-adc)

Since you're doing a table lookup anyway, why not just go directly from adc to R in the table?

Lawson

I did see Phil's Umath object. It would fit the bill, but I'm wondering why it does a long division instead of something using spin's "/" and "//" operators?

I'm doing the table interpolation separate from the resistance conversion because I may substitute in one of the standard thermistor calibration equations in the future. Also, the linear interpolation equation in my table lookup suffers from this same precision challenge, so I'd have to solve it anyway.

Lawson

Tracy Allen

A rep from Vishay Precision came by our shop a couple of weeks ago and showed off eye-opening demos of resistor stability. TC, overload, ESD. We had purchased a few of their 0.01% 0.2ppm bulk foil resistors to use for calibration purposes. Is that the kind of reference R you are using? Great resistors, but we won't be purchasing them in quantity at $15 a pop!

I have heard that Steinhart-Hart can give accuracy to better than 0.01 K. But doesn't it take individual calibration to achieve that kind of accuracy? A live 3-point or regression?

A long division process is necessary in general to find the double precision result when you are dealing with numbers that push up against $7FFFFFFF. It might be possible to do something with / and // alone in a couple of steps if there were more "headroom" available. The division algorithm that I suggested is not really double precision. It takes a proper fraction and converts it to a number that can be used with the ** operator.

Lawson

Tracy Allen wrote: »

A rep from Vishay Precision came by our shop a couple of weeks ago and showed off eye-opening demos of resistor stability. TC, overload, ESD. We had purchased a few of their 0.01% 0.2ppm bulk foil resistors to use for calibration purposes. Is that the kind of reference R you are using? Great resistors, but we won't be purchasing them in quantity at $15 a pop!

I have heard that Steinhart-Hart can give accuracy to better than 0.01 K. But doesn't it take individual calibration to achieve that kind of accuracy? A live 3-point or regression?

A long division process is necessary in general to find the double precision result when you are dealing with numbers that push up against $7FFFFFFF. It might be possible to do something with / and // alone in a couple of steps if there were more "headroom" available. The division algorithm that I suggested is not really double precision. It takes a proper fraction and converts it to a number that can be used with the ** operator.

Yea I'm using 3 resistors similar to those Vishay ones. (one for each thermistor) Don't remember the part number off hand but I remember 0.05% and 0.2ppm specifications and "only" $5.50 each. (may well have been Vishay parts too)

I've been very careful to say that I need crazy precision. I need the precision and low noise so I can make the thermal controller on this laser as stiff and stable as practical without making the output too noisy. High accuracy will be nice but ultimately I'll get that when I lock the laser to a molecular absorption cell. So within 0.2 [K] is plenty for now. Hence why I'm using a lookup table pulled from the thermistor's data-sheet, but interpolating it out to 23+ bits. (and thinking I might have to switch to Steinhart-Hart to avoid the discontinuous derivatives of temperature due to interpolation)

For long division I was thinking of an extension of base 10 long division. Only in this case it would be base 2^16 and work on the four 16-bit "symbols" that make up the 64-bit numerator. (or base 2^21 using 3 symbols and seperately handling the sign bit?)

Lawson

Tracy Allen

Is that because the modes of the laser are so temperature dependent that 0.001 °C of stability is required? Intriguing! If the goal is a crazy precision set point based on feedback, there would be no need to convert from adc counts to resistance to temperature. The set point could be in terms of the raw count. Is it a ratiometric converter?