Suppose you need to filter your noisy measurement a bit. Digital filter, oh, and those Z transforms... Most engineers I know would roll out a moving average filter (beware of overflow and aliasing!), few would pull their Mitra and Proakis from the bookshelves, looking up Butterworth again. But there is a simpler way.
Ignoring computer arithmetic imperfections for now, these two lines are equivalent (thinking in C, x
is an input expression, y
output variable, say double
, float
or int
):
y = x;
y += x - y;
Now suppose we add just a fraction of a difference:
y += k * (x - y); // 0 < k < 1
This is a beautiful first order IIR (recursive) low-pass filter. y
is its state, x - y
innovation in each iteration.
Consider the step answer. The difference to 1 vanishes like 1, 1 - k, (1 - k)^2, ...
, so the time constant T
and cutoff frequency fc
can be calculated like this:
// (1 - k)^(T/dt) = exp(-1) =>
T = - dt / log(1 - k) =~ dt / k
fc = 1/(2*pi*T) = -log(1 - k) * fs / (2*pi) =~ k * fs / (2*pi)
// fs - sampling frequency, dt = 1/fs - sampling period
// log(1 + x) =~ x for small x
So it needs about 1 / k
samples to reach 63% of constant input.
In integer case:
y += (x - y) >> shift;
With rounding:
y += (x - y + (1 << (shift-1))) >> shift;
One line, no Z transforms, no memory. And if you need a high-pass filter, x - y
is your filter output.
The Matlab form is
y = filter(k, [1, k - 1], x);