Questions conerning fixed point implementation of a moving window Newton-Raphson Method

Question

Here is the reference post I am drawing upon: Any Faster RMS Value Calculation in C?

#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;
}

1. Is the parameter sample already a fixed point value and is already scaled by S?

2. Is sum_squares being changed to a fixed point value in the following line?

static uint32_t sum_squares = 1UL * SAMPLES * INITIAL * INITIAL;

3. Is the following line to counteract the the squaring of sample above it? In addition, this is integer division which means the fractional portion will be truncated. Is this okay? Aren't we losing precision?

sum_squares / SAMPLES

4. If sum_squares is a fixed point value, shouldn't the 2 in the second to last line also be changed to a fixed point value?

Answer 1

1. Is the parameter sample already a fixed point value and is already scaled by S?

There is no evident fixed-point representation in the code shown other than the trivial fixed-point format of ordinary integers.

2. Is sum_squares being changed to a fixed point value in the following line?

static uint32_t sum_squares = 1UL * SAMPLES * INITIAL * INITIAL;

No. The filter is initialized as if prior samples had all had the value INITIAL. In this case, the square of one sample is INITIAL*INITIAL, and the sum of the squares of SAMPLES samples is SAMPLES * INITIAL * INITIAL. The 1UL is present to ensure the operands are converted to at least an unsigned long width before multiplication proceeds, to avoid some overflow.

This gives us a clue that there is no fixed-point representation here. A fixed-point representation essentially represents some value x with some scaled value sx, where s is a fixed scale chosen for the particular fixed-point format. (Often s is a power of two or ten, but any number can be used.) In a fixed-point format, if we multiplied the representative of x, sx, by the representative of y, sy, using integer arithmetic, we would get s²xy. But what we want is the representative of xy, which is sxy. Thus, we would have to divide by s. So, if there were fixed-point arithmetic here, we would expect to see something like INITIAL * INITIAL / SCALE. Since there is no division by a scale (or multiplication by the inverse of a scale), there is no fixed-point arithmetic.

3. Is the following line to counteract the the squaring of sample above it?

sum_squares / SAMPLES.

No, this line does not counteract the squaring. It calculates the mean of the squares, as part of the calculation of the square root of the mean of the squares. Since sum_squares has the sum of the squares and the number of them is SAMPLES, sum_squares / SAMPLES is the mean.

In addition, this is integer division which means the fractional portion will be truncated. Is this okay? Aren't we losing precision?

Some accuracy is lost. Whether it is okay or not depends on the application. If it is okay to get approximate results, and the errors caused by this truncation are not too great for the application, then it is okay. (By and large, signal processing people seem fine with losing lots of accuracy.)

4. If sum_squares is a fixed point value, shouldn't the 2 in the second to last line also be changed to a fixed point value?

No, the 2 is used to take an average of the earlier root-mean-square estimate and the new root-mean-square estimate.

The idea here is that rms contains the root-mean-square of the prior samples, or something close to it. Then, since sum_squares / SAMPLES is the mean of the squares of the most recently samples, dividing that new mean by the square root of the old mean (rms) gives something that is approximately the new root-mean-square. That is added to the old mean and averaged by dividing by two. They are playing fast and loose with the mathematics here, but that is what signal processing people do. Some analysis of how well this approximation tracks the true root-mean-square might be interesting, but I will not tackle it here.