First-Order Smoothing Filter
If you are dealing with analog measurements in your embedded applications, you might have also had to deal with filters. The analog anti-aliasing one is of course obligatory, however sometimes it is useful to filter the obtained measurements digitally too. This can either be to remove natural noise in the measurements, or even to trade between response time of the measured values and precision of the measurement.
This is really easy, right? Well I thought so, until I spent half a day debugging a bug in my filter that was actually a feature, I had completely blanked on. Therefore I decided to write down all the basics on this type of filter, if for nothing else, at least for my future reference.
Discrete-time first order smoothing filters considered in this posts are of the form below, where y(k) is the output of the filter at instant k, u(k) is the input of the filter at time k. Coefficients a and b specify the properties of the filter, like gain and cut-off frequency.
y(k) = b u(k) + a y(k-1)
The z-transform of the filter has the following form:
Y(z) = b U(z) + a Y(z) z^{-1}
From the above, the transfer function can be derived. This function has the following form:
F(z) = \frac{Y(z)}{U(z)} = \frac{b}{1 - a z^{-1}}
Generally, we want a filter to have gain of 1, in which case the filter can be characterized by single parameter known as forget factor \gamma. The two parameters of the filter can then be characterized using this factor as follows.
\begin{aligned} a &= 1 - \gamma \\ b &= \gamma \end{aligned}
The forget factor can be calculated, for given cut-off frequency f_c and sampling frequency f_s, as shown below.
\gamma = 1 - e^{-2 \pi f_c / f_s}
The filter can be on embedded devices implemented using fixed-point arithmetic
as shown in the example code below. The proposed filter's update() method
needs to be called with new input value with frequency given by sampling
frequency f_s.
#include <cstdint>
class Filter
{
public:
static const int FRACTION_BITS = 16;
struct Parameters
{
std::int32_t input_gain;
std::int32_t feedback_gain;
};
Filter():
previous_output_(0)
{ }
std::int16_t update(const Parameters & parameters, std::int16_t value)
{
std::int64_t accumulator;
//
// y(k) = B * u(k) + A * y(k-1)
//
// B * u(k)
accumulator = static_cast<std::int64_t>(parameters.input_gain) *
(static_cast<std::int64_t>(value) * (1ll << FRACTION_BITS));
// B * u(k) + A * y(k-1)
accumulator += static_cast<std::int64_t>(parameters.feedback_gain) *
static_cast<std::int64_t>(previous_output_);
// k <- k + 1
previous_output_ = static_cast<std::int32_t>(
accumulator / (1ll << FRACTION_BITS));
return static_cast<std::int16_t>(
previous_output_ / (1l << FRACTION_BITS));
}
private:
std::int32_t previous_output_;
};To verify the proper behavior of the smoothing filter, it is possible to measure the time constant \tau of the filter's step response. The expected time constant can be calculated using equation below.
\tau = \frac{1}{2 \pi f_c}
Update
As it turned out this topic is more interesting than I originally thought. As pointed out by my colleague, the implementation above can be greatly improved, if we constrain it slightly. What is wrong with the above implementation? Firstly, it uses often non-native 64-bit arithmetic, and secondly, the output of such filter will never actually reach its input value.
The latter is what we would expect in ideal contiguous filter - the output of the filter asymptotically approaching the input value (assuming the input value is constant), but never actually reaching it. This is, however, not the case in the discrete space. Thanks to the quantization error the output value shall reach the constant input value in finite time.
First constraint of the improved filter implementation is that it can only implement the filters of the following time-domain form:
y(k) = \gamma u(k) + (1 - \gamma) y(k-1)
This, however, is not at all limiting, as this is the form of the filter we wanted to implement in the first place - a filter with gain value of 1. This filter can also be represented with the following diagram:
+----------------+
| |
+-------+ + v + |
u(k) --->| gamma |--->(_)------------------------------+-> y(k)
+-------+ ^ - | |
| | |
| +-------+ | +------+ |
+--| gamma |<----+----| z^-1 |<-+
+-------+ y(k-1) +------+The equation above can be further transformed into more simplified form, which removes one multiplication operation:
y(k) = \Big( u(k) - y(k-1) \Big) \gamma + y(k-1)
The second constraint deals with the representation of the filtered signal value. In the implementation below the filtered signal value is represented as an unsigned 16-bit integer. This again is not significant limitation, as most likely your ADC will generate unsigned number of less than 16 bit length (typically 8, 10, or 12 bits).
#include <cstdint>
class Filter
{
public:
static const int FRACTION_BITS = 16;
struct Parameters
{
/**
* @brief Filter parameter represented as fixed-point number
*
* @see @ref Filter::FRACTION_BITS
*/
std::uint32_t forget_factor;
};
FirstOrderFilter():
previous_output_(0)
{ }
std::uint16_t output() const
{
return previous_output_ >> FRACTION_BITS;
}
void step(std::uint16_t value, const Parameters & parameters)
{
std::uint32_t output;
output = value;
output -= output();
output *= parameters.forget_factor;
output += previous_output_;
previous_output_ = output;
}
private:
std::uint32_t previous_output_;
};Further benefit of the improved implementation is the fact that it does not utilize division operation (only division-like operation is bit-shift). This is significant for platforms which do not have division instruction (namely ARM Cortex-M0 and Cortex-M0+). Furthermore the multiplication can also be done away with by replacing it with bit-shift. That is only possible if we further constraint forget factor \gamma to be (negative) power of two, e.g. 1, 0.5, 0.25, 0.125, etc.
I would like to give credit to my colleague Miroslav Stehlik for this update.