tape_machine.hpp

Go to the documentation of this file.

#pragma once

 
// `character_delay/tape_physical.hpp` MUST come before `tape_loss.hpp` here,
// which is why this block is not alphabetized. `hysteresis.hpp` (pulled in by
// tape_physical) names `detail::LinearPhaseOversamplingStage2x` unqualified from
// inside `pulp::signal::chardelay`, and that only resolves to
// `pulp::signal::detail` while `pulp::signal::chardelay::detail` does not exist
// yet. `tape_loss.hpp` opens `chardelay::detail`, so including it first makes the
// nearer namespace win and the build fails on a name that is plainly there.
#include <pulp/signal/tape_machine_components.hpp>
#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 {

template <typename SampleType = float>
class TapeMachineT {
public:
    using Curve = TapeCurve;
    using Archetype = TapeArchetype;
 
    // ── Design parameters (the complete roster) ───────────────────────────

    static constexpr double kBiasMin = -1.0;
    static constexpr double kBiasMax = 1.0;
    static constexpr double kBiasDefault = 0.0;

    static constexpr double kOverBiasHfCutDb = 6.0;
    static constexpr double kUnderBiasHfBoostDb = 4.0;
    static constexpr double kBiasShelfHz = 8000.0;

    static constexpr double kBiasDistK = 0.6;

    static constexpr double kBiasErasureDb = 3.0;

    static constexpr double kDriveDbMax = 24.0;
    static constexpr double kDriveDefault = 0.3;

    static constexpr double kCrosstalkDbMin = -45.0;
    static constexpr double kCrosstalkDbMax = -22.0;

    static constexpr double kCrosstalkTiltHz = 5000.0;
    static constexpr double kCrosstalkTiltMaxDb = 6.0;

    static constexpr double kPrintThroughDbMin = -80.0;
    static constexpr double kPrintThroughDbMax = -38.0;

    static constexpr double kPrintThroughOffsetMsDefault = 700.0;
    static constexpr double kPrintThroughOffsetMsMin = 200.0;
    static constexpr double kPrintThroughOffsetMsMax = 2000.0;

    static constexpr double kPreEchoAsymmetryDb = 4.0;

    static constexpr double kAlignmentCeilingDb = 24.0;

    static constexpr double kInstabilityNominalMs = 2.0;

    static int oversampler_latency_samples() noexcept {
        return chardelay::TapePhysicalChannel::oversampler_latency_samples();
    }
 
    // ── Lifecycle ─────────────────────────────────────────────────────────
 
    TapeMachineT() { set_archetype(Archetype::ampex_350_440); }

    void prepare(double sample_rate) {
        if (std::isfinite(sample_rate) && sample_rate > 0.0)
            sample_rate_ = sample_rate;
 
        for (int index = 0; index < 2; ++index)
            channels_[static_cast<std::size_t>(index)].prepare(sample_rate_, index);
 
        instability_nominal_samples_ = static_cast<int>(std::llround(
            units::ms_to_samples(kInstabilityNominalMs, sample_rate_)));
        print_max_samples_ = static_cast<int>(std::llround(
            units::ms_to_samples(kPrintThroughOffsetMsMax, sample_rate_)));
 
        compander_l_.prepare(sample_rate_);
        compander_r_.prepare(sample_rate_);
 
        update_speed(speed_ips_, /*redesign=*/true);
        update_eq();
        apply_age();
        update_crosstalk();
        update_print_through();
        reset();
    }

    void reset() {
        for (auto& channel : channels_) channel.reset();
        compander_l_.reset();
        compander_r_.reset();
    }
 
    // ── Parameters ────────────────────────────────────────────────────────

    void set_archetype(Archetype archetype) {
        archetype_ = archetype;
        const tape::ArchetypePreset preset = tape::archetype_preset(archetype);
        curve_ = preset.default_curve;
        crosstalk_db_ = preset.crosstalk_db;
        companding_ = preset.companding;
        age01_ = preset.age01;
        set_speed_ips(preset.default_speed_ips);
        update_eq();
        apply_age();
        update_crosstalk();
        update_print_through();
    }
 
    Archetype archetype() const noexcept { return archetype_; }

    void set_speed_ips(double ips) {
        if (!std::isfinite(ips)) return;
        update_speed(tape::snap_speed_ips(archetype_, ips), /*redesign=*/false);
        update_eq();
    }
 
    double speed_ips() const noexcept { return speed_ips_; }

    void set_eq_curve(Curve curve) {
        curve_ = curve;
        update_eq();
    }
 
    Curve eq_curve() const noexcept { return curve_; }
 
    void set_bias(float bias_c) {
        if (!std::isfinite(static_cast<double>(bias_c))) return;
        bias_request_ = std::clamp(static_cast<double>(bias_c), kBiasMin, kBiasMax);
        update_bias();
        update_print_through();
    }

    double effective_bias() const noexcept { return bias_; }
 
    void set_drive(float drive01) {
        if (!std::isfinite(static_cast<double>(drive01))) return;
        drive01_ = std::clamp(static_cast<double>(drive01), 0.0, 1.0);
        update_bias();
    }
 
    double drive() const noexcept { return drive01_; }
 
    void set_age(float age01) {
        if (!std::isfinite(static_cast<double>(age01))) return;
        age01_ = std::clamp(static_cast<double>(age01), 0.0, 1.0);
        apply_age();
        update_print_through();
    }
 
    double age() const noexcept { return age01_; }
 
    void set_crosstalk_db(float db) {
        if (!std::isfinite(static_cast<double>(db))) return;
        crosstalk_db_ = std::clamp(static_cast<double>(db), kCrosstalkDbMin, kCrosstalkDbMax);
        update_crosstalk();
    }
 
    double crosstalk_db() const noexcept { return crosstalk_db_; }
 
    void set_companding_enabled(bool enabled) noexcept { companding_ = enabled; }
    bool companding_enabled() const noexcept { return companding_; }

    void set_print_through(float level_db, float offset_ms, bool pre_echo_enabled) {
        if (!std::isfinite(static_cast<double>(level_db)) ||
            !std::isfinite(static_cast<double>(offset_ms)))
            return;
        print_through_db_ =
            std::clamp(static_cast<double>(level_db), kPrintThroughDbMin, kPrintThroughDbMax);
        print_offset_ms_ = std::clamp(static_cast<double>(offset_ms),
                                      kPrintThroughOffsetMsMin, kPrintThroughOffsetMsMax);
        pre_echo_ = pre_echo_enabled;
        update_print_through();
    }
 
    double print_through_db() const noexcept { return print_through_db_; }
    double print_offset_ms() const noexcept { return print_offset_ms_; }
    bool pre_echo_enabled() const noexcept { return pre_echo_; }
 
    void set_mix(float wet01) noexcept {
        if (!std::isfinite(static_cast<double>(wet01))) return;
        mix_ = std::clamp(static_cast<double>(wet01), 0.0, 1.0);
    }
 
    // ── Audio ─────────────────────────────────────────────────────────────
 
    void process(const SampleType* in_l, const SampleType* in_r, SampleType* out_l,
                 SampleType* out_r, int n_frames) {
        for (int n = 0; n < n_frames; ++n) {
            const double dry_l = static_cast<double>(in_l[n]);
            const double dry_r = static_cast<double>(in_r[n]);
 
            // Inputs immediately enter several recursive stages. Treat a
            // non-finite frame as a recovery boundary before it can latch
            // those states; reset is allocation-free by contract.
            if (!std::isfinite(dry_l) || !std::isfinite(dry_r)) {
                recover_audio_fault();
                out_l[n] = SampleType{0};
                out_r[n] = SampleType{0};
                continue;
            }
 
            // Record chain. Companding first, so the compressor sees the
            // programme rather than the tape's own colouration.
            double left = companding_ ? compander_l_.encode(dry_l) : dry_l;
            double right = companding_ ? compander_r_.encode(dry_r) : dry_r;
 
            left = channels_[0].record(left, drive_gain_, inverse_drive_gain_);
            right = channels_[1].record(right, drive_gain_, inverse_drive_gain_);
 
            // Playback chain up to the point where the two channels meet.
            left = channels_[0].playback(left);
            right = channels_[1].playback(right);
 
            // Crosstalk: each channel leaks a tilted copy of the other. Read
            // both before writing either, or the second leak would be computed
            // from a channel that already contains the first.
            const double leak_l = crosstalk_gain_ * channels_[0].crosstalk_tilt.process(right);
            const double leak_r = crosstalk_gain_ * channels_[1].crosstalk_tilt.process(left);
            left += leak_l;
            right += leak_r;
 
            left = channels_[0].print_through(left, print_taps_);
            right = channels_[1].print_through(right, print_taps_);
 
            // Reproduce EQ, then the bias-dependent high shelf cascaded with it,
            // then the expander.
            left = channels_[0].bias_shelf.process(channels_[0].playback_eq.process(left));
            right = channels_[1].bias_shelf.process(channels_[1].playback_eq.process(right));
 
            if (companding_) {
                left = compander_l_.decode(left);
                right = compander_r_.decode(right);
            }
 
            const double output_l = dry_l + mix_ * (left - dry_l);
            const double output_r = dry_r + mix_ * (right - dry_r);
            if (!std::isfinite(output_l) || !std::isfinite(output_r)) {
                recover_audio_fault();
                out_l[n] = SampleType{0};
                out_r[n] = SampleType{0};
                continue;
            }
            out_l[n] = static_cast<SampleType>(output_l);
            out_r[n] = static_cast<SampleType>(output_r);
        }
    }

    int latency_samples() const noexcept {
        return oversampler_latency_samples() + instability_nominal_samples_ +
               (pre_echo_ ? print_taps_.offset_samples : 0);
    }

    double worst_case_insertion_gain() const noexcept {
        double worst = 1.0;
        const double nyquist = 0.5 * sample_rate_;
        for (const Curve curve : {Curve::nab, Curve::iec_ccir, Curve::cassette_type1,
                                  Curve::cassette_type2}) {
            const tape::EqTimeConstants tc = tape::eq_time_constants(curve, speed_ips_);
            TapeEqSection record, playback;
            record.set(tc, sample_rate_, /*inverse=*/false);
            playback.set(tc, sample_rate_, /*inverse=*/true);
            for (double hz = 20.0; hz < nyquist; hz *= 1.005) {
                worst = std::max(worst,
                                 units::db_to_linear(record.response_db(hz, sample_rate_)));
                worst = std::max(worst,
                                 units::db_to_linear(playback.response_db(hz, sample_rate_)));
            }
        }
        return worst * units::db_to_linear(kUnderBiasHfBoostDb);
    }
 
    // ── Test surface ──────────────────────────────────────────────────────
    //
    // The acceptance criteria for the EQ networks, the gap filter and the
    // compander are stated as tolerances on THOSE STAGES (±0.5 dB against an
    // analytic prototype, a ≥40 dB extinction null, ±0.1 dB of round-trip
    // transparency). None of them is measurable end-to-end through a chain that
    // also contains Wallace loss, a head bump, hiss and pitch modulation, so the
    // stages are reachable directly. See the suite's notes.
 
    const TapeEqSection& record_eq(int channel = 0) const noexcept {
        return channels_[static_cast<std::size_t>(channel)].record_eq;
    }
 
    const TapeEqSection& playback_eq(int channel = 0) const noexcept {
        return channels_[static_cast<std::size_t>(channel)].playback_eq;
    }

    const std::vector<double>& gap_fir() const noexcept {
        return channels_[0].gap_taps;
    }

    double reproduce_gap_m() const noexcept {
        return tape::head_gap_geometry(speed_ips_).reproduce_m;
    }

    double reproduce_alignment_db() const noexcept {
        return units::linear_to_db(channels_[0].alignment_gain);
    }
 
    double sample_rate() const noexcept { return sample_rate_; }
 
private:
    void recover_audio_fault() noexcept {
        for (auto& channel : channels_) channel.recover_audio_fault();
        compander_l_.reset();
        compander_r_.reset();
    }

    struct PrintTaps {
        int offset_samples = 0;
        double post_gain = 0.0;
        double pre_gain = 0.0;
        bool pre_echo = false;
    };

    struct Channel {
        chardelay::TapePhysicalChannel physical{};
        chardelay::TapeInstability instability{};
        chardelay::FractionalDelayLine instability_line{};
        chardelay::FixedFir gap_fir{};
        chardelay::FirstOrderShelf bias_shelf{};
        chardelay::FirstOrderShelf crosstalk_tilt{};
        TapeEqSection record_eq{};
        TapeEqSection playback_eq{};
 
        std::vector<double> gap_taps{};
        std::vector<float> print_line{};
        std::size_t print_write = 0;
        std::size_t print_valid = 0;
 
        double wow_depth_ms = 0.0;
        double flutter_depth_ms = 0.0;
        double nominal_samples = 0.0;
        double sample_rate = 48000.0;
        double alignment_gain = 1.0;

        static double fir_magnitude(const std::vector<double>& taps, double hz,
                                    double fs) noexcept {
            if (taps.empty()) return 1.0;
            double re = 0.0, im = 0.0;
            for (std::size_t n = 0; n < taps.size(); ++n) {
                const double theta = 2.0 * chardelay::kPi * hz * static_cast<double>(n) / fs;
                re += taps[n] * std::cos(theta);
                im -= taps[n] * std::sin(theta);
            }
            return std::hypot(re, im);
        }

        void align(double ceiling_db) noexcept {
            const double tilt_db = chardelay::tape_loss_iir_magnitude_db(
                physical.loss_parameters(), tape::kEqReferenceHz);
            const double reused_gap =
                fir_magnitude(physical.gap_coefficients(), tape::kEqReferenceHz, sample_rate);
            const double own_gap = fir_magnitude(gap_taps, tape::kEqReferenceHz, sample_rate);
            const double loss = units::db_to_linear(tilt_db) *
                                std::max(reused_gap, 1e-6) * std::max(own_gap, 1e-6);
            alignment_gain = std::min(1.0 / std::max(loss, 1e-6),
                                      units::db_to_linear(ceiling_db));
        }
 
        void prepare(double fs, int index) {
            sample_rate = fs;
            physical.prepare(fs);
            instability.prepare(fs);
            // Distinct seeds per channel, DERIVED from the shared constant
            // rather than drawn, so the two channels decorrelate while both
            // renders stay bit-identical across a reset (series law 2). The
            // right channel also starts its flutter stack at the reused
            // fixed phase offset.
            const auto seed = static_cast<std::uint32_t>(
                chardelay::kPrngSeed ^ (0x9E3779B9u * static_cast<std::uint32_t>(index)));
            instability.configure(seed, index == 0 ? 0.0 : chardelay::kFlutterPhaseRight);
            physical.set_seeds(seed);
 
            nominal_samples = units::ms_to_samples(kInstabilityNominalMs, fs);
            const auto capacity = static_cast<std::size_t>(
                std::llround(nominal_samples + units::ms_to_samples(4.0, fs))) + 8u;
            instability_line.prepare(capacity);
 
            const auto print_capacity = static_cast<std::size_t>(std::llround(
                2.0 * units::ms_to_samples(kPrintThroughOffsetMsMax, fs))) + 4u;
            print_line.assign(print_capacity, 0.0f);
            print_write = 0;
            print_valid = 0;
        }
 
        void reset() {
            physical.reset();
            instability.reset();
            instability_line.reset();
            gap_fir.reset();
            bias_shelf.reset();
            crosstalk_tilt.reset();
            record_eq.reset();
            playback_eq.reset();
            std::fill(print_line.begin(), print_line.end(), 0.0f);
            print_write = 0;
            print_valid = 0;
        }
 
        void recover_audio_fault() noexcept {
            physical.recover_audio_fault();
            instability.reset();
            instability_line.discard_history();
            gap_fir.discard_history();
            bias_shelf.reset();
            crosstalk_tilt.reset();
            record_eq.reset();
            playback_eq.reset();
            print_write = 0;
            print_valid = 0;
        }
 
        void set_gap_taps(std::vector<double> taps) {
            gap_taps = std::move(taps);
            gap_fir.prepare(gap_taps.size());
        }

        double record(double x, double drive_gain, double inverse_drive_gain) noexcept {
            const double emphasized = record_eq.process(x);
            return physical.pre_process(emphasized * drive_gain) * inverse_drive_gain;
        }

        double playback(double x) noexcept {
            gap_fir.push(x);
            const double gapped = gap_taps.empty()
                                      ? x
                                      : gap_fir.convolve(gap_taps.data(), gap_taps.size());
            const double lossy = physical.post_process(gapped) * alignment_gain;
 
            instability.tick();
            instability_line.push(lossy);
            const double offset_ms = instability.offset_ms(wow_depth_ms, flutter_depth_ms);
            const double read =
                nominal_samples + units::ms_to_samples(offset_ms, sample_rate);
            return instability_line.read(read);
        }

        double print_through(double x, const PrintTaps& taps) noexcept {
            const std::size_t size = print_line.size();
            if (size == 0 || taps.offset_samples <= 0) return x;
            const auto offset = static_cast<std::size_t>(taps.offset_samples);
 
            print_line[print_write] = static_cast<float>(x);
            const auto at = [&](std::size_t back) {
                if (back >= print_valid + 1u) return 0.0;
                return static_cast<double>(print_line[(print_write + size - back) % size]);
            };
 
            double y = 0.0;
            if (taps.pre_echo) {
                // Programme delayed by one offset; the print arrives at the head
                // of the stream (early) and again two offsets back (late).
                y = at(offset) + taps.pre_gain * at(0) + taps.post_gain * at(2u * offset);
            } else {
                y = at(0) + taps.post_gain * at(offset);
            }
            print_write = (print_write + 1u) % size;
            if (print_valid < size) ++print_valid;
            return y;
        }
    };
 
    void update_speed(double ips, bool redesign) {
        const bool changed = std::abs(ips - speed_ips_) > 1e-9;
        speed_ips_ = ips;
        if (!changed && !redesign) return;
 
        const tape::HeadGapGeometry geometry = tape::head_gap_geometry(speed_ips_);
        chardelay::TapeLossGeometry gap_geometry;
        gap_geometry.speed_ips = speed_ips_;
        gap_geometry.gap_m = geometry.reproduce_m;
        // The reused minimum-phase design helper, given THIS module's speed-class
        // reproduce gap rather than the delay tier's fixed 3 µm head. Same
        // construction (cepstral fold), a different filter.
        std::vector<double> taps = chardelay::design_tape_gap_fir(
            sample_rate_, chardelay::tape_gap_fir_taps(sample_rate_), gap_geometry);
 
        for (auto& channel : channels_) {
            channel.set_gap_taps(taps);
            channel.physical.set_speed_ips(speed_ips_);
            channel.align(kAlignmentCeilingDb);
        }
    }
 
    void update_eq() {
        const tape::EqTimeConstants tc = tape::eq_time_constants(curve_, speed_ips_);
        for (auto& channel : channels_) {
            channel.record_eq.set(tc, sample_rate_, /*inverse=*/false);
            channel.playback_eq.set(tc, sample_rate_, /*inverse=*/true);
        }
    }

    void apply_age() {
        for (auto& channel : channels_) {
            channel.physical.update(age01_);
            // Age moves the reused tier's spacing, so realignment belongs here
            // and at every speed change — the same two events that move it on a
            // real machine.
            channel.align(kAlignmentCeilingDb);
        }
        print_through_db_ = tape::age_print_through_db(age01_);
        update_bias();
    }
 
    void update_bias() {
        // The age table's drift is ADDED to what the user asked for, then
        // clamped — a tired machine biases hotter than its front panel says.
        bias_ = std::clamp(bias_request_ + tape::age_bias_drift(age01_), kBiasMin, kBiasMax);
 
        const double hf_shelf_db = -kOverBiasHfCutDb * std::max(0.0, bias_) +
                                   kUnderBiasHfBoostDb * std::max(0.0, -bias_);
        const double shelf_gain = units::db_to_linear(hf_shelf_db);
        for (auto& channel : channels_)
            channel.bias_shelf.set(kBiasShelfHz, shelf_gain, sample_rate_);
 
        const double distortion_mult = 1.0 + kBiasDistK * bias_ * bias_;
        drive_gain_ = units::db_to_linear(drive01_ * kDriveDbMax * distortion_mult);
        inverse_drive_gain_ = 1.0 / drive_gain_;
 
        // Wow and flutter depths come from the reused table, driven by this
        // module's own age axis.
        const double wow = chardelay::interpolate_knots(chardelay::kTapeAxis,
                                                        chardelay::kTapeWowDepthMs, age01_);
        const double flutter = chardelay::interpolate_knots(
            chardelay::kTapeAxis, chardelay::kTapeFlutterDepthMs, age01_);
        for (auto& channel : channels_) {
            channel.wow_depth_ms = wow;
            channel.flutter_depth_ms = flutter;
        }
    }
 
    void update_crosstalk() {
        const double tilt = units::db_to_linear(kCrosstalkTiltMaxDb);
        for (auto& channel : channels_)
            channel.crosstalk_tilt.set(kCrosstalkTiltHz, tilt, sample_rate_);
        // The tilt shelf is normalized to unity at the same 1 kHz reference the
        // EQ networks use, so a crosstalk measurement at 1 kHz reads the
        // configured figure EXACTLY rather than the figure plus half a decibel
        // of shelf skirt. Normalizing against the realized digital shelf, not
        // its analog prototype, so the bilinear's frequency warping does not
        // leak into that measurement. The stated ceiling becomes
        // `crosstalk_db + kCrosstalkTiltMaxDb − (skirt at 1 kHz)`, which is
        // strictly below the ceiling the spec asks for.
        crosstalk_gain_ = units::db_to_linear(
            crosstalk_db_ - tape::first_order_shelf_db(kCrosstalkTiltHz, tilt,
                                                       tape::kEqReferenceHz, sample_rate_));
    }
 
    void update_print_through() {
        print_taps_.offset_samples = std::min(
            static_cast<int>(std::llround(units::ms_to_samples(print_offset_ms_, sample_rate_))),
            print_max_samples_);
        const double reduction = kBiasErasureDb * std::max(0.0, bias_);
        const double post_db = print_through_db_ - reduction;
        print_taps_.post_gain = units::db_to_linear(post_db);
        print_taps_.pre_gain = units::db_to_linear(post_db - kPreEchoAsymmetryDb);
        print_taps_.pre_echo = pre_echo_;
    }
 
    std::array<Channel, 2> channels_{};
    TapeCompanderT<SampleType> compander_l_{};
    TapeCompanderT<SampleType> compander_r_{};
 
    double sample_rate_ = 48000.0;
    Archetype archetype_ = Archetype::ampex_350_440;
    Curve curve_ = Curve::nab;
    double speed_ips_ = 15.0;
    double bias_request_ = kBiasDefault;
    double bias_ = kBiasDefault;
    double drive01_ = kDriveDefault;
    double age01_ = 0.0;
    double crosstalk_db_ = -38.0;
    bool companding_ = false;
    double print_through_db_ = kPrintThroughDbMin;
    double print_offset_ms_ = kPrintThroughOffsetMsDefault;
    bool pre_echo_ = false;
    double mix_ = 1.0;
 
    double drive_gain_ = 1.0;
    double inverse_drive_gain_ = 1.0;
    double crosstalk_gain_ = 0.0;
    PrintTaps print_taps_{};
    int instability_nominal_samples_ = 0;
    int print_max_samples_ = 0;
};
 
using TapeMachine = TapeMachineT<float>;
using TapeMachine64 = TapeMachineT<double>;
using TapeCompander = TapeCompanderT<float>;
 
}  // namespace pulp::signal