Any Faster RMS Value Calculation in C?

Question

I am writing a software for a small 8-bit microcontroller in C. Part of the code is to read the ADC value of a current transformer (ZCT), and then calculate the RMS value. The current flowing through the ZCT is sinusoidal but it can be distorted. My code as follow:

float       adc_value, inst_current;
float       acc_load_current;           // accumulator = (I1*I1 + I2*I2 + ... + In*In)
double      rms_current;

// Calculate the real instantanous value from the ADC reading
inst_current = (adc_value/1024)*2.5;    // 10bit ADC, Voltage ref. 2.5V, so formula is: x=(adc/1024)*2.5V                           

// Update the RMS value with the new instananous value:
// Substract 1 sample from the accumulator (sample size is 512, so divide accumulator by 512 and substract it from the accumulator)
acc_load_current -= (acc_load_current / 512);       
inst_current *= inst_current;          // square the instantanous current
acc_load_current += inst_current;      // Add it to the accumulator

rms_current = (acc_load_current / 512);  // Get the mean square value. (sample size is 512)
rms_current = sqrt(rms_current);         // Get RMS value

// Now the < rms_current > is the real RMS current

However, it has many floating point calculations. This adds a large burden to my small MCU. And I found that the sqrt() function does not work in my compiler.

Is there any code that could run faster?

Answer 1

1. divisions/multiplications by power of 2

   can be done by changing the exponent only via bit mask operations and +,- so mask/extract the exponent to integer value then apply bias. After that add/sub the value log2(operand) and encode back to your double value

2. sqrt

   how fast and accurate it should be? You can use binary search on fixed point or use sqrt(x)=pow(x,0.5)=exp2(0.5*log2(x)). Again on fixed point it is quite easy to implement. You can temporarily make double a fixed point by taking the mantissa and bit shift it to some known exponent around your used values + handle the offset or to 2^0 if you have enough bits ...

   compute sqrt and then store back to double. If you do not need too big precision then you can stay on operand exponent and do the binary search directly on mantissa only.

[edit1] binary search in C++

//---------------------------------------------------------------------------
double f64_sqrt(double x)
    {
    const int h=1;      // may be platform dependent MSB/LSB order
    const int l=0;
    DWORD b;            // bit mask
    int e;              // exponent
    union               // semi result
        {
        double f;       // 64bit floating point
        DWORD u[2];     // 2x32 bit uint
        } y;
    // fabs
    y.f=x;
    y.u[h]&=0x7FFFFFFF; x=y.f;
    // set safe exponent (~ abs half)
    e=((y.u[h]>>20)&0x07FF)-1023;
    e/=2;               // here can use bit shift with copying MSB ...
    y.u[h]=((e+1023)&0x07FF)<<20;
    // mantisa=0
    y.u[l] =0x00000000;
    y.u[h]&=0xFFF00000;
    // correct exponent
    if (y.f*y.f>x) { e--; y.u[h]=((e+1023)&0x07FF)<<20; }
    // binary search
    for (b =0x00080000;b;b>>=1) { y.u[h]|=b; if (y.f*y.f>x) y.u[h]^=b; }
    for (b =0x80000000;b;b>>=1) { y.u[l]|=b; if (y.f*y.f>x) y.u[l]^=b; }
    return y.f;
    }
//---------------------------------------------------------------------------

it returns sqrt(abs(x)) the results match "math.h" implementation from mine C++ IDE (BDS2006 Turbo C++). Exponent starts at half of the x value and is corrected by 1 for values x>1.0 if needed. The rest is pretty obvious

Code is in C++ but it is still not optimized it can be surely improved ... If your platform does not know DWORD use unsigned int instead. If your platform does not support 32 bit integer types then chunk it to 4 x 16 bit values or 8 x 8 bit values. If you have 64 bit then use single value to speed up the process

Do not forget to handle exponent also as 11 bit .... so for 8 bit registers only use 2 ... The FPU operations can be avoided if you multiply and compare just mantissas as integers. Also the multiplication itself is cumulative so you can use previous sub-result.

[notes]

For the bit positions see wiki double precision floating point format

Answer 2

To find the square root use the below app note from microchip for 8 bit controllers

Fast Integer Square Root

which is much fast and can find square root in just 9 loops

Answer 3

When you need to get faster on an processor that lacks an FPU, your best bet is to replace floating point calculations with fixed point. Combine this with joop's suggestion (a one Newton-Raphson sqrt) and you get something like this:

#define INITIAL 512  /* Initial value of the filter memory. */
#define SAMPLES 512

uint16_t rms_filter(uint16_t sample)
{
    static uint16_t rms = INITIAL;
    static uint32_t sum_squares = 1UL * SAMPLES * INITIAL * INITIAL;

    sum_squares -= sum_squares / SAMPLES;
    sum_squares += (uint32_t) sample * sample;
    if (rms == 0) rms = 1;    /* do not divide by zero */
    rms = (rms + sum_squares / SAMPLES / rms) / 2;
    return rms;
}

Just run your raw ADC samples through this filter. You may add a few bit-shifts here and there to get more resolution, but you have to be careful not to overflow your variables. And I doubt you really need the extra resolution.

The output of the filter is in the same unit as its input. In this case, it is the unit of your ADC: 2.5 V / 1024 ≈ 2.44 mV. If you can keep this unit in subsequent calculations, you will save cycles by avoiding unnecessary conversions. If you really need the value to be in volts (it may be an I/O requirement), then you will have to convert to floating point. If you want millivolts, you can stay in the integer realm:

uint16_t rms_in_mV = rms_filter(raw_sample) * 160000UL >> 16;

Answer 4

Since your sum-of-squares value acc_load_current does not vary very much between iterations, its square root will be almost constant. A Newton-Raphson sqrt() function normally converges in only a few iterations. By using one iteration per step, the computation is smeared out.

static double one_step_newton_raphson_sqrt(double val, double hint)
{
double probe;
if (hint <= 0) return val /2;
probe = val / hint;
return (probe+hint) /2;
}

static double      acc_load_current = 0.0;           // accumulator = (I1*I1 + I2*I2 + ... + In*In)
static double      rms_current = 1.0;
float       adc_value, inst_current;
double      tmp_rms_current;

// Calculate the real instantanous value from the ADC reading
inst_current = (adc_value/1024)*2.5;    // 10bit ADC, Voltage ref. 2.5V, so formula is: x=(adc/1024)*2.5V                           

// Update the RMS value with the new instananous value:
// Substract 1 sample from the accumulator (sample size is 512, so divide accumulator by 512 and substract it from the accumulator)
acc_load_current -= (acc_load_current / 512);
inst_current *= inst_current;          // square the instantanous current
acc_load_current += inst_current;      // Add it to the accumulator

tmp_rms_current = (acc_load_current / 512);
rms_current = one_step_newton_raphson_sqrt(tmp_rms_current, rms_current);         // Polish RMS value

// Now the <rms_current> is the APPROXIMATE RMS current

---

Notes:

• I changed some of the data types from float to double (which is normal on a general purpose machine/desktop) If double is very expensive on your microcomputer you could change them back.
• I also added static, because I did not know if the code was from a function or from a loop.
• I made the function static to force it to be inlined. If the compiler does not inline static functions, you should inline it manually.

Answer 5

Hopefully your project is for measuring "Big" AC voltages ( and not something like 9v motor control. ) If this happens to be the case then you can cheat because your error can then be within accepted limits..

Do all of the maths in integer, and use a simple lookup map for the sqrt operation. ( which you can pre-calculate at startup, you would only REALLY need about ~600 odd values if you are doing 3 phase..

Also this begs the question do you ACTUALLY need VAC RMS or some other measure of power ? (for example can you get away with a simple box mean algorythm? )