0

I need to implement an RMS calculations of sine wave in MCU (microcontroller, resource constrained). MCU lacks FPU (floating point unit), so I would prefer to stay in integer realm. Captures are discrete via 10 bit ADC.

Looking for a solution, I've found this great solution here by Edgar Bonet: https://stackoverflow.com/a/28812301/8264292

Seems like it completely suits my needs. But I have some questions.

  1. Input are mains 230 VAC, 50 Hz. It's transformed & offset by hardware means to become 0-1V (peak to peak) sine wave which I can capture with ADC getting 0-1023 readings. Hardware are calibrated so that 260 VRMS (i.e. about -368:+368 peak to peak) input becomes 0-1V peak output. How can I "restore" back original wave RMS value providing I want to stay in integer realm too? Units can vary, mV will do fine also. My first guess was subtracting 512 from the input sample (DC offset) and later doing this "magic" shift as in Edgar Bonet answer. But I've realized it's wrong because DC offset aren't fixed. Instead it's biased to start from 0V. I.e. 130 VAC input would produce 0-500 mV peak to peak output (not 250-750 mV which would've worked so far). With real RMS to subtract the DC offset I need to subtract squared sum of samples from the sum of squares. Like in this formula: RMS formula

So I've ended up with following function:

#define INITIAL 512
#define SAMPLES 1024
#define MAX_V 368UL // Maximum input peak in V ( 260*sqrt(2) )
/* K is defined based on equation, where 64 = 2^6,
 * i.e. 6 bits to add to 10-bit ADC to make it 16-bit 
 * and double it for whole range in -peak to +peak
 */
#define K (MAX_V*64*2)

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

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

...
// Somewhere in a loop
getSample(&sample);
rms = rms_filter(sample);
...
// After getting at least N samples (SAMPLES * X?)
uint16_t vrms = (uint32_t)(rms*K) >> 16;
printf("Converted Vrms = %d V\r\n", vrms);

Does it looks fine? Or am I doing something wrong like this?

  1. How does SAMPLES (window size?) number relates to F (50Hz) and my ADC capture rate (samples per second)? I.e. how much real samples do I need to feed to rms_filter() before I can get real RMS value providing my capture speed are X sps? At least how to evaluate required minimum N of samples?
NStorm
  • 140
  • 7

1 Answers1

2

I did not test your code, but it looks to me like it should work fine. Personally, I would not have implemented the function this way. I would instead have removed the DC part of the signal before trying to compute the RMS value. The DC part can be estimated by sending the raw signal through a low pass filter. In pseudo-code this would be

rms = sqrt(low_pass(square(x - low_pass(x))))

whereas what you wrote is basically

rms = sqrt(low_pass(square(x)) - square(low_pass(x)))

It shouldn't really make much of a difference though. The first formula, however, spares you a multiplication. Also, by removing the DC component before computing the square, you end up multiplying smaller numbers, which may help in allocating bits for the fixed-point implementation.

In any case, I recommend you test the filter on your computer with synthetic data before committing it to the MCU.

How does SAMPLES (window size?) number relates to F (50Hz) and my ADC capture rate (samples per second)?

The constant SAMPLES controls the cut-off frequency of the low pass filters. This cut-off should be small enough to almost completely remove the 50 Hz part of the signal. On the other hand, if the mains supply is not completely stable, the quantity you are measuring will slowly vary with time, and you may want your cut-off to be high enough to capture those variations.

The transfer function of these single-pole low-pass filters is

H(z) = z / (SAMPLES * z + 1 − SAMPLES)

where

  • z = exp(i 2 π f / f₀),
  • i is the imaginary unit,
  • f is the signal frequency and
  • f₀ is the sampling frequency

If f₀ ≫ f (which is desirable for a good sampling), you can approximate this by the analog filter:

H(s) = 1/(1 + SAMPLES * s / f₀)

where s = i2πf and the cut-off frequency is f₀/(2π*SAMPLES). The gain at f = 50 Hz is then

1/sqrt(1 + (2π * SAMPLES * f/f₀)²)

The relevant parameter here is (SAMPLES * f/f₀), which is the number of periods of the 50 Hz signal that fit inside your sampling window. If you fit one period, you are letting about 15% of the signal through the filter. Half as much if you fit two periods, etc.

You could get perfect rejection of the 50 Hz signal if you design a filter with a notch at that particular frequency. If you don't want to dig into digital filter design theory, the simplest such filter may be a simple moving average that averages over a period of exactly 20 ms. This has a non trivial cost in RAM though, as you have to keep a full 20 ms worth of samples in a circular buffer.

Edgar Bonet
  • 3,416
  • 15
  • 18