devol.dev

Filtering a noisy sensor

A sensor does not report the quantity you asked for. It reports that quantity plus everything else it is sensitive to: thermal drift in the bias, a tone coupled in from the mains, and the electrical noise floor of the part itself.

Those three corruptions look alike in a plot of the raw stream and are nothing alike in frequency. Drift lives at zero. Hum is a single line at 60 Hz. Noise is spread flat across the whole band. Filtering is choosing which parts of the frequency axis to keep, so the fix follows from where the problem lives, and a filter aimed at the wrong one does nothing but add delay.

Sensing / signal conditioning

Four filters on one noisy sensor

A 1 Hz reading arrives buried under white noise, mains hum, and a drifting bias. Pick a filter and watch both lanes: the stream on top, and what the filter does to every frequency underneath. Drag along the lower lane to move the corner.

RMS error
White noise
Lag at 1 Hz
Hum at 60 Hz

One equation, four filters

Samples arrive at a fixed rate fsf_s, so every filter here is the same kind of object: a weighted sum of recent inputs and recent outputs.

y[n]=kbkx[nk]    k1aky[nk]y[n] = \sum_{k} b_k\,x[n-k] \;-\; \sum_{k \ge 1} a_k\,y[n-k]

Write z1z^{-1} for one sample of delay and that becomes a ratio of polynomials, the transfer function

H(z)=b0+b1z1+b2z2+1+a1z1+a2z2+.H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots} {1 + a_1 z^{-1} + a_2 z^{-2} + \cdots}.

Feed in a sine wave at frequency ff and a sine wave at ff comes out, changed in two ways only: scaled, and shifted in time. Both are read off HH evaluated at z=ejωz = e^{j\omega}, where ω=2πf/fs\omega = 2\pi f / f_s is radians per sample. The magnitude H|H| is the scaling, plotted in the lower lane. The slope of the phase is the delay,

τ(ω)=dargHdω,\tau(\omega) = -\frac{d\,\arg H}{d\omega},

which is the number in the readout. There is no fifth thing a linear filter can do to a signal.

The four below differ only in where they put the roots of those two polynomials. Roots of the numerator are zeros, and a zero on the unit circle kills its frequency outright. Roots of the denominator are poles, and a pole near the circle lifts and prolongs its frequency.

Low pass

The exponential moving average, and the filter most likely to already be in your code:

y[n]=y[n1]+α(x[n]y[n1])y[n] = y[n-1] + \alpha\,\bigl(x[n] - y[n-1]\bigr)

One multiply, one state variable, no buffer. It is a single pole at r=1αr = 1 - \alpha, and placing that pole at

r=e2πfc/fsr = e^{-2\pi f_c / f_s}

puts the half-power point at fcf_c for any cutoff well below Nyquist. Gain at DC is exactly 1, so a constant reading passes through untouched.

Drag the cutoff down and the trace gets smoother and later. That is the whole trade, and the readout prices it: noise falls, lag rises.

One pole rolls off at 20 dB per decade, which is gentle, and it is why a low pass is a poor answer to hum specifically. Twenty dB down at 60 Hz is nearly free here. The next twenty costs a factor of ten in bandwidth, which puts the corner below the reading, and the one after that is not available at any price. Depth against a single tone is exactly what this filter cannot sell cheaply.

Moving average

The unweighted average of the last NN samples. Same job, different mechanism, and the differences are worth knowing.

y[n]=1Nk=0N1x[nk]H(ejω)=1Nsin(Nω/2)sin(ω/2)  ejω(N1)/2y[n] = \frac{1}{N}\sum_{k=0}^{N-1} x[n-k] \qquad\Longrightarrow\qquad H(e^{j\omega}) = \frac{1}{N}\, \frac{\sin(N\omega/2)}{\sin(\omega/2)}\;e^{-j\omega (N-1)/2}

Three consequences, all visible in the tool.

The phase term is exactly linear. Every frequency is delayed by the same (N1)/2(N-1)/2 samples. Nothing is distorted in shape, only moved, which is why this filter is popular for anything a human will look at.

There are exact nulls at every multiple of fs/Nf_s/N, where the window spans a whole number of periods and the sine cancels itself. Sampling at 200 Hz with N=10N = 10 puts nulls at 20, 40, 60 and 80 Hz, and the one at 60 removes mains hum completely. Choosing the window to span a whole number of hum cycles is a free notch, and it is the reason so many data acquisition systems average over exactly one mains period.

Between the nulls it leaks. The first sidelobe of a flat window is only about 13 dB down, so a moving average is a poor low pass anywhere it is not sitting on a null. Compare its curve to the one-pole above and the ripple is obvious.

Noise falls by 1/N1/\sqrt{N}: four samples halve it, and a hundred samples divide it by ten. Diminishing returns paid for in linearly growing delay.

High pass

The same pole as the low pass, plus a zero at z=1z = 1:

y[n]=r(y[n1]+x[n]x[n1])y[n] = r\,\bigl(y[n-1] + x[n] - x[n-1]\bigr)

That zero sits exactly on DC, so the gain there is exactly zero. Whatever standing bias the sensor carries is gone within a few time constants of startup, and the slow wander goes with it. Turn drift on and switch between low pass and high pass to see them separate.

Two things to note. The corner of a one-pole high pass lands somewhat above the fcf_c that placed the pole, because the zero at DC lifts the skirt on the way up, unlike the low pass built from the same pole. And a high pass does nothing about noise or hum, which are above it: it is a fix for one problem, the slow one.

This is what an accelerometer bias estimator, an AC-coupled input, and the complementary filter in an attitude estimator all have at their core.

Notch

Two zeros on the unit circle at the offending frequency ω0=2πf0/fs\omega_0 = 2\pi f_0/f_s, and two poles just inside it at the same angle:

H(z)=g12cosω0z1+z212rcosω0z1+r2z2,r=eπBW/fsH(z) = g\, \frac{1 - 2\cos\omega_0\,z^{-1} + z^{-2}} {1 - 2r\cos\omega_0\,z^{-1} + r^{2}z^{-2}}, \qquad r = e^{-\pi\,\mathrm{BW}/f_s}

The zeros make the gain at f0f_0 identically zero. The poles, at radius rr just inside them, pull the response back to unity a short way either side, and how close they sit is what sets the width. The constant gg normalises DC to 1.

Narrow is not free. A pole at radius rr decays with a time constant of 1/lnr-1/\ln r samples, or about

τ1πBW\tau \approx \frac{1}{\pi\,\mathrm{BW}}

seconds, so a 1 Hz notch rings for about a third of a second after any disturbance. Wide settles fast and takes more of the signal with it.

A notch has to be aimed. Move it off 60 Hz by more than its own width and the readout goes flat: it is doing nothing. That is the tradeoff of a filter this selective, and it is why a mains notch is usually paired with the knowledge of which country the robot is in.

Which one

Match the filter to where the corruption lives, not to how bad the trace looks.

ProblemWhere it livesFilter
Electrical noiseeverywherelow pass, or average
Thermal drift, biasat DChigh pass
Mains hum, a known toneone linenotch, or an average with a null on it
Aliased vibrationwherever it folded tonone of these, fix the sampling

The last row matters. Nothing after the sampler can undo aliasing, because a tone above Nyquist arrives already folded down to a frequency it does not actually have. That is an analogue filter’s job, before the converter.

Everything else on this page costs delay, and delay inside a control loop is phase margin spent. A filter that makes the trace look calm and the PID loop oscillate has made things worse. Filter as little as the estimator can live with.

The code

Four filters, one loop, and the difference between them is a table of numbers:

constexpr double kPi = 3.14159265358979323846;
constexpr double kFs = 200.0;   // sample rate, Hz

// Every filter here is this object. Two coefficient lists, a[0] = 1:
//
//   y[n] = sum_k b[k] x[n-k] - sum_{k>=1} a[k] y[n-k]
//
// Nothing downstream asks which filter a pair came from.
struct Coef {
    std::vector<double> b;
    std::vector<double> a;
};

// Pole radius for a corner at fc, clamped just short of the unit
// circle. At fc = 0 the exponential is exactly 1, and a pole sitting on
// the circle is an integrator that never forgets.
double poleRadius(double fc, double fs) {
    return std::min(0.99999, std::exp(-2 * kPi * fc / fs));
}

// One pole, the exponential moving average:
//   y[n] = y[n-1] + alpha (x[n] - y[n-1]),   alpha = 1 - r
Coef lowPass(double fc, double fs = kFs) {
    const double r = poleRadius(fc, fs);
    return {{1 - r}, {1, -r}};
}

// The same pole, plus a zero at z = 1:
//   y[n] = r (y[n-1] + x[n] - x[n-1])
// The zero is on DC, so the gain there is exactly zero.
Coef highPass(double fc, double fs = kFs) {
    const double r = poleRadius(fc, fs);
    return {{r, -r}, {1, -r}};
}

// n taps of 1/n, and no feedback at all: the impulse response is n
// samples long and then nothing.
Coef movingAverage(int n) {
    const int len = std::max(1, n);
    return {std::vector<double>(len, 1.0 / len), {1}};
}

// Two zeros on the unit circle at f0, two poles just inside it at the
// same angle. bw sets how far inside, and closer means narrower and
// ringing for longer.
Coef notch(double f0, double bw, double fs = kFs) {
    const double w0 = 2 * kPi * f0 / fs;
    const double r = std::min(0.9999, std::exp(-kPi * bw / fs));
    const double c = std::cos(w0);

    // Normalised so DC comes through at unity rather than at whatever
    // the pole pair happens to leave there.
    const double g = (1 - 2 * r * c + r * r) / (2 - 2 * c);

    return {{g, -2 * c * g, g}, {1, -2 * r * c, r * r}};
}

// Direct form I, zero initial state: the k <= n guards read history
// before the first sample as zero. That makes the first outputs a
// start-up transient rather than the filter's steady behaviour, so run
// a second of pre-roll through and throw it away.
std::vector<double> filter(const Coef& f, const std::vector<double>& x) {
    std::vector<double> y(x.size());

    for (std::size_t n = 0; n < x.size(); n++) {
        double acc = 0;
        for (std::size_t k = 0; k < f.b.size() && k <= n; k++)
            acc += f.b[k] * x[n - k];
        for (std::size_t k = 1; k < f.a.size() && k <= n; k++)
            acc -= f.a[k] * y[n - k];
        y[n] = acc / f.a[0];
    }

    return y;
}

// The moving average is the one filter that does not want to be run
// that way. As coefficients it is n taps, so the loop above spends n
// multiplies a sample and gets slower as the window grows. The
// difference equation cannot see that consecutive windows share n-1 of
// their terms. A ring buffer can, and costs one add and one subtract
// whatever n is.
//
// Same output as filter(movingAverage(n), x), same nulls at every
// multiple of fs/n. It sits outside the frame because it is an
// identity about the coefficients, not another set of them.
class RunningMean {
public:
    explicit RunningMean(int n) : buf_(std::max(1, n), 0.0) {}

    double operator()(double x) {
        sum_ += x - buf_[i_];
        buf_[i_] = x;
        i_ = (i_ + 1) % buf_.size();
        return sum_ / buf_.size();
    }

private:
    std::vector<double> buf_;
    double sum_ = 0;
    std::size_t i_ = 0;
};

The running sum carries its rounding error forward where the tap-by-tap version recomputes from scratch each sample, so the two drift apart in the last bits over a long stream. In double precision that is a few parts in 101610^{16} after a million samples, standing bias included, which is well under anything a sensor knows.

Explore further

  • Set the moving average to 10 samples and watch the hum vanish completely with no notch at all. Then move to 11 and watch total rejection become 21 dB.
  • Take a low pass down until 60 Hz is 40 dB down, and read what it cost the 1 Hz signal.
  • Turn drift on. Nothing but the high pass touches it, and the high pass has to go below 1 Hz before it is removing drift rather than signal.
  • Aim the notch at 45 Hz with the hum still at 60, and watch a filter do nothing at all.
  • Find the cutoff that minimises RMS error. It is not the lowest one: past a point, lag costs more than noise does.