analog_filter.hpp

#ifndef CONSTFILT_ANALOG_FILTER_HPP
#define CONSTFILT_ANALOG_FILTER_HPP
 
#include "discretize.hpp"
#include "filter.hpp"
#include "stability.hpp"
 
namespace constfilt
{
 
struct LowPass
{
};
 
struct HighPass
{
};
 
// Discretize an analog (continuous-time, s-domain) transfer function into a
// digital Filter.
//
// Coefficients are in descending power order:
//   coeff[0]*s^N + coeff[1]*s^{N-1} + ... + coeff[N]
//
// Template parameters:
//   T         - numeric type (float, double, ...)
//   N         - filter order (degree of denominator)
//   Method    - TustinNW (default), TustinPW, ZOH, or MatchedZ
//   CheckStab - reserved for future use; currently has no effect (see below).
//
// Constructors:
//   AnalogFilter(b_c, a_c, sample_rate_hz)
//     b_c            - s-domain numerator array [N+1], descending order
//     a_c            - s-domain denominator array [N+1], descending order
//     sample_rate_hz - sample rate in Hz
//
//   AnalogFilter(b_c, a_c, sample_rate_hz, method_tag)
//     method_tag     - method instance; for TustinPW supply
//                      constfilt::prewarp(warp_hz)
//
//   AnalogFilter(continuous_tf, sample_rate_hz)
//     continuous_tf  - s-domain TransferFunction (e.g. from a subclass)
//     sample_rate_hz - sample rate in Hz
//
//   AnalogFilter(continuous_tf, sample_rate_hz, method_tag)
//     method_tag     - method instance; for TustinPW supply
//                      constfilt::prewarp(warp_hz)
//
//   AnalogFilter(continuous_tf, factored_tf, sample_rate_hz, method_tag)
//     factored_tf    - the poles, zeros, and gain of continuous_tf in
//                      factored form. Required alongside continuous_tf.
//                      For ZOH, the poles are used to build the Vandermonde
//                      eigenvector matrix analytically. For MatchedZ, the
//                      poles and zeros are mapped directly to the z-domain.
template <typename T, consteig::Size N, typename Method = TustinNW,
          bool CheckStab = true>
class AnalogFilter : public Filter<T, N + 1u, N + 1u>
{
    static_assert(N >= 1u, "Filter order must be at least 1");
 
    using BoundMethod = typename bind_method<T, Method>::type;
 
  public:
    constexpr AnalogFilter(const T (&b_c)[N + 1u], const T (&a_c)[N + 1u],
                           T sample_rate_hz)
        : AnalogFilter(
              checked_discretize(b_c, a_c, sample_rate_hz, BoundMethod{}))
    {
        static_assert(!is_tustinpw_tag<Method>::value,
                      "TustinPW requires a warp frequency; use the "
                      "method-tag constructor with constfilt::prewarp(wc_hz)");
    }
 
    constexpr AnalogFilter(const T (&b_c)[N + 1u], const T (&a_c)[N + 1u],
                           T sample_rate_hz, BoundMethod method_tag)
        : AnalogFilter(checked_discretize(b_c, a_c, sample_rate_hz, method_tag))
    {
    }
 
    constexpr AnalogFilter(TransferFunction<T, N + 1u, N + 1u> continuous_tf,
                           T sample_rate_hz)
        : AnalogFilter(checked_discretize(continuous_tf.b, continuous_tf.a,
                                          sample_rate_hz, BoundMethod{}))
    {
        static_assert(!is_tustinpw_tag<Method>::value,
                      "TustinPW requires a warp frequency; use the "
                      "method-tag constructor with constfilt::prewarp(wc_hz)");
    }
 
    constexpr AnalogFilter(TransferFunction<T, N + 1u, N + 1u> continuous_tf,
                           T sample_rate_hz, BoundMethod method_tag)
        : AnalogFilter(checked_discretize(continuous_tf.b, continuous_tf.a,
                                          sample_rate_hz, method_tag))
    {
    }
 
    constexpr AnalogFilter(TransferFunction<T, N + 1u, N + 1u> continuous_tf,
                           const FactoredTF<T, N> &factored_tf,
                           T sample_rate_hz, BoundMethod method_tag)
        : AnalogFilter(checked_discretize_factored(continuous_tf.b,
                                                   continuous_tf.a, factored_tf,
                                                   sample_rate_hz, method_tag))
    {
    }
 
  private:
    constexpr explicit AnalogFilter(
        TransferFunction<T, N + 1u, N + 1u> digital_tf)
        : Filter<T, N + 1u, N + 1u>(digital_tf.b, digital_tf.a)
    {
    }
 
    static constexpr TransferFunction<T, N + 1u, N + 1u> checked_discretize(
        const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], T sample_rate_hz,
        BoundMethod method_tag)
    {
        // Stability check is currently disabled. See:
        // https://github.com/MitchellThompkins/constfilt/issues/14
        //
        // The original implementation used:
        //   throw "constfilt: unstable analog filter";
        // which serves dual purpose in C++17: a runtime exception and, when
        // evaluated in a constexpr context, a compile-time error. This is
        // incompatible with -fno-exceptions and freestanding (no-stdlib)
        // builds.
        //
        // Candidate replacements, each with trade-offs:
        //   __builtin_trap()          - preserves compile-time error via the
        //                               non-constexpr call rule, but is
        //                               GCC/Clang specific.
        //   [[noreturn]] customization point (declared, not defined) -
        //   portable,
        //                               user supplies the handler; linker error
        //                               if forgotten; compile-time error still
        //                               works via non-constexpr call rule.
        //   Remove entirely           - current state; no check at runtime or
        //                               compile time; caller's responsibility.
        //
        // if (CheckStab &&
        //     check_stability(tf_to_ss<T, N>(b_c, a_c)) == Stability::Unstable)
        // {
        //     throw "constfilt: unstable analog filter";
        // }
        return analog_to_digital<T, N>(
            b_c, a_c, static_cast<T>(1) / sample_rate_hz, method_tag);
    }
 
    static constexpr TransferFunction<T, N + 1u, N + 1u>
    checked_discretize_factored(const T (&b_c)[N + 1u], const T (&a_c)[N + 1u],
                                const FactoredTF<T, N> &factored_tf,
                                T sample_rate_hz, BoundMethod method_tag)
    {
        return discretize_with_factored<T, N>(
            b_c, a_c, factored_tf, static_cast<T>(1) / sample_rate_hz,
            method_tag);
    }
};
 
} // namespace constfilt
 
#endif // CONSTFILT_ANALOG_FILTER_HPP