Moving average filter implementation in the floating point [closed]

Question

I have following C++ implementation of the moving average digital filter.

template<uint32_t WINDOW_SIZE>
class MovingAverage{
    public:
        MovingAverage(){
            for(uint32_t elem = 0; elem < WINDOW_SIZE; elem++){
                buffer[elem] = 0;
            }
            input_value = 0;
            sum = 0;
            average = 0;
            index = 0;
        }
        
        void setInput(float input){
            input_value = input;
        }
        
        void calculate(){
            sum -= buffer[index];
            buffer[index] = input_value;
            sum += buffer[index];
            average = sum/WINDOW_SIZE;
            index = (index + 1) % WINDOW_SIZE;
        }
        
        float getOutput(){
            return average;
        }
        
    private:
        float buffer[WINDOW_SIZE];
        float input_value;
        float sum;
        float average;
        uint32_t index;
};

I have encountered a weird behavior on the ARM Cortex A9 platform. In case I at first pass the values coming from the adc (those values fluctuates about 514) into the moving average with window size equal to 32 samples and then I force the value 128 to the filter then the filter outputs value about 126 instead of 128. At this moment the buffer contains value 128 in all its items but the sum contains different value than the expected 32*128. In case I run the filter algorithm on my PC it works as expected. Does anybody have any idea what could be the reason of this behavior?

Answer 1

Well, 514 is close enough to a power of two, when a floating point number crosses a power of two it looses one digit accuracy.

The floating point number will be stored with 24-bit mantissa.

Your sum will be around $2^{15}$, when the sum is less than $2^{15}$, your input will be accumulated with precision $2^{14-24}=2^{-10}$, if the sum accumulator is greater or equal $2^{15}$, then the exponent changes from $14$ to $15$ and the the bit $2^{-10}$ is rounded of, doing this systematically for a few thousands of iterations will accumulate to the error you are observing.

Let $E=2^{-11}$, and $-1/2 \le {\tt rand()} < 1/2$ a uniform random number, and avg.step set the input and calculate, and return the difference of the accumulated sum and the direct sum of the values in the buffer. Repeat the following.

    delta = E if s < 512 else -E
    s = avg.step(512 + delta)
    s = avg.step(512 + (-0.5 <= rand() <0.5))

The accumulated errors evolves like this

You should not see the same behavior if you are not crossing a power of two boundary. You could sacrifice a few bits precision to avoid this accumulated errors.

A quick solution is initialize your sum to a different value, then you subtract that value when you calculate the average. For signed inputs sum = 3 * WINDOW_SIZE * MAX_INPUT, then for inputs between -MAX_INPUT and MAX_INPUT the sum will stay in a range of numbers represented with the same exponent. For unsigned MAX_INPUT * WINDOW_SIZE gives you one additional bit accuracy. This is effectively making your computations fixed point.