thiran_allpass.hpp

Go to the documentation of this file.

// sbl/dsp/primitives/thiran_allpass.hpp — Maximally-flat group-delay allpass
//
// Fractional delay for a fixed delay amount: the phase delay is flat to
// order 2N at DC, so a waveguide loop can be tuned to a fraction of a sample
// without the comb-filtering a linear interpolator introduces. Hermite
// interpolation (DelayLine::read_cubic) remains the tool for a delay that
// moves per sample; this is the tool for one that is set and held
// (RPT-029 Part 3).
//
// The delay line carries the integer part; this carries the fraction. Total
// phase delay D is clamped to [N - 0.5, N + 0.5], the range where the design
// is both accurate and stable (Thiran is stable for D > N - 1). D = N gives
// every coefficient as zero — a pure N-sample delay, exactly.
//
// Energy: an allpass has unit magnitude at every frequency, so it stores and
// returns energy but never creates it (RPT-029 Part 4 rule 1). Its state is
// a delay line's worth of past input and output; the energy function is the
// usual sum of squares over that state.
//
// Usage:
//   sbl::dsp::primitives::ThiranAllpass<1> tuner;
//   tuner.set_delay(1.0f + frac);       // frac in [-0.5, 0.5]
//   float y = tuner.process(x);
//   float d = tuner.group_delay_samples();
 
#ifndef SBL_DSP_PRIMITIVES_THIRAN_ALLPASS_HPP_
#define SBL_DSP_PRIMITIVES_THIRAN_ALLPASS_HPP_
 
#include <cstdint>
 
#include <sbl/dsp/math/clamp.hpp>
 
namespace sbl::dsp::primitives {
 
namespace detail {
 
/// Binomial coefficient C(n, k) for the small orders this filter allows.
inline constexpr float binomial(uint8_t n, uint8_t k) {
    float result = 1.0f;
    for (uint8_t i = 0; i < k; ++i) {
        result *= static_cast<float>(n - i) / static_cast<float>(i + 1);
    }
    return result;
}
 
} // namespace detail
 
/**
 * @brief Order-N Thiran allpass: fractional delay by phase, not interpolation
 *
 * @tparam N Filter order; the delay it can express is N ± 0.5 samples.
 *
 * @note All public methods are ISR-safe — bounded computation, no I/O.
 */
template<uint8_t N = 1>
class ThiranAllpass {
    static_assert(N >= 1, "Thiran allpass needs at least first order");
    static_assert(N <= 4, "high orders cost accuracy in float; use a longer delay line");
 
public:
    ThiranAllpass() { set_delay(static_cast<float>(N)); }
 
    /**
     * @brief Set the total phase delay in samples
     *
     * Clamped to [N - 0.5, N + 0.5]. Coefficients are recomputed here, so
     * call it when the pitch changes (block rate), not per sample.
     */
    void set_delay(float delay_samples) {
        constexpr float lo = static_cast<float>(N) - 0.5f;
        constexpr float hi = static_cast<float>(N) + 0.5f;
        const float d = math::clamp(delay_samples, lo, hi);
        delay_ = d;
 
        // a_k = (-1)^k C(N,k) prod_{n=0..N} (D - N + n) / (D - N + k + n)
        for (uint8_t k = 1; k <= N; ++k) {
            float coeff = detail::binomial(N, k);
            for (uint8_t n = 0; n <= N; ++n) {
                coeff *= (d - static_cast<float>(N) + static_cast<float>(n)) /
                         (d - static_cast<float>(N) + static_cast<float>(k + n));
            }
            a_[k - 1] = (k & 1) ? -coeff : coeff;
        }
    }
 
    /**
     * @brief Process one sample
     *
     * Numerator coefficients are the denominator's, reversed — that is what
     * makes the magnitude response unity.
     */
    float process(float x) {
        float y = a_[N - 1] * x;                       // b_0 = a_N
        for (uint8_t k = 1; k < N; ++k) {
            y += a_[N - 1 - k] * x_hist_[k - 1];       // b_k = a_(N-k)
        }
        y += x_hist_[N - 1];                           // b_N = a_0 = 1
        for (uint8_t k = 1; k <= N; ++k) {
            y -= a_[k - 1] * y_hist_[k - 1];
        }
 
        for (uint8_t i = N - 1; i > 0; --i) {
            x_hist_[i] = x_hist_[i - 1];
            y_hist_[i] = y_hist_[i - 1];
        }
        x_hist_[0] = x;
        y_hist_[0] = y;
        return y;
    }
 
    /// Phase delay the filter is tuned to, in samples (flat at DC by design).
    float group_delay_samples() const { return delay_; }
 
    void reset() {
        for (uint8_t i = 0; i < N; ++i) {
            x_hist_[i] = 0.0f;
            y_hist_[i] = 0.0f;
        }
    }
 
    static constexpr uint8_t order() { return N; }
 
    /// Smallest and largest delay this order can express.
    static constexpr float min_delay() { return static_cast<float>(N) - 0.5f; }
    static constexpr float max_delay() { return static_cast<float>(N) + 0.5f; }
 
private:
    float a_[N]{};            ///< a_1 .. a_N
    float x_hist_[N]{};       ///< x[n-1] .. x[n-N]
    float y_hist_[N]{};       ///< y[n-1] .. y[n-N]
    float delay_ = static_cast<float>(N);
};
 
} // namespace sbl::dsp::primitives
 
#endif // SBL_DSP_PRIMITIVES_THIRAN_ALLPASS_HPP_