Simple smoothing filters

A kitchen scale’s readings jitter by about 0.2 kg, and a 1 kg parcel lands on it. Below, the readings go through averages of 1, 4 and 16 points. Watch the line grow steadier while its jump comes later.

A longer average: less noise, more delay

Picture a kitchen scale. Its load cell is read 50 times a second, so a new reading arrives every 20 ms and $f_s=50$ Hz. Each reading is off by a little random noise, about $x_\text{noise,rms}=0.2$ kg, using the RMS of How big is a signal (1.3). At $t=2$ s a 1 kg parcel lands on the scale.

Nobody wants a display that jumps around by 0.2 kg. The first fix most people try is the moving average of Difference equations (6.1): show the mean of the last $N$ readings.

Why does averaging help? The noise in one reading has nothing to do with the noise in the next. For unrelated noises like these, the powers add: the power of a sum is the sum of the powers. Power here is the mean square of 1.3, the square of the RMS.

So the sum of $N$ readings carries $N$ times the noise power of one reading. Dividing the sum by $N$ divides its power by $N^2$, which leaves $1/N$ of the power of one reading. The RMS is the square root of the power, so the average keeps

of the noise. Four readings halve it, and sixteen leave a quarter.

The same argument works for any filter that adds up weighted readings, $y[n]=\sum_k h[k]\,{x[n-k]}$. Each term carries $h[k]^2$ times the noise power of one reading, and the terms add. So the noise that comes out has RMS $x_\text{noise,rms}\sqrt{\sum_n h[n]^2}$.

I call $\sqrt{\sum_n h[n]^2}$ the filter’s noise gain. For the average every tap is $1/N$, so $\sum_n h[n]^2=N/N^2=1/N$, and the noise gain is $1/\sqrt N$ as before.

Averaging has a price. The average’s $h$ is $N$ equal taps, symmetric about the centre $(N-1)/2$. By Linear-phase systems (17.3), symmetric taps delay every frequency by that many samples. So the average reports each change $(N-1)/2$ samples late, and it spreads a sudden jump over $N$ samples.

The picture at the top of the page runs the scale’s readings through averages of 1, 4 and 16 points.

Watch the noise readout fall while the delay readout grows. Sixteen points keep a quarter of the noise, and the jump is 150 ms late. It is like a phone poll averaged over a week: it is steadier than one day’s poll, but it tells you about the past week, not about today.

The noise left in the picture comes from the formula. The noise on this particular trace measures 0.197, 0.098 and 0.050 kg, close to it.

You can see the delay in the numbers too. The true weight jumps between $n=99$ and $n=100$, at $n=99.5$. The 4-point average reads 0.25, 0.5, 0.75 and 1 kg at $n=100$ to $103$, so it passes the halfway mark at $n=101$, 1.5 samples late.

After the clip, set $N$ yourself with the slider, from 1 to 32.

A 9-point average and its one-multiply twin

A 16-point average has to keep the last 16 readings in memory. Sensor code often uses the leaky integrator of 6.1 instead, $y[n]=a\,{y[n-1]}+(1-a)\,x[n]$ with $0<a<1$. In sensor code and in statistics it is called the exponential moving average, or EMA.

It stores one number, the last output. Written as $y[n]={y[n-1]}+(1-a)\bigl(x[n]-{y[n-1]}\bigr)$, it costs one multiply per reading.

Its impulse response is $h[n]=(1-a)a^n$ for $n\ge0$, and it never ends. So there is no centre tap to read a delay from.

Instead I use the balance point of $h$, the centre $\sum_n n\,y[n]\big/\sum_n y[n]$ of 17.3. The taps of both filters add up to 1, so the balance point is $\sum_n n\,h[n]$. I call it the mean delay. For the $N$-point average it is $(N-1)/2$ again.

For the EMA, let $S=a+2a^2+3a^3+\dots$, the sum of $n\,a^n$. Subtracting $aS$ leaves $a+a^2+a^3+\dots$, which is $a/(1-a)$ by 6.1’s geometric sum, because $a^{n+1}$ shrinks to 0. So $S=a/(1-a)^2$, and

The noise gain comes from the same kind of sum. Each $h[n]^2$ is $(1-a)^2a^{2n}$, and the powers $a^{2n}$ add up to $1/(1-a^2)$. Since $1-a^2=(1-a)(1+a)$,

Now ask for an EMA that matches an $N$-point average. Equal noise needs $(1-a)/(1+a)=1/N$. Equal mean delay needs $a/(1-a)=(N-1)/2$. Solve either one and you get the same answer:

I call this the twin rule. For $N=9$ it gives $a=0.8$.

Think of guessing when your bus will come. You could keep a notebook of the last nine trips and average them. Or you could nudge your guess a little after each trip, which needs no notebook. Below, the two work side by side on the scale’s readings.

Notice that the four readouts come in two equal pairs. The shapes of the jumps differ, though.

The average’s ramp is straight and finishes. The EMA’s step response is $1-0.8^{n+1}$ of 6.1: it is 0.2, 0.36 and 0.488 of the way after one, two and three samples, and it never quite arrives.

After the clip, set the EMA’s $a$ yourself with its slider. Try 0.9: that EMA is the twin of a 19-point average.

The two twins remove the same amount of noise, but not from the same frequencies. The average’s gain is zero at $\Omega=2\pi k/N$ (Frequency response of discrete-time systems, 12.4). In hertz, by Frequency in discrete time (12.1), these nulls sit at

For 9 points at 50 Hz they are the multiples of 5.56 Hz. The EMA’s gain follows from 12.4’s method, as for the leaky integrator there:

The denominator is never zero, so the EMA has no nulls.

Choosing a

There are two common ways to pick $a$. The first starts from a time constant. In Differential equations and analog systems (6.2) the leaky integrator is the step-by-step copy of an RC circuit, with $a=e^{-T_s/RC}$. Here $RC$ is the time constant $\tau$, so

At 50 Hz, a time constant of 0.2 s gives $a=e^{-0.02/0.2}=e^{-0.1}=0.905$. That EMA leaves 0.045 kg of noise and has a mean delay of 9.5 samples, 190 ms. It is the twin of a 20-point average.

The second starts from an equivalent $N$ and uses the twin rule. The pandas library’s ewm(span=N) sets $1-a=2/(N+1)$, which is the same rule. The internet’s TCP protocol tracks the round-trip time of its packets with an EMA of $a=7/8$ (RFC 6298). That is the twin of a 15-point average.

The average’s nulls are useful by themselves. Mains power puts a hum at 50 Hz, and often at 100 Hz, 150 Hz and so on. Sample at 1 kHz and average 20 points, exactly one 20 ms mains period. The nulls fall at every multiple of 50 Hz, so the hum and all its harmonics are removed.

In a country with 60 Hz mains, the same average keeps 0.157 of the hum. Averaging 100 ms instead, 100 points at 1 kHz, puts a null at every multiple of 10 Hz, which covers both 50 and 60 Hz. That is why bench multimeters measure over a whole number of mains cycles.

A median ignores spikes

Not every problem is steady noise. Sometimes a single reading is far off, because a wire made a bad contact or someone bumped the scale. I call such a reading a glitch. An average spreads a glitch into a bump that lasts $N$ samples.

There is a smoother that ignores glitches. Sort the last 5 readings and take the middle one, their median. For 0, 0, 2, 0, 0 the sorted list is 0, 0, 0, 0, 2, and the middle value is 0. A median filter outputs the median of the last $N$ readings at every sample.

With 5 readings the middle value is the third. One glitch, or two, sort to an end of the list and never reach the middle. A jury’s middle score works the same way: one judge’s wild mark does not move it. Below, a 5-point average and a 5-point median work on the same readings.

Notice that the two glitches side by side give the average its biggest bump, 0.8. Each glitch adds $2/5=0.4$ to every average that contains it, and two together add 0.8. The median does not move at all, and its step is as sharp as the input’s.

The median pays for this by giving up linearity. Apply the scale-and-add test of System properties (4.2, “Does scaling and adding pass through?”) to a 3-point median.

Take two windows of readings, $\{1,0,0\}$ and $\{0,0,1\}$. Each has median 0, so the medians add to 0. But the sum of the windows is $\{1,0,1\}$, whose median is 1.

So a median filter is not linear, and it has no impulse response that describes it. Feed it a single 1 and it outputs only zeros, yet its output for a step is a step. Noise gain, nulls and frequency response do not apply to it. On a flat stretch, a median of $2K+1$ readings removes bursts of up to $K$ glitches in a row; a burst of $K+1$ outweighs the rest and survives.

Three smoothers, one trace

Here are all three on the scale’s readings from the first instrument, with four glitches added. Three are 1.5 kg up, at samples 40, 41 and 150, and one is 1.5 kg down, at sample 170. Each smoother has 9 points, or the 9-point twin, so all three lag by about 4 samples.

Look at the glitches at 0.8 s. The average turns the pair into a flat bump of 0.33 kg that lasts 9 samples. The EMA jumps by 0.54 kg and then creeps back.

The median barely moves, and its step at 2 s is the sharpest of the three. But it keeps the most noise, 0.42 of it, because for noise like this the middle value of 9 readings is a rougher estimate than their mean.

So which one should you choose? If the problem is noise and the delay does not matter, use an average, or an EMA when memory or computing time is short, as on a small microcontroller. If the problem is a fixed hum, average over exactly its period. If the problem is glitches or dropouts, use a median, often followed by an average to take out the noise that is left.

Worked example

Let’s work out the scale’s 9-point average and its twin by hand, at $f_s=50$ Hz with $x_\text{noise,rms}=0.2$ kg.

Noise and delay of the average. The noise left is $0.2/\sqrt9=0.2/3=0.067$ kg. The delay is $(9-1)/2=4$ samples, which is $4\times20=80$ ms.

Nulls of the average. They sit at multiples of $50/9$ Hz: 5.56, 11.11, 16.67 and 22.22 Hz, all below the top frequency of 25 Hz.

At 25 Hz itself, $\Omega=\pi$ and $e^{-j\pi k}=(-1)^k$. The nine taps alternate in sign, five plus and four minus, so the gain is $1/9=0.111$. Solving for the gain $1/\sqrt2$ numerically puts the −3 dB point at 2.47 Hz.

The twin. The twin rule gives $a=8/10=0.8$. Its noise gain is $\sqrt{0.2/1.8}=1/3$, so it also leaves 0.067 kg, and its mean delay is $0.8/0.2=4$ samples. At 25 Hz, $\cos\Omega=-1$, so the EMA’s squared gain is $(1-a)^2/(1+a)^2$. Its gain is $(1-a)/(1+a)=0.2/1.8=0.111$, the same as the average’s.

The EMA’s −3 dB point. Set the squared gain to $\tfrac12$: $2(1-a)^2=1-2a\cos\Omega+a^2$. Rearranged, this gives

So $\Omega=\cos^{-1}0.975=0.224$ rad/sample, and $f=0.224\times50/2\pi=1.78$ Hz. The EMA’s gain falls faster at first than the average’s, but it never reaches zero.

From a time constant, and TCP’s twin. For $\tau=0.2$ s, $a=e^{-0.1}=0.9048$, the noise gain is $\sqrt{0.0952/1.9048}=0.2235$ (0.045 kg), and the mean delay is $0.9048/0.0952=9.51$ samples (190 ms). Its twin has $N=1.9048/0.0952=20.0$ points. For TCP’s $a=7/8$, the twin has $N=15$ and the noise gain is $\sqrt{1/15}=0.258$.

Median against average. On the closing figure’s readings, the three smoothers compare like this. The raw readings have 0.203 kg of noise on the flat stretches.

The average and the EMA keep about a third of the noise, as the formula $1/\sqrt9$ predicts. The median keeps more noise and moves least at the glitches.

Where you’ll meet this

Sensor firmware is full of these three. A temperature or weight reading on a small microcontroller is often an EMA, because it needs one stored number. A display that must not flicker uses a short average. A reading that sometimes glitches gets a median first.

The internet’s TCP protocol estimates round-trip time with an EMA of $a=7/8$, as above. Rolling means of daily counts, such as a 7-day average of cases or sales, are moving averages whose nulls remove the weekly pattern. Median filters clean isolated bad pixels from images, which is the subject of Image filtering (28.2).

Smoothers with unequal weights, such as the Savitzky–Golay filter, which fits a polynomial to each window, come up with Special FIR filters (19.4). Moving averages built in hardware, called CIC filters, are in Polyphase structures (22.4). The Kalman filter, whose steady state is an EMA, is in RLS and the Kalman filter (26.4). A full sensor chain is in Sensor signal conditioning (32.2).

Reference card