tape_machine_components.hpp

Go to the documentation of this file.

#pragma once

 
#include <pulp/signal/character_delay/tape_physical.hpp>
#include <pulp/signal/ballistics_filter.hpp>
#include <pulp/signal/character_delay/primitives.hpp>
#include <pulp/signal/character_delay/tables.hpp>
#include <pulp/signal/character_delay/tape.hpp>
#include <pulp/signal/character_delay/tape_loss.hpp>
#include <pulp/signal/denormal.hpp>
#include <pulp/signal/units.hpp>
#include <algorithm>
#include <array>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <vector>
 
namespace pulp::signal {

enum class TapeCurve : std::uint8_t {
    nab,             
    iec_ccir,        
    cassette_type1,  
    cassette_type2,  
};

enum class TapeArchetype : std::uint8_t {
    ampex_350_440,  
    studer_a800,    
    cassette_deck,  
};

namespace tape {

inline constexpr std::array<double, 5> kSpeedsIps = {1.875, 3.75, 7.5, 15.0, 30.0};

inline constexpr double kMetresPerInch = 0.0254;

inline constexpr double kEqReferenceHz = 1000.0;

struct EqTimeConstants {
    double bass_s = 0.0;
    double treble_s = 0.0;
 
    bool has_bass_shelf() const noexcept { return bass_s > 0.0; }
};

inline EqTimeConstants eq_time_constants(TapeCurve curve, double speed_ips) noexcept {
    switch (curve) {
        case TapeCurve::cassette_type1: return {3180e-6, 120e-6};
        case TapeCurve::cassette_type2: return {3180e-6, 70e-6};
        case TapeCurve::iec_ccir:
            // CCIR adopted 70 µs (from 100 µs) in 1966; 35 µs at 15 ips and
            // 17.5 µs at 30 ips per the A800 technical data. No bass shelf at
            // any reel speed — the structural difference from NAB.
            if (speed_ips >= 22.5) return {0.0, 17.5e-6};
            if (speed_ips >= 11.25) return {0.0, 35e-6};
            return {0.0, 70e-6};
        case TapeCurve::nab:
        default:
            // The professional convention extends NAB to 30 ips by dropping the
            // bass shelf and halving the treble constant again.
            if (speed_ips >= 22.5) return {0.0, 17.5e-6};
            return {3180e-6, 50e-6};
    }
}

inline double eq_corner_hz(double tau_s) noexcept {
    return tau_s > 0.0 ? 1.0 / (2.0 * chardelay::kPi * tau_s) : 0.0;
}

inline double eq_record_response_db(const EqTimeConstants& tc, double hz) noexcept {
    auto magnitude = [&tc](double f) {
        const double w = 2.0 * chardelay::kPi * f;
        const double numerator = tc.has_bass_shelf() ? std::hypot(1.0, w * tc.bass_s) : 1.0;
        return numerator / std::hypot(1.0, w * tc.treble_s);
    };
    return 20.0 * std::log10(magnitude(hz) / magnitude(kEqReferenceHz));
}

inline double eq_playback_response_db(const EqTimeConstants& tc, double hz) noexcept {
    return -eq_record_response_db(tc, hz);
}
 
// ── Reproduce and record head gaps ────────────────────────────────────────

struct HeadGapGeometry {
    double reproduce_m = 3.5e-6;
    double record_m = 10e-6;
};

inline HeadGapGeometry head_gap_geometry(double speed_ips) noexcept {
    if (speed_ips >= 11.25) return {3.5e-6, 10e-6};   // studio 15 / 30 ips
    if (speed_ips >= 2.8125) return {5e-6, 12e-6};    // consumer reel 7.5 / 3.75 ips
    return {0.8e-6, 3e-6};                            // cassette 1.875 ips
}

inline double speed_m_per_s(double speed_ips) noexcept {
    return speed_ips * kMetresPerInch;
}

inline double gap_loss_magnitude(double hz, double gap_m, double speed_ips) noexcept {
    const double x = chardelay::kPi * gap_m * hz / speed_m_per_s(speed_ips);
    if (x < 1e-9) return 1.0;
    return std::abs(std::sin(x) / x);
}
 
inline double gap_loss_db(double hz, double gap_m, double speed_ips) noexcept {
    return 20.0 * std::log10(std::max(gap_loss_magnitude(hz, gap_m, speed_ips), 1e-12));
}

inline double gap_null_hz(double gap_m, double speed_ips) noexcept {
    return speed_m_per_s(speed_ips) / std::max(gap_m, 1e-12);
}
 
// ── The age axis this module adds on top of the reused one ────────────────

inline constexpr std::array<double, 3> kAgeAxis = {0.0, 0.5, 1.0};

inline constexpr std::array<double, 3> kAgePrintThroughDb = {-80.0, -55.0, -38.0};

inline constexpr std::array<double, 3> kAgeBiasDriftUnits = {0.0, 0.15, 0.4};
 
inline double age_print_through_db(double age01) noexcept {
    return chardelay::interpolate_knots(kAgeAxis, kAgePrintThroughDb,
                                        std::clamp(age01, 0.0, 1.0));
}
 
inline double age_bias_drift(double age01) noexcept {
    return chardelay::interpolate_knots(kAgeAxis, kAgeBiasDriftUnits,
                                        std::clamp(age01, 0.0, 1.0));
}

inline double first_order_shelf_db(double corner_hz, double gain_linear, double hz,
                                   double sample_rate) noexcept {
    const double fc = std::clamp(corner_hz, 0.1, 0.49 * sample_rate);
    const double t = std::tan(chardelay::kPi * fc / sample_rate);
    const double g = t / (1.0 + t);
    const double theta = 2.0 * chardelay::kPi * hz / sample_rate;
    const std::complex<double> z_inv(std::cos(theta), -std::sin(theta));
    const std::complex<double> lp = g * (1.0 + z_inv) / (1.0 - (1.0 - 2.0 * g) * z_inv);
    const std::complex<double> shelf = gain_linear + (1.0 - gain_linear) * lp;
    return 20.0 * std::log10(std::max(std::abs(shelf), 1e-12));
}
 
// ── Archetype presets ─────────────────────────────────────────────────────

struct ArchetypePreset {
    std::array<double, 3> legal_speeds_ips{};
    int legal_speed_count = 0;
    double default_speed_ips = 15.0;
    TapeCurve default_curve = TapeCurve::nab;
    double crosstalk_db = -40.0;
    bool companding = false;
    double age01 = 0.0;
};

inline ArchetypePreset archetype_preset(TapeArchetype archetype) noexcept {
    ArchetypePreset p;
    switch (archetype) {
        case TapeArchetype::studer_a800:
            p.legal_speeds_ips = {7.5, 15.0, 30.0};
            p.legal_speed_count = 3;
            p.default_speed_ips = 15.0;
            p.default_curve = TapeCurve::nab;
            p.crosstalk_db = -40.0;
            p.companding = false;
            p.age01 = 0.05;
            return p;
        case TapeArchetype::cassette_deck:
            p.legal_speeds_ips = {1.875, 0.0, 0.0};
            p.legal_speed_count = 1;
            p.default_speed_ips = 1.875;
            p.default_curve = TapeCurve::cassette_type1;
            p.crosstalk_db = -30.0;
            // Companding on by default: dbx-equipped decks were a documented
            // cassette-era convention, and the reel machines of the other two
            // classes predate consumer companding entirely.
            p.companding = true;
            p.age01 = 0.3;
            return p;
        case TapeArchetype::ampex_350_440:
        default:
            p.legal_speeds_ips = {7.5, 15.0, 0.0};
            p.legal_speed_count = 2;
            p.default_speed_ips = 15.0;
            p.default_curve = TapeCurve::nab;
            p.crosstalk_db = -38.0;
            p.companding = false;
            p.age01 = 0.2;
            return p;
    }
}

inline double snap_speed_ips(TapeArchetype archetype, double requested_ips) noexcept {
    const ArchetypePreset preset = archetype_preset(archetype);
    double best = preset.legal_speeds_ips[0];
    double best_distance = std::abs(std::log(std::max(requested_ips, 1e-6) / best));
    for (int i = 1; i < preset.legal_speed_count; ++i) {
        const double candidate = preset.legal_speeds_ips[static_cast<std::size_t>(i)];
        // Nearest in RATIO, not in difference: the ladder is geometric, so
        // 5 ips is equidistant from 3.75 and 7.5 the way an ear hears it, and a
        // linear metric would bias every midpoint upward.
        const double distance = std::abs(std::log(std::max(requested_ips, 1e-6) / candidate));
        if (distance < best_distance) {
            best = candidate;
            best_distance = distance;
        }
    }
    return best;
}
 
}  // namespace tape

class TapeEqSection {
public:
    static constexpr double kInversePoleDamping = 0.90;

    static constexpr double kPrewarpCeiling = 0.49;
 
    void set(const tape::EqTimeConstants& tc, double sample_rate, bool inverse) noexcept {
        const double fs = sample_rate > 0.0 ? sample_rate : 48000.0;
        const double treble_hz =
            std::min(tape::eq_corner_hz(tc.treble_s), kPrewarpCeiling * fs);
        const double u_treble = std::tan(chardelay::kPi * treble_hz / fs);
 
        double b0 = 0.0, b1 = 0.0;
        const double a1 = (u_treble - 1.0) / (u_treble + 1.0);
        if (tc.has_bass_shelf()) {
            const double bass_hz =
                std::min(tape::eq_corner_hz(tc.bass_s), kPrewarpCeiling * fs);
            const double u_bass = std::tan(chardelay::kPi * bass_hz / fs);
            const double scale = (u_treble / u_bass) / (u_treble + 1.0);
            b0 = scale * (u_bass + 1.0);
            b1 = scale * (u_bass - 1.0);
        } else {
            b0 = u_treble / (u_treble + 1.0);
            b1 = b0;
        }
 
        normalize(b0, b1, a1, fs);
 
        if (inverse) {
            const double inv = 1.0 / b0;
            const double damping = tc.has_bass_shelf() ? 1.0 : kInversePoleDamping;
            const double new_b0 = inv;
            const double new_b1 = a1 * inv;
            const double new_a1 = damping * b1 * inv;
            b0 = new_b0;
            b1 = new_b1;
            normalize(b0, b1, new_a1, fs);
            b0_ = b0;
            b1_ = b1;
            a1_ = new_a1;
        } else {
            b0_ = b0;
            b1_ = b1;
            a1_ = a1;
        }
    }
 
    void reset() noexcept {
        x1_ = 0.0;
        y1_ = 0.0;
    }
 
    double process(double x) noexcept {
        const double y = b0_ * x + b1_ * x1_ - a1_ * y1_;
        x1_ = x;
        y1_ = snap_to_zero(y);
        return y;
    }

    double response_db(double hz, double sample_rate) const noexcept {
        return 20.0 * std::log10(std::max(magnitude(b0_, b1_, a1_, hz, sample_rate), 1e-12));
    }
 
private:
    static double magnitude(double b0, double b1, double a1, double hz,
                            double sample_rate) noexcept {
        const double theta = 2.0 * chardelay::kPi * hz / sample_rate;
        const double c = std::cos(theta), s = std::sin(theta);
        // e^{-jθ} = c − js.
        const double num = std::hypot(b0 + b1 * c, -b1 * s);
        const double den = std::hypot(1.0 + a1 * c, -a1 * s);
        return den > 1e-30 ? num / den : 0.0;
    }
 
    static void normalize(double& b0, double& b1, double a1, double fs) noexcept {
        const double gain = magnitude(b0, b1, a1, tape::kEqReferenceHz, fs);
        if (gain > 1e-30) {
            b0 /= gain;
            b1 /= gain;
        }
    }
 
    double b0_ = 1.0, b1_ = 0.0, a1_ = 0.0;
    double x1_ = 0.0, y1_ = 0.0;
};

template <typename SampleType = float>
class TapeCompanderT {
public:
    static constexpr double kAttackMs = 1.0;

    static constexpr double kReleaseMs = 100.0;

    static constexpr double kRefLevelDbfs = -20.0;
 
    void prepare(double sample_rate) {
        for (auto* follower : {&encode_env_, &decode_env_}) {
            follower->prepare(static_cast<SampleType>(sample_rate));
            follower->set_attack_ms(static_cast<SampleType>(kAttackMs));
            follower->set_release_ms(static_cast<SampleType>(kReleaseMs));
        }
        reference_ = units::db_to_linear(kRefLevelDbfs);
        reset();
    }
 
    void reset() noexcept {
        encode_env_.reset();
        decode_env_.reset();
    }

    double encode(double x) noexcept { return x * gain(encode_env_, x, -0.5); }

    double decode(double x) noexcept { return x * gain(decode_env_, x, +1.0); }
 
private:
    double gain(BallisticsFilterT<SampleType>& follower, double x, double slope) noexcept {
        const auto envelope =
            static_cast<double>(follower.process(static_cast<SampleType>(x)));
        const double level_db = units::linear_to_db(envelope / reference_);
        return units::db_to_linear(slope * level_db);
    }
 
    BallisticsFilterT<SampleType> encode_env_{};
    BallisticsFilterT<SampleType> decode_env_{};
    double reference_ = 0.1;
};
 
}  // namespace pulp::signal