bridged_t_resonator.hpp

Go to the documentation of this file.

#pragma once
 
#include <pulp/signal/denormal.hpp>
 
#include <algorithm>
#include <cmath>
 
namespace pulp::signal {

struct BridgedTComponents {
    double c41 = 15e-9;      
    double c42 = 15e-9;      
    double r161 = 1.0e6;     
    double r165 = 47.0e3;    
    double r166 = 6.8e3;     
    double r167 = 1.0e6;     
    double r170 = 470.0e3;   
};

struct Q43Leakage {
    double alpha = 14.3150;  
    double v0 = -0.5560;     
    double m = 1.4765e-5;    
};

inline double q43_collector_current(double vcomm, const Q43Leakage& leak) noexcept {
    const double e = -leak.alpha * (vcomm - leak.v0);
    const double scale = leak.m / leak.alpha;
    if (e > 30.0) return -e * scale;                 // log1p(exp(e)) == e here
    if (e < -30.0) return -std::exp(e) * scale;      // log1p(x) == x here
    return -std::log1p(std::exp(e)) * scale;
}

inline double q43_branch_resistance(double vcomm,
                                    const BridgedTComponents& c,
                                    const Q43Leakage& leak,
                                    double guard_v = -0.01) noexcept {
    const double series = c.r165 + c.r166;
    if (vcomm > guard_v) return series;
    const double ic = q43_collector_current(vcomm, leak);
    const double den = vcomm + c.r165 * ic;
    if (den == 0.0) return series;
    return std::clamp(series * vcomm / den, c.r166, series);
}

inline double bridged_t_shunt_resistance(double r_e,
                                         const BridgedTComponents& c) noexcept {
    return 1.0 / (1.0 / c.r161 + 1.0 / c.r170 + 1.0 / r_e);
}

inline double bridged_t_center_frequency(double r_eff,
                                         const BridgedTComponents& c) noexcept {
    constexpr double two_pi = 6.28318530717958647692;
    return 1.0 / (two_pi * std::sqrt(r_eff * c.r167 * c.c41 * c.c42));
}

inline BridgedTComponents bridged_t_retuned_to(const BridgedTComponents& c,
                                               double target_hz) noexcept {
    constexpr double kMinFarads = 1.0e-10;
    constexpr double kMaxFarads = 1.0e-6;
    if (!(target_hz > 0.0)) return c;
    const double r_rest = c.r165 + c.r166;
    const double current =
        bridged_t_center_frequency(bridged_t_shunt_resistance(r_rest, c), c);
    if (!(current > 0.0)) return c;
 
    BridgedTComponents out = c;
    const double scale = current / target_hz;
    out.c41 = std::clamp(c.c41 * scale, kMinFarads, kMaxFarads);
    out.c42 = std::clamp(c.c42 * scale, kMinFarads, kMaxFarads);
    return out;
}
 
inline double bridged_t_q(double r_eff, const BridgedTComponents& c) noexcept {
    const double alpha2 = r_eff * c.r167 * c.c41 * c.c42;
    const double alpha1 = r_eff * (c.c41 + c.c42);
    return std::sqrt(alpha2) / alpha1;
}

class BridgedTResonator {
public:
    struct Output {
        double vbt = 0.0;    
        double vcomm = 0.0;  
    };
 
    void set_components(const BridgedTComponents& c) noexcept { comps_ = c; }
    const BridgedTComponents& components() const noexcept { return comps_; }
 
    void set_leakage(const Q43Leakage& l) noexcept { leak_ = l; }
    const Q43Leakage& leakage() const noexcept { return leak_; }

    void set_sigh_enabled(bool enabled) noexcept { sigh_enabled_ = enabled; }
    bool sigh_enabled() const noexcept { return sigh_enabled_; }

    void set_rail_voltage(double v) noexcept { rail_v_ = v; }
 
    void prepare(double sample_rate) noexcept {
        h_ = 1.0 / sample_rate;
        reset();
    }

    void reset() noexcept {
        x1_ = 0.0;
        x2_ = 0.0;
        vcomm_delayed_ = 0.0;
        r_e_ = comps_.r165 + comps_.r166;
    }

    void set_attack_shunt(bool active) noexcept { attack_shunt_ = active; }

    Output process(double v_plus, double v_fb, double v_rp) noexcept {
        // Shift logic reads a delayed Vcomm, which is what breaks the
        // Vcomm -> iC -> R_eff -> Vcomm loop into something computable.
        if (attack_shunt_) {
            r_e_ = comps_.r166;
        } else if (!sigh_enabled_) {
            r_e_ = comps_.r165 + comps_.r166;
        } else {
            r_e_ = q43_branch_resistance(vcomm_delayed_, comps_, leak_);
        }
 
        const double g = 1.0 / bridged_t_shunt_resistance(r_e_, comps_);
 
        // States are the capacitor voltages: x1 = v_C41 = V+ - Vcomm,
        // x2 = v_C42 = Vbt - Vcomm. Only a21 moves with R_eff.
        const double inv_r167_c41 = 1.0 / (comps_.r167 * comps_.c41);
        const double a11 = -inv_r167_c41;
        const double a12 = inv_r167_c41;
        const double a21 = (1.0 / comps_.r167 - g) / comps_.c42;
        const double a22 = -1.0 / (comps_.r167 * comps_.c42);
        const double u2 = (g * v_plus - v_rp / comps_.r161 - v_fb / comps_.r170) / comps_.c42;
 
        const double hh = 0.5 * h_;
        const double m11 = 1.0 - hh * a11;
        const double m12 = -hh * a12;
        const double m21 = -hh * a21;
        const double m22 = 1.0 - hh * a22;
        const double r1 = (1.0 + hh * a11) * x1_ + hh * a12 * x2_;
        const double r2 = hh * a21 * x1_ + (1.0 + hh * a22) * x2_ + h_ * u2;
        const double det = m11 * m22 - m12 * m21;
 
        x1_ = (m22 * r1 - m12 * r2) / det;
        x2_ = (m11 * r2 - m21 * r1) / det;
 
        Output out;
        out.vcomm = v_plus - x1_;
        out.vbt = std::clamp(out.vcomm + x2_, -rail_v_, rail_v_);
        vcomm_delayed_ = out.vcomm;
        return out;
    }

    void snap_denormals() noexcept {
        x1_ = snap_to_zero(x1_);
        x2_ = snap_to_zero(x2_);
        vcomm_delayed_ = snap_to_zero(vcomm_delayed_);
    }

    double branch_resistance() const noexcept { return r_e_; }

    double center_frequency_hz() const noexcept {
        return bridged_t_center_frequency(bridged_t_shunt_resistance(r_e_, comps_), comps_);
    }
 
private:
    BridgedTComponents comps_{};
    Q43Leakage leak_{};
    double h_ = 1.0 / 48000.0;
    double rail_v_ = 15.0;
    bool sigh_enabled_ = true;
    bool attack_shunt_ = false;
 
    double x1_ = 0.0;
    double x2_ = 0.0;
    double vcomm_delayed_ = 0.0;
    double r_e_ = 47.0e3 + 6.8e3;
};
 
}  // namespace pulp::signal