How to implement a moving average in C without a buffer?

Question

Is it possible to implement a moving average in C without the need for a window of samples?

I've found that I can optimize a bit, by choosing a window size that's a power of two to allow for bit-shifting instead of dividing, but not needing a buffer would be nice. Is there a way to express a new moving average result only as a function of the old result and the new sample?

Define an example moving average, across a window of 4 samples to be:

ma <= (a + b + c + d) / 4

Add new sample e:

ma_new <= (a + b + c + d) / 4 - (a / 4) + (e / 4)
ma_new(ma, oldest_sample, new_sample) <= ma - (a / 4) + (e / 4)

Answer 1

A moving average can be implemented recursively, but for an exact computation of the moving average you have to remember the oldest input sample in the sum (i.e. the a in your example). For a length $N$ moving average you compute:

$$y[n]=\frac{1}{N}\sum_{k=n-N+1}^nx[k]\tag{1}$$

where $y[n]$ is the output signal and $x[n]$ is the input signal. Eq. (1) can be written recursively as

$$y[n]=y[n-1]+\frac{1}{N}(x[n]-x[n-N])\tag{2}$$

So you always need to remember the sample $x[n-N]$ in order to compute (2).

As pointed out by Conrad Turner, you can use an (infinitely long) exponential window instead, which allows you to compute the output only from the past output and the current input:

$$y[n]=\alpha x[n]+(1-\alpha)y[n-1]\tag{3}$$

but this is not a standard (unweighted) moving average but an exponentially weighted moving average, where samples further in the past get a smaller weight, but (at least in theory) you never forget anything (the weights just get smaller and smaller for samples far in the past).

Answer 2

What is wrong with a fading memory (exponential) moving average:

ma_new = alpha * new_sample + (1-alpha) * ma_old

Answer 3

code for (3) above , simple recursive filter $y(n)=αx(n) + (1-α)y(n-1)$. Where $α<1$ . The smaller $α$ the smoother the filter. I use it to smooth a data stream with no buffers.

uint32 filter(uint32 x, uint32 y)
{
    return (x + 9*y) / 10;  //( α=1/10 )
}

main()
{
  angle = filter(new_angle,angle);   // y(n) = filter( x(n), y(n-1) )
}

Answer 4

This method gives you an approximation of the moving average by basically assuming that the value of the sample window_size samples ago is equal to the previous moving average, which is updated every window_size samples.

It works well if your values are randomly distributed, but outliers will skew it more than the exact moving average.

previous_average = 0
total = 0
for count, sample in enumerate(samples):
    if count % window_size == 0: # Update previous_average every window_size samples
        previous_average = total/window_size
    total += sample
    if count > window_size:
        total -= previous_average
    current_average = total/window_size

With the 4 point example we can estimate the error:

$\begin{align} y_n = \frac{(a + b + c + d)}{4} &, y_{n+1} = \frac{(a+b+c+d)}{4} + \frac{e}{4} - \frac{a}{4} \\ y*_{n+1} &= \frac{(a + b + c + d) - y_{n} + e}{4} \\ &= \frac{(a+b+c+d)}{4} - \frac{(a+b+c+d)}{16} + \frac{e}{4} \\ &= \frac{3(a+b+c+d)}{16} + \frac{e}{4} \\ \end{align}$

$\begin{align} E &= y*_{n+1} - y_{n+1} \\ &= \frac{3a}{16} - \frac{a+b+c+d}{16} \end{align}$

I am lazily going to assume that with an arbitrary window size $W$

$\begin{align} E &= \frac{(W-1)a}{W^2} - \frac{y_n}{W} \\ \end{align}$

If $W$ is large enough, then:

$\begin{align} E &\approx \frac{a - y_n}{W} \end{align}$

Which is the residual of point $a$ divided by the total number of residuals...

This is just an approximation of the error for the first count % window_size == 0 iteration, and I'm already trying to go further than my maths abilities, but since we can start at any $x_n = a$ it suggests that as long as $W$ is sufficiently large or each $x_{n} - y_{n}$ is sufficiently small then this average will produce the same average residual as a fit with the exact formula.

I've attached an image from my test script showing the exact 100 point moving average compared to this method:

Answer 5

I implemented a moving average without individual item memory for a GPS tracking program I wrote.

I start with 1 sample and divide by 1 to get the current avg.

I then add anothe sample and divide by 2 to the the current avg.

This continues until I get to the length of the average.

Each time afterwards, I add in the new sample, get the average and remove that average from the total.

I am not a mathematician but this seemed like a good way to do it. I figured it would turn the stomach of a real math guy but, it turns out it is one of the accepted ways of doing it. And it works well. Just remember that the higher your "length" the slower it is following what you want to follow. That may not matter most of the time but when following satellites, if you are slow, the trail could be far from the actual position and it will look bad. You could have a gap between the sat and the trailing dots. I chose a length of 15 updated 6 times per minute to get adequate smoothing and not get too far from the actual sat position with the smoothed trail dots.

Answer 6

#define COEF 32

uint32_t filter(uint32_t raw)
{
    static uint64_t accum = 0;
    accum = accum - accum / COEF + raw;
    return accum/COEF;
}
```

Answer 7

Instigated by the szmoore answer 1

What is the best approximation to a moving average of length $W$ with an estimator?

$\mu_{k+1} = a \mu_{k} + b x_{k+1}$

With $a=1-\alpha$, $b=\alpha$, we get to the repeatedly mentioned exponentially weighted moving average 2, 3 4, with $\alpha=1/W$ we end up with 1 (except for the complications added to the algorithm).

Steady state performance For a step signal we expect $\mu_{k} = x_k$, let's assume $x_k = C$. Then we have $C = a C + b C$, i.e. $a = 1-b$, so we can restate the equation as $\mu_{k+1} = (\mu_k - x_{k+1}) a$, so the parameterization

Immediate response: Suppose an impulsive $x$, where $x_{k}=C$, and $x_j=0$ for $j \ne k$, then the average $\mu_{k+1}=C/W$, this requires $\alpha = 1/W$.

At this point I was satisfied that I was able to defend that answer as the optimal solution in some aspects. But then I noted that he was doing some complicated stuff there, combining some subsampling and exponential averaging, and actually did not satisfy the "no-buffer" requirement.