Help avoiding floating point math

Question

I'm using a Texas Instruments DRV2605 haptic controller in my project to drive a LRA vibration motor. I would like to calculate the LRA rated voltage during run time using the formula below as provided in the datasheet:

I'm successfully calculating the rated voltage, but the calculation requires floating point calculations and I would like to avoid that because I'm low on available flash space on my ATTiny85 micro controller. Here is how I'm currently calculating the rated voltage value:

uint8_t LRA_closed_loop_rated_voltage( uint16_t voltage ) {
   float RatedVoltage;
   uint16_t SampleTime = 300;
   uint16_t frequency = 175;

   RatedVoltage = ( (float) voltage ) * sqrt( 1 - ( ( 4 * SampleTime + 300 ) * 0.000001 ) * frequency );
   RatedVoltage = ( RatedVoltage * 255 ) / 5300;

   return RatedVoltage;
}

The function above is just a simplified version showing the calculation and the SampleTime and frequency variables are set to their defaults values here. The voltage argument is the actuator's rated voltage in millivolt.

EDIT: Even though SampleTime and frequency have fixed values in the function above, their values will be determined by reading the relevant registers in the final implementation. SampleTime may have one of the following values: 150, 200, 250, 300 and the frequency will depend on the manufacturers specifications for the actuator.

It was easy to overcome the floating point issue on the other calculations by simply scaling up the values by 1000 (fixed point arithmetic):

ODClamp = ( Vpeak * 255) / 5.6 /* Vpeak in volts */

TO:

ODClamp = ( Vpeak * 255) / 5600 /* Vpeak in millivolts */

How can I do the calculation in question without having to use floating point?

Answer 1

As I said in my comment, it is all about computing an approximation of sqrt(1-x). Assuming the typical resonance frequency is between 175 and 235 Hz, then we have

(4 × 150 µs + 300 µs) 175 Hz ≤ x ≤ (4 × 300 µs + 300 µs) 235 Hz

i.e.

0.1575 ≤ x ≤ 0.3525

This range is small enough to try a simple polynomial approximation. A linear approximation is just a bit too crude for 8 bit resolution. Playing with quadratic polynomials, I found this approximation which is good to 3.1e-5 in the relevant range:

√(1−x) ≈ 0.9983795 − x (0.480515 + x 0.1955)

This can be easily implemented in fixed point by scaling everything by a factor 2¹⁶ = 65536. One has to take care that any product of two such fixed point numbers must be computed in 32 bits, and the result subsequently divided by 65536. Which gives the following implementation:

/*
 * Compute a fixed point approximation of f(x) = sqrt(1-x).
 * Valid for x in [.1575, .3525].
 * Argument and result scaled by 2^16.
 */
uint16_t f(uint16_t x)
{
    return 65430 - x * (31491 + x * 12812UL / 65536) / 65536;
}