5
\$\begingroup\$

Exponential damping is a type of "moving average", but it's not an arithmetic mean. The latter has the disadvantage that it requires storage of the last N samples in a std::deque or similar.

Exponential damping has constant, and very small, storage costs of just 3 or 4 primitive types, no matter how "slow moving" or "aggressively smoothing" the damping is. The algorithm is also trivially simple, so this type of "moving" averaging is suitable for use in microcontrollers or very high speed systems. In fact, some microcontrollers such as Microchip's PIC series implement this algorithm in hardware, to smooth noisy analog input signals. (Their hardware version doesn't provide the "startup bridge" and the time_constant is in integer powers of 2, saving the "expensive" integer division. See below).

This code review proposes a very simple generic class to make exponential damping available with flexible types to tailor it to the platform and data types. I have called this class damper. Its only constructor takes a single parameter which I have called time_constant for reasons, which are explained below.

Here is damper smoothing the FramesPerSecond readout of a real time physics simulation (Note: these are not "visual frames" sent to the video card, but merely "frames" in the simulation").

Frames per second damped with damper(300)

That is probably "over smoothing", so lets drop the time_constant to 30.

Frames per second damped with damper(30)

In order to gain some intuition of what this algorithm does, we can feed it with a "step function". We can see that the response is very similar to an exponential decay function, and is in fact very closely related, mathematically. We use this property below in our unit tests. This similarity is also the origin of the name of the constructor parameter: time_constant. The higher this value, the more aggressive the smoothing.

Testing damper(10) with a step function

The algorithm is so simple, that it's all done in 6 or 7 lines of code, and it really only has a single useful member function, for which we have used operator(). We have included an internal count_ to bridge the startup phase. This avoids damper only slowly ramping up to a steady state set of input samples.

To make it slightly more challenging I have tried to make the class completely flexible on types, and this did indeed require some interesting testing regards "signedness", "narrowing", and warnings of "loss of precision". So I have turned the compiler warnings way up and included a set of unit tests below. We are using:

clang++-13 -std=c++20 -O3 -Wall -Wextra -Wpedantic -Wconversion -Wshadow -Wextra-semi -Werror

There is also fairly heavy commenting of what I found to be the tricky areas in dealing with the tested range of types. I chose to support and test the full range of c++ primitive integral and fp types, but nothing more; ie not user defined types, to keep it simple, fast and small memory footprint.

I attempted to algebraically remove the need for an internal sum_ member, because this is in danger of overflowing, if the Sum type is inappropriately chosen by the user - she must choose Sum to be large enough to hold average_sample * time_constant. I was unable to find a way for the class implementation to protect the user from this danger. There is a commented out test illustrating the problem.

It would be possible to omit the damped_value_ member but I decided for clarity and speed that it was better in. reset() and current() are trivial, and probably will only see rare usage.

The final TEST shows a minimal damper, which occupies just 4 bytes of storage, for use in a uC, perhaps.

damper.hpp

#pragma once

#include <concepts>
#include <type_traits>

// Dampens a "noisy" value using "approx exponential damping".
// Submit samples using `operator()`, which returns the damped value.
// MUST choose a SUM type which can hold: sample avg * time_constant
template <typename Value, typename Sum = Value,
          // type switch is required to prevent warnings AND incorrect arithmetic with -ves
          typename Count = std::conditional_t<std::is_signed_v<Value>, short, unsigned short>>
requires std::floating_point<Value> || std::integral<Value>
class damper {
  public:
    explicit damper(Count time_constant) : time_constant_(time_constant) {}

    Value operator()(Value sample) {
        if (count_ != time_constant_) {
            // branch avoids distortions for the first time_constant samples
            sum_ += sample;
            ++count_;
        } else {
            // Main running branch: once system is "primed" with time_constant number of samples
            // provides "approximately exponential damping" with the given time constant

            // correct, and well defined, even with unsigned types
            // but there remains is potential for overflow, if user choses a Sum type which cannot
            // hold sample avg * time_constant
            sum_ += sample - damped_value_;
        }
        // first static_cast suppresses -Wimplicit-int-conversion about very small types being
        // implicitly promoted to integer then demoted back again to be assigned to a small Value
        // type.

        // second static_cast is about a -Wimplicit-int-float-conversion of count_ from
        // eg (unsigned) int to float which is never likely to be relevant
        damped_value_ = static_cast<Value>(sum_ / static_cast<Sum>(count_));

        return damped_value_;
    }

    Value current() { return damped_value_; }

    void reset() {
        sum_          = Sum{};
        count_        = Count{};
        damped_value_ = Value{};
    }

  private:
    Sum         sum_{};
    Value       damped_value_{};
    Count       count_{};
    const Count time_constant_;
};

damper_tests.cpp

#include "damper.hpp"
#include "gtest/gtest.h"
#include <cmath>
#include <concepts>
#include <limits>
#include <type_traits>

// for testing purposes only
template <typename Value, typename Sum = Value,
          // type switch is required to prevent warnings AND incorrect arithmetic with -ves
          typename Count = std::conditional_t<std::is_signed_v<Value>, int, unsigned>>
requires std::floating_point<Value> || std::integral<Value>
void test_damper(Value pre_step = 100, Value post_step = 0, Count tc = 10) {

    // modelling a step response from Value{pre_step} down to Value{post_step} after tc samples
    // expecting flat Value{100} and then "approximately exponential decay" towards Value{0}
    // `tc` is aprox equiv to the "time constant" in exponential damping

    // test empty damper
    auto d = damper<Value, Sum, Count>(tc);
    static_assert(std::is_same_v<decltype(d.current()), Value>);
    EXPECT_EQ(d.current(), Value{});

    // tc samples with pre_step value
    for (auto i = Count{1}; i <= tc; ++i) {
        auto dv = d(pre_step);
        static_assert(std::is_same_v<decltype(dv), Value>);
        EXPECT_EQ(dv, pre_step); // check at every step that the `damped_value_` is "flat"
    }
    // `damper` is now fully "primed" with pre_step values

    // tc further samples of `post_step` values
    for (auto i = Count{1}; i <= tc; ++i) {
        auto dv = d(post_step);
        static_assert(std::is_same_v<decltype(dv), Value>);

        // check at every sample that we are matching the exponential decay curve
        // which predicts "approximately exponential damping"
        auto expected = post_step + (pre_step - post_step) * std::pow((tc - 1.0L) / tc, i);

        if constexpr (std::is_integral_v<Value>) {
            auto predicted = static_cast<Value>(std::round(expected));
            // allow an integer result to be "off-by-one" after rounding
            // std::abs can't be used here, because programm would be ill-formed for some types
            EXPECT_TRUE(dv - predicted == Value{0} || dv - predicted == Value{1} ||
                        predicted - dv == Value{1});
        } else {
            auto predicted = static_cast<Value>(expected);
            // for FP we expect the prediction to be within appropiately scaled epsilon
            EXPECT_LT(std::fabs(dv - predicted),
                      std::fabs(pre_step - post_step) * std::numeric_limits<Value>::epsilon());
        }
    }
}

// signed integer types
TEST(damper, short) { test_damper<short>(); }
TEST(damper, int) { test_damper<int>(); }
TEST(damper, long) { test_damper<long>(); }
TEST(damper, long_long) { test_damper<long long>(); }

// signed integer types with a step from negative to positive
TEST(damper, short_negstep) { test_damper<short>(-100, 100); }
TEST(damper, int_negstep) { test_damper<int>(-100, 100); }
TEST(damper, long_negstep) { test_damper<long>(-100, 100); }
TEST(damper, long_long_negstep) { test_damper<long long>(-100, 100); }

// unsigned integer types
TEST(damper, unsigned_short) { test_damper<unsigned short>(); }
TEST(damper, unsigned) { test_damper<unsigned>(); }
TEST(damper, unsigned_long) { test_damper<unsigned long>(); }
TEST(damper, unsigned_long_long) { test_damper<unsigned long long>(); }

// FP types
TEST(damper, float) { test_damper<float>(); }
TEST(damper, double) { test_damper<double>(); }
TEST(damper, long_double) { test_damper<long double>(); }

// FP types with a step from negative to positive
TEST(damper, float_negstep) { test_damper<float>(-100.0, 100.0); }
TEST(damper, double_negstep) { test_damper<double>(-100.0, 100.0); }
TEST(damper, long_double_negstep) { test_damper<long double>(-100.0, 100.0); }

// tiny integer types (one use case is microcontroller analog input damping)
TEST(damper, uint8_t__uint16_t) { test_damper<std::uint8_t, std::uint16_t>(); }
TEST(damper, uint8_t__uint16_t__uint8_t) {
    test_damper<std::uint8_t, std::uint16_t, std::uint8_t>();
}

// this test fails, because the sum overflows
// I have tried to algebraically eliminate the sum, but it seems that, for integer arithmetic,
// we always need at least a temporary result which can hold sample avg * time_constant
// TEST(damper, uint8_t__uint8_t__uint8_t) { test_damper<std::uint8_t, std::uint8_t,
// std::uint8_t>(); }

// however, it is still a valid set of template params when sample avg * tc is small
TEST(damper, uint8_t__uint8_t__uint8_t_small_values) {
    test_damper<std::uint8_t, std::uint8_t, std::uint8_t>(0, 40, 5);
}

\$\endgroup\$

3 Answers 3

1
\$\begingroup\$

Focusing in on the type of Count, following comments on my other answer.

I see why you made it conditionally signed. I think that could be simplified. I was able to keep it unsigned, but still work with signed Sum type (also defining a concept so we can use abbreviated template constraints):

template<typename T>
concept arithmetic = std::is_arithmetic_v<T>;

// Dampens a "noisy" value using "approx exponential damping".
// Submit samples using `operator()`, which returns the damped value.
// MUST choose a SUM type which can hold: sample avg * time_constant
template <arithmetic Value,
          arithmetic Sum = Value,
          std::unsigned_integral Count = unsigned short>
class damper

We then check that that we don't use a time constant that's too large to convert to the sum type:

explicit damper(Count time_constant)
    : time_constant{time_constant}
{
    if (time_constant <= 0) {
        throw std::invalid_argument("damper needs positive time constant");
    }
    if constexpr (std::is_integral_v<Sum> && std::is_signed_v<Sum>) {
        if (!std::in_range<Sum>(time_constant)) {
            throw std::invalid_argument("time constant too large for signed arithmetic");
        }
    }
}

And the final calculation divides using the sum type:

    return damped_value = static_cast<Value>(sum / static_cast<Sum>(count));

This passes the full provided test suite, and this additional test:

TEST(Damper, MixedSignedness)
{
    test_damper<int, int, unsigned int>();
    EXPECT_THROW((damper<int, int, unsigned>(std::numeric_limits<unsigned>::max())),
                 std::invalid_argument);
}
\$\endgroup\$
6
  • 1
    \$\begingroup\$ Yes, these are good points and the types of aspects I have been considering. I will think about them more carefully a bit later. I extracted the signedness problem for an SO question: stackoverflow.com/questions/25609091/… and the majority of comments coming in seem to suggest "prevent the mixed signedness situation in the first place". \$\endgroup\$ Commented Feb 18, 2022 at 13:07
  • 1
    \$\begingroup\$ Yes, that would be the other alternative - but if so, then really prevent it (with a constraint) rather than just defaulting it. \$\endgroup\$ Commented Feb 18, 2022 at 18:59
  • \$\begingroup\$ Yes, i have some code playing with static_asserts or concept requires clauses or a mixture \$\endgroup\$ Commented Feb 18, 2022 at 19:04
  • 1
    \$\begingroup\$ I have accepted this answer rather than the other one, because I think the type juggling is the most challenging bit. I think I am ultimately going to go with some static_asserts and I might get away with just one single static_cast (for assigning back into damped_value_ because that is a genuine narrowing). gcc has bugs in this area I believe, see here: stackoverflow.com/questions/71181566/… Thanks for your input. \$\endgroup\$ Commented Feb 19, 2022 at 12:02
  • 1
    \$\begingroup\$ See here: codereview.stackexchange.com/a/274262/212940 \$\endgroup\$ Commented Feb 19, 2022 at 13:19
1
\$\begingroup\$

Pretty good code. And I appreciated the clear explanation!

We have a reset() function that's not tested. However, I don't think we need that, since simply assigning a newly-constructed damper could be made to have the same effect:

damper d;
d = damper(d.time_constant());  // same as d.reset()

We'd need to remove the const qualifier on the member object, and make it visible for this. Or just test the function.

We carefully select a default Count type of same signedness as Value (shouldn't it be the same as Sum instead?), but that's pointless since we cast it where we use it. I think we can just use unsigned short for that:

template <typename Value, typename Sum = Value,
          std::unsigned_integral Count = unsigned short>
    requires std::floating_point<Value> || std::integral<Value>

I got a warning from the += in the main running branch; that's easily fixed:

        sum_ += static_cast<Sum>(sample - damped_value_);

If we push a NaN value, that poisons the damper forever; we should reject or ignore them:

Value operator()(Value sample)
{
    if (std::isnan(sample)) {
        // ignore it
        return damped_value_;
    }

The current() member doesn't change object state:

Value current() const { return damped_value_; }

The long_double tests fail for me (GCC, Linux, amd64):

274171.cpp:109: Failure
Expected: (std::fabs(dv - predicted)) < (std::fabs(pre_step - post_step) * std::numeric_limits<Value>::epsilon()), actual: 1.42109e-14 vs 1.0842e-17
274171.cpp:109: Failure
Expected: (std::fabs(dv - predicted)) < (std::fabs(pre_step - post_step) * std::numeric_limits<Value>::epsilon()), actual: 1.42109e-14 vs 1.0842e-17
274171.cpp:109: Failure
Expected: (std::fabs(dv - predicted)) < (std::fabs(pre_step - post_step) * std::numeric_limits<Value>::epsilon()), actual: 2.13163e-14 vs 1.0842e-17
274171.cpp:109: Failure
Expected: (std::fabs(dv - predicted)) < (std::fabs(pre_step - post_step) * std::numeric_limits<Value>::epsilon()), actual: 2.13163e-14 vs 1.0842e-17

We probably met the limits of std::pow(), I think.

It might be a good idea to push more than tc initial values, to confirm that the damped signal remains flat as we cross the count boundary.

We don't need to repeatedly static_assert the type of d.current() during the test - types can't change within a template instantiation. One assertion at the beginning is enough.

Instead of having a separate test path for integer and floating-point dampers, we can use a single path, by always doing the comparison using long double:

auto const maxerr = allowed_error(pre_step, post_step);

// tc further samples of `post_step` values
for (auto i = Count{1}; i <= tc; ++i) {
    auto dv = d(post_step);
    // check at every sample that we are matching the exponential decay curve
    // which predicts "approximately exponential damping"
    auto expected = post_step + (pre_step - post_step) * std::pow((tc - 1.0L) / tc, i);

    EXPECT_LT(std::fabs(static_cast<long double>(dv) - expected), maxerr);
}

That's using a helper function which is overloaded to select the appropriate error limit:

template<std::integral T>
auto constexpr allowed_error(T, T)
{
    // integer dampers are allowed to be out by up to 1 unit
    return 1.0L;
}
template<std::floating_point T>
auto constexpr allowed_error(T from, T to)
{
    return std::fabs(to - from) * std::numeric_limits<T>::epsilon();
}

I thought we should be using EXPECT_NEAR here, but that seems not to be implemented for long double. That said, I'd be happy testing with double and leaving long double less tested.

We can and should use SCOPED_TRACE to identify which iteration(s) the assertions fail:

for (auto i = Count{1}; i <= tc; ++i) {
    SCOPED_TRACE(std::to_string(i));

When we write tests against a reference implementation, it's clearer to have the reference stand as a separate function, rather than being intertwined with the tests themselves:

auto reference_damping(double pre_step, double post_step,
                       double smoothing, unsigned iterations)
{
    auto constexpr f = (smoothing - 1) / smoothing;
    return post_step
        + (pre_step - post_step) * std::pow(f, iterations);
}

Finally, Google Test discourages test cases with underscore in the name (as that's what it uses to construct its class names).


Modified code

#include <cmath>
#include <concepts>
#include <stdexcept>
#include <type_traits>

// Dampens a "noisy" value using "approx exponential damping".
// Submit samples using `operator()`, which returns the damped value.
// MUST choose a SUM type which can hold: sample avg * time_constant
template <typename Value, typename Sum = Value,
          std::unsigned_integral Count = unsigned short>
    requires std::floating_point<Value> || std::integral<Value>
class damper
{
    Sum         sum{};
    Value       damped_value{};
    Count       count{};
    const Count time_constant;

  public:
    explicit damper(Count time_constant)
        : time_constant{time_constant}
    {
        if (time_constant <= 0) {
            throw std::invalid_argument("damper needs positive time constant");
        }
    }

    Value operator()(Value sample)
    {
        if (std::isnan(sample)) {
            // ignore it
            return damped_value;
        }
        if (count != time_constant) {
            // branch avoids distortions for the first time_constant samples
            sum += sample;
            ++count;
        } else {
            // Main running branch: once system is "primed" with
            // time_constant number of samples provides "approximately
            // exponential damping" with the given time constant

            // correct, and well defined, even with unsigned types but
            // there remains is potential for overflow, if user choses
            // a Sum type which cannot hold sample avg * time_constant
            sum += static_cast<Sum>(sample - damped_value);
        }
        // first static_cast suppresses -Wimplicit-int-conversion
        // about very small types being implicitly promoted to integer
        // then demoted back again to be assigned to a small Value
        // type.

        // second static_cast is about a -Wimplicit-int-float-conversion
        // of count from eg (unsigned) int to float which is never likely
        // to be relevant
        return damped_value = static_cast<Value>(sum / static_cast<Sum>(count));
    }

    Value current() const { return damped_value; }
};

Modified tests

#include "damper.hpp"
#include <gtest/gtest.h>
#include <cmath>
#include <concepts>
#include <limits>
#include <string>
#include <type_traits>

// Reference implementation to compare against: an exponential decay
// curve which predicts "approximately exponential damping"
auto reference_damping(double pre_step, double post_step,
                       double smoothing, unsigned iterations)
{
    auto constexpr f = (smoothing - 1) / smoothing;
    return post_step
        + (pre_step - post_step) * std::pow(f, iterations);
}

// Overload set to determine how close the result needs to match the
// reference implementation:
template<std::integral T>
constexpr double allowed_error(T, T)
{
    // integer dampers are allowed to be out by up to 1 unit
    return 1.0;
}
template<std::floating_point T>
constexpr double allowed_error(T from, T to)
{
    return std::fabs(to - from) * std::numeric_limits<T>::epsilon();
}
constexpr double allowed_error(long double from, long double to)
{
    // long double only needs to be as precise as double
    return allowed_error(static_cast<double>(from), static_cast<double>(to));
}


// The test function
template<typename Value, typename Sum = Value, typename Count = unsigned>
void test_damper(Value pre_step = 100, Value post_step = 0, Count tc = 10)
{
    // This test exercises the step response from pre_step to to
    // post_step.  It fills the damper before switching, and then
    // expects approximately exponential decay towards post_step.
    // `tc` approximates to the "time constant" in exponential
    // damping.

    // test empty damper
    auto d = damper<Value, Sum, Count>(tc);
    static_assert(std::is_same_v<decltype(d.current()), Value>);
    EXPECT_EQ(d.current(), Value{});

    // tc samples with pre_step value
    for (Count i = 1;  i <= tc + 10;  ++i) {
        SCOPED_TRACE(std::to_string(i));
        auto dv = d(pre_step);

        EXPECT_EQ(dv, pre_step); // damped value mustn't change
    }

    // The damper is now fully "primed" with pre_step values.

    auto const maxerr = allowed_error(pre_step, post_step);
    auto const v0 = static_cast<double>(pre_step);
    auto const v1 = static_cast<double>(post_step);
    // tc further samples of `post_step` values
    for (Count i = 1;  i <= tc;  ++i) {
        SCOPED_TRACE(std::to_string(i));
        auto dv = static_cast<double>(d(post_step));
        double expected = reference_damping(v0, v1, tc, i);

        EXPECT_NEAR(dv, expected, maxerr);
    }
}

// signed integer types
TEST(Damper, Short) { test_damper<short>(); }
TEST(Damper, Int) { test_damper<int>(); }
TEST(Damper, Long) { test_damper<long>(); }
TEST(Damper, LongLong) { test_damper<long long>(); }

// signed integer types with a step from negative to positive
TEST(Damper, ShortNegstep) { test_damper<short>(-100, 100); }
TEST(Damper, IntNegstep) { test_damper<int>(-100, 100); }
TEST(Damper, LongNegstep) { test_damper<long>(-100, 100); }
TEST(Damper, LongLongNegstep) { test_damper<long long>(-100, 100); }

// unsigned integer types
TEST(Damper, UnsignedShort) { test_damper<unsigned short>(); }
TEST(Damper, Unsigned) { test_damper<unsigned>(); }
TEST(Damper, UnsignedLong) { test_damper<unsigned long>(); }
TEST(Damper, UnsignedLongLong) { test_damper<unsigned long long>(); }

// FP types
TEST(Damper, Float) { test_damper<float>(); }
TEST(Damper, Double) { test_damper<double>(); }
TEST(Damper, LongDouble) { test_damper<long double>(); }

// FP types with a step from negative to positive
TEST(Damper, FloatNegstep) { test_damper<float>(-100.0, 100.0); }
TEST(Damper, DoubleNegstep) { test_damper<double>(-100.0, 100.0); }
//TEST(damper, long_double_negstep) { test_damper<long double>(-100.0, 100.0); }

// tiny integer types (one use case is microcontroller analog input damping)
TEST(Damper, Uint8Uint16) { test_damper<std::uint8_t, std::uint16_t>(); }
TEST(Damper, Uint8Uint16Uint8) {
    test_damper<std::uint8_t, std::uint16_t, std::uint8_t>();
}

// Standard test fails with 8-bit sum which overflows.
// I have tried to algebraically eliminate the sum, but it seems that,
// for integer arithmetic, we always need at least a temporary result
// which can hold sample avg * time_constant.
//
// TEST(damper, uint8_t__uint8_t__uint8_t) {
//     test_damper<std::uint8_t, std::uint8_t, std::uint8_t>();
// }

// However, it is still a valid instantiation when sample avg * tc is small.
TEST(Damper, Uint8Uint8Uint8SmallValues) {
    test_damper<std::uint8_t, std::uint8_t, std::uint8_t>(0, 40, 5);
}
\$\endgroup\$
7
  • \$\begingroup\$ Thanks. Some interesting comments on the testing techniques, especially SCOPE_TRACE. Interesting that the long double tests failed for you - looks like "just more rounding error". I used clang-13 on linux. The earlier version of my tests (see CR hisoiry) had a separate fully templated function for the model. I chose not to test against the long double because I wanted to be sure that the casting didn't affect things. \$\endgroup\$ Commented Feb 17, 2022 at 12:56
  • \$\begingroup\$ The "signedness switch" in the template was there to address a specific problem, before I added the static_casts to suppress the final remaining warnings. Interestingly and confusingly, before I had the static cast and the signedness switch I was actually getting BADLY failing tests for the cases with step from -100 => 100. Like the answer was TOTALLY wrong, the dv seemed to overflowing a signed integer which is UB and the UBSan was kicking off. I put the type switch for Count in to solve this. BUT NOW when I take out the type switch and the static_casts it seems fine. I am not sure..? \$\endgroup\$ Commented Feb 17, 2022 at 13:01
  • \$\begingroup\$ Further thought on both your and my code: Should the concept constraints be applied to all three types: Value, Sum, and Count? And is it bad form to declare and use template <typename T> concept arithmetic = std::integral<T> || std::floating_point<T>; because the library people chose not to, for some obscure reason? \$\endgroup\$ Commented Feb 17, 2022 at 13:13
  • \$\begingroup\$ I think I refound my formerly failing test (the one that motivated the std::conditional on Count): TEST(Damper, IntNegstep) { test_damper<int>(-100, 100); } and then change the default type for Count to unsigned in both damper and test_damper and remove all the static casts. Boom... all the EXPECTs fail. gcc is silent, clang says: "implicit coversion changes signedness". The reason unsigned short works is due to the "auto promotion to int"? And of course a static cast fixes it. But still nasty. Especially on gcc with no warning. And I don't fully get it. Can you explain? \$\endgroup\$ Commented Feb 17, 2022 at 17:19
  • \$\begingroup\$ gcc needs an additional -Wsign-conversion becuase that is not included for integers under -Wconversion. It then warns similarly to clang... So that's a good learning point. Both compilers report signedness change twice. First sum is "promoted" from int to unsigned and then the result from sum_ / count_ is sign converted from unsigned to int. What is a clean way to "generically" fix this. Just clobber the problem with a static_cast? \$\endgroup\$ Commented Feb 17, 2022 at 17:34
0
\$\begingroup\$

The final choices I went with. See comments

damper.hpp

#pragma once

#include <cmath>
#include <concepts>
#include <cstdint>
#include <stdexcept>
#include <type_traits>

// Dampens a "noisy" value using "approx exponential damping".
// Submit samples using `operator()`, which returns the damped value.
// MUST choose a SUM type which can hold: sample avg * time_constant
template <typename Value, typename Sum = Value,
          typename Count = std::conditional_t<std::is_signed_v<Sum>, short, unsigned short>>
requires ((std::floating_point<Value> || std::integral<Value>) &&
          (std::floating_point<Sum> || std::integral<Sum>)
          && std::integral<Count>)
class damper {
  public:
    // having given the user a lot of choice on types, we now apply some fairly
    // strict rules to check their choices were sensible
    static_assert(std::is_signed_v<Value> == std::is_signed_v<Sum>,
                  "Value and Sum types must have the same signedness");

    static_assert(std::is_signed_v<Sum> == std::is_signed_v<Count>,
                  "Sum and Count types must have the same signedness");

    static_assert(std::is_floating_point_v<Value> == std::is_floating_point_v<Sum>,
                  "Sum and Value types must both be floating point or both be integral.");

    static_assert(sizeof(Sum) >= sizeof(Value), "Sum type should be at least as big as Value type");
    static_assert(sizeof(Sum) >= sizeof(Count), "Sum type should be at least as big as Count type");

    explicit damper(Count time_constant) : time_constant_(time_constant) {
        if (time_constant <= 0) throw std::invalid_argument("damper needs positive time constant");
    }

    Value operator()(Value sample) noexcept {
        if (std::isnan(sample)) return damped_value_; // ignore it

        if (count_ != time_constant_) {
            // branch avoids distortions for the first time_constant samples
            sum_ += sample;
            ++count_;
        } else {
            // Main running branch: once system is "primed" with time_constant number of samples
            // provides "approxinately exponential damping" with the given time constant

            // correct, and well defined, even with unsigned types
            // but is potential remains for overflow, if user choses a Sum type which cannot
            // hold sample avg * time_constant

            // working around gcc vageries with -Wconversion by doing this in two steps
            // https://stackoverflow.com/a/71183117/1087626
            sum_ += sample;
            sum_ -= damped_value_;
        }

        // only use the static_cast when it should be legitimately needed
        if constexpr (sizeof(Sum) > sizeof(Value)) {
            // static_cast reduces Sum back to Value if different. "Should fit" after the
            // division.
            damped_value_ = static_cast<Value>(sum_ / count_);
        } else {
            // this may cause a -Wconversion warning on gcc if -fsanitize=undefined is passed
            // https://stackoverflow.com/q/71181566/1087626
            damped_value_ = sum_ / count_;
        }
        return damped_value_;
    }

    [[nodiscard]] Value current() const { return damped_value_; }

    void reset() noexcept {
        sum_          = Sum{};
        count_        = Count{};
        damped_value_ = Value{};
    }

  private:
    Sum         sum_{};
    Value       damped_value_{};
    Count       count_{};
    const Count time_constant_;
};

damper_tests.cpp

#include "damper.hpp"
#include "gtest/gtest.h"
#include <cmath>
#include <concepts>
#include <cstdint>
#include <limits>
#include <string>
#include <type_traits>

// for testing purposes only
template <typename Value, typename Sum = Value,
          typename Count = std::conditional_t<std::is_signed_v<Sum>, short, unsigned short>>
void test_damper(Value pre_step = 100, Value post_step = 0, Count tc = 10) {

    // modelling a step response from Value{pre_step} down to Value{post_step} after tc samples
    // expecting flat Value{100}
    // then we feed another tc samples at the pre_step value to ensure "main branch" of `damper`
    // continues to return a flat profile
    // Then we start feeding `post_step` values and expect "approximately exponential decay" towards
    // Value{0}. `tc` is aprox equiv to the "time constant" in exponential damping

    // test empty damper
    auto d = damper<Value, Sum, Count>(tc);
    static_assert(std::is_same_v<decltype(d.current()), Value>);
    EXPECT_EQ(d.current(), Value{});

    // tc samples with pre_step value
    for (auto i = Count{1}; i <= tc; ++i) {
        auto dv = d(pre_step);
        SCOPED_TRACE("Stage 1: feeding " + std::to_string(i) + "th value of " +
                     std::to_string(pre_step));
        EXPECT_EQ(dv, pre_step); // check at every step that the `damped_value_` is "flat"
    }
    // `damper` is now fully "primed" with pre_step values

    // continue to push another tc pre_step values to ensure the output doesn't change
    for (auto i = Count{1}; i <= tc; ++i) {
        auto dv = d(pre_step);
        SCOPED_TRACE("Stage 2: feeding " + std::to_string(i) + "th value of " +
                     std::to_string(pre_step));
        EXPECT_EQ(dv, pre_step); // check at every step that the `damped_value_` is "flat"
    }

    // tc further samples of `post_step` values
    for (auto i = Count{1}; i <= tc; ++i) {
        auto dv = d(post_step);
        SCOPED_TRACE("Stage 3: feeding " + std::to_string(i) + "th value of " +
                     std::to_string(post_step));

        // check at every sample that we are matching the exponential decay curve
        // which predicts "approximately exponential damping"
        auto expected = post_step + (pre_step - post_step) * std::pow((tc - 1.0L) / tc, i);

        if constexpr (std::is_integral_v<Value>) {
            auto predicted = static_cast<Value>(std::round(expected));
            // allow an integer result to be "off-by-one" after rounding
            // std::abs can't be used here, because programm would be ill-formed for some types
            EXPECT_TRUE(dv - predicted == Value{0} || dv - predicted == Value{1} ||
                        predicted - dv == Value{1});
        } else {
            auto predicted = static_cast<Value>(expected);
            // for FP we expect the prediction to be within appropiately scaled epsilon
            EXPECT_LT(std::fabs(dv - predicted),
                      std::fabs(pre_step - post_step) * std::numeric_limits<Value>::epsilon());
        }
    }
}


// signed integer types
TEST(Damper, Short) { test_damper<short>(); }
TEST(Damper, Int) { test_damper<int>(); }
TEST(Damper, Long) { test_damper<long>(); }
TEST(Damper, LongLong) { test_damper<long long>(); }

// signed integer types with a step from negative to positive
TEST(Damper, ShortNegstep) { test_damper<short>(-100, 100); }
TEST(Damper, IntNegstep) { test_damper<int>(-100, 100); }
TEST(Damper, LongNegstep) { test_damper<long>(-100, 100); }
TEST(Damper, LongLongNegstep) { test_damper<long long>(-100, 100); }

// unsigned integer types
TEST(Damper, UnsignedShort) { test_damper<unsigned short>(); }
TEST(Damper, Unsigned) { test_damper<int>(-100, 100); }
TEST(Damper, UnsignedLong) { test_damper<unsigned long>(); }
TEST(Damper, UnsignedLongLong) { test_damper<unsigned long long>(); }

// FP types
TEST(Damper, Float) { test_damper<float>(); }
TEST(Damper, Double) { test_damper<double>(); }
TEST(Damper, LongDouble) { test_damper<long double>(); }

// FP types with a step from negative to positive
TEST(Damper, FloatNegstep) { test_damper<float>(-100.0, 100.0); }
TEST(Damper, DoubleNegstep) { test_damper<double>(-100.0, 100.0); }
TEST(Damper, LongDoubleNegstep) { test_damper<long double>(-100.0, 100.0); }

// tiny integer types (one use case is microcontroller analog input damping)
TEST(Damper, Uint8Uint16) { test_damper<std::uint8_t, std::uint16_t>(); }
TEST(Damper, Int8Int16) { test_damper<std::int8_t, std::int16_t>(); }

// this test fails, because the sum overflows
// I have tried to algebraically eliminate the sum, but it seems that, for integer arithmetic,
// we always need at least a temporary result which can hold sample avg * time_constant
// TEST(damper, uint8_t__uint8_t__uint8_t) { test_damper<std::uint8_t, std::uint8_t, std::uint8_t>(); }

// however, it is still a valid set of template params when sample avg * tc is small
TEST(Damper, Uint8Uint8SmallValues) {
    test_damper<std::uint8_t, std::uint8_t, std::uint8_t>(0, 40, 5);
}

TEST(Damper, Int8Int8Int8Negstep) { test_damper<std::int8_t, std::int8_t, std::int8_t>(-10, 10, 10); }
TEST(Damper, IntIntNegstep) { test_damper<int, int>(-100'000, 100'000, 10); }
TEST(Damper, Int8Int16Negstep) { test_damper<std::int8_t, std::int16_t>(-10, 10, short(10)); }

```
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.