Retro 80s Tape Wow & Flutter with faust

griffinboy

I've been doing some work to extract the real transport from multiple machines. I'm not terribly convinced by the quasi-random shapes most people use to model wow.

In terms of hysteresis and the frequency response of tape, everyone else here is correct: Jatin's work (chow tape) is a good place to start.

HISEnberg

@Morphoice That reflects my use case as well, I am more interested in experimenting with (and pushing) the potentials of tape styled-modulation.

I agree with you that the saturation element is queue here, I will let you know if I make any headway with the hysteresis (though I remember reading its full implementation is super CPU intensive).

I suppose a nifty feature to explore could be tape startup, rewind, and different tape speeds.

Ben Catman

Good work !

Morphoice

@griffinboy I agree. As the wow results from the revolving reel or cassette it's rather periodic than random, the random wow makes more sense to me on a tape delay which has a bunch of lose tapeloop in a cartridge

Morphoice

@HISEnberg sadly I'm far from being able to implement the hysteresis myself or port it to faust somehow, but if you ever manage to do so please let me know. the question then again is, though, does the end-user even hear the difference of is a simple saturator sufficient. I'm not even sure all the big companies like Waves tape emulation plugins even have something like exact hysteresis. Jatin points out the various solvers to calculate a hysteresis loop and they are indeed quite cpu heavy. maybe the UAD plugins did such on thing on their DSP

griffinboy

@Morphoice

If you want a c++ node for it (works in scriptnode) I can lend you mine.
3% cpu usage when you oversample it, if I remember correctly. It's a clone of Chow tape's hysteresis, with some minor alterations and no dependancies.

Morphoice

@griffinboy I'd love to take a look at it, though I've never used a C++ node before and will probably fail terribly at getting it running lol

griffinboy

@Morphoice

Nope, it's incredibly easy to use a c++ node!
This one is stereo only though, not multichannel but apart from that it's stable.
It requires a clipper to be placed before it. The simulation is unstable when driven very hard. To find the right place to clip, use a hise limit node before the hysteresis, and also a gain node before that. Now drive the gain node into the hysteresis at somepoint it will collapse. Lower the limiter until it no longer collapses.

To import the c++ node, create this file named JG_Tape_Model.h

Here is the file

#pragma once
#include <JuceHeader.h>
#include <cmath>
#include <algorithm>
#include <cassert>


namespace project
{
    using namespace juce;
    using namespace hise;
    using namespace scriptnode;

    // ==========================| Constants |==========================
    constexpr double ONE_THIRD = 1.0 / 3.0;
    constexpr double alpha = 1.6e-3;
    double upperLim = 20.0;

    // ==========================| Utility Functions |==========================
    inline int sign(double x) {
        return int(x > 0.0) - int(x < 0.0);
    }

    inline double tanh_approx(double x) {
        if (x < -3.0) return -1.0;
        if (x > 3.0) return 1.0;
        double x2 = x * x;
        return x * (135135.0 + x2 * (17325.0 + x2 * (378.0 + x2))) /
            (135135.0 + x2 * (62370.0 + x2 * (3150.0 + x2 * 28.0)));
    }

    inline double coth_approx(double x) {
        double tanh_x = tanh_approx(x);
        return 1.0 / tanh_x;
    }


    // ==========================| HysteresisState Struct |==========================
    struct HysteresisState {
        double M_s = 1.0;
        double a = M_s / 4.0;
        double k = 0.47875;
        double c = 1.7e-1;
        double nc = 1 - c;
        double M_s_oa = M_s / a;
        double M_s_oa_talpha = alpha * M_s / a;
        double M_s_oa_tc = c * M_s / a;
        double M_s_oa_tc_talpha = alpha * c * M_s / a;
        double M_s_oaSq_tc_talpha = alpha * c * M_s / (a * a);
        double M_s_oaSq_tc_talphaSq = alpha * alpha * c * M_s / (a * a);

        double Q, M_diff, L_prime, kap1, f1Denom, f1, f2, f3;
        double coth = 0.0;
        bool nearZero = false;
        double oneOverQ, oneOverQSq, oneOverQCubed, cothSq, oneOverF3, oneOverF1Denom;
    };

    // ==========================| Langevin Functions |==========================
    template <typename Float>
    Float langevin(const HysteresisState& hp) noexcept {
        constexpr double threshold = 1e-4;
        return std::abs(hp.Q) > threshold ? (hp.coth) - (hp.oneOverQ) : hp.Q * ONE_THIRD;
    }

    template <typename Float>
    Float langevinD(const HysteresisState& hp) noexcept {
        constexpr double threshold = 1e-4;
        return std::abs(hp.Q) > threshold ? hp.oneOverQSq - hp.cothSq + 1.0 : ONE_THIRD;
    }

    // ==========================| Hysteresis Function |==========================
    template <typename Float>
    Float hysteresisFunc(double M, double H, double H_d, HysteresisState& hp) noexcept {
        hp.Q = H - hp.k * M;
        hp.oneOverQ = 1.0 / hp.Q;
        hp.oneOverQSq = hp.oneOverQ * hp.oneOverQ;
        hp.oneOverQCubed = hp.oneOverQ * hp.oneOverQSq;

        // Use the approximations
        hp.coth = coth_approx(hp.Q);
        constexpr double threshold = 1e-4;
        hp.nearZero = std::abs(hp.Q) < threshold;

        hp.cothSq = hp.coth * hp.coth;
        hp.M_diff = langevin<Float>(hp) * hp.M_s - M;

        const auto delta = (H_d >= 0.0) - (H_d < 0.0);
        const auto delta_M = sign(delta) == sign(hp.M_diff);
        hp.kap1 = hp.nc * delta_M;

        hp.L_prime = langevinD<Float>(hp);

        hp.f1Denom = (hp.nc * delta) * hp.k - alpha * hp.M_diff;
        hp.oneOverF1Denom = 1.0 / hp.f1Denom;
        hp.f1 = hp.kap1 * hp.M_diff * hp.oneOverF1Denom; // Use cached value
        hp.f2 = hp.L_prime * hp.M_s_oa_tc;
        hp.f3 = 1.0 - (hp.L_prime * hp.M_s_oa_tc_talpha);
        hp.oneOverF3 = 1.0 / hp.f3;

        return H_d * (hp.f1 + hp.f2) * hp.oneOverF3;
    }

    // ==========================| HysteresisProcessing Class |==========================
    class HysteresisProcessing {
    public:
        HysteresisProcessing() {
            reset();
        }

        void reset() {
            M_n1 = 0.0;
            H_n1 = 0.0;
            H_d_n1 = 0.0;
            hpState.coth = 0.0;
            hpState.nearZero = false;
        }

        void softReset() {
            M_n1 = 0.0;
            H_n1 = 0.0;
            H_d_n1 = 0.0;
            hpState = HysteresisState(); // Reset to default values
        }

        void setSampleRate(double newSR) {
            fs = newSR;
            T = 1.0 / fs;
        }

        void configure(double drive, double width, double sat, bool v1) {
            hpState.M_s = 0.5 + 1.5 * (1.0 - sat);
            hpState.a = hpState.M_s / (0.01 + 6.0 * drive);
            hpState.c = std::sqrt(1.0f - width) - 0.01;
            hpState.k = 0.47875;

            if (v1) {
                hpState.k = 27.0e3;
                hpState.c = 1.7e-1;
                hpState.M_s *= 50000.0;
                hpState.a = hpState.M_s / (0.01 + 40.0 * drive);
            }

            hpState.nc = 1.0 - hpState.c;
            hpState.M_s_oa = hpState.M_s / hpState.a;
            hpState.M_s_oa_talpha = alpha * hpState.M_s_oa;
            hpState.M_s_oa_tc = hpState.c * hpState.M_s_oa;
            hpState.M_s_oa_tc_talpha = alpha * hpState.M_s_oa_tc;
            hpState.M_s_oaSq_tc_talpha = hpState.M_s_oa_tc_talpha / hpState.a;
            hpState.M_s_oaSq_tc_talphaSq = alpha * hpState.M_s_oaSq_tc_talpha;
        }

        template <typename Float>
        Float process(Float H) noexcept {
            // Limit input to prevent extreme values
            H = std::clamp(H, static_cast<Float>(-upperLim), static_cast<Float>(upperLim));

            // Initialize H_d before using it
            Float H_d = deriv<Float>(H, H_n1, H_d_n1, static_cast<Float>(T));

            // RK4 solver
            Float k1 = T * hysteresisFunc<Float>(M_n1, H, H_d, hpState);
            Float k2 = T * hysteresisFunc<Float>(M_n1 + 0.5 * k1, H, H_d, hpState);
            Float k3 = T * hysteresisFunc<Float>(M_n1 + 0.5 * k2, H, H_d, hpState);
            Float k4 = T * hysteresisFunc<Float>(M_n1 + k3, H, H_d, hpState);

            Float M = M_n1 + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;

            bool illCondition = std::isnan(M) || std::isinf(M) || M > upperLim || M < -upperLim;
            if (illCondition) {
                softReset();
                return 0.0;
            }

            // Add a limiter to prevent extreme values
            M = std::clamp(M, static_cast<Float>(-upperLim), static_cast<Float>(upperLim));

            M_n1 = M;
            H_n1 = H;
            H_d_n1 = H_d;

            return M;
        }

        // Ensure the deriv function is correctly defined and accessible
        template <typename Float>
        static Float deriv(Float H, Float H_n1, Float H_d_n1, Float T) noexcept {
            return (H - H_n1) / T;
        }

    private:
        double fs = 48000.0;
        double T = 1.0 / fs;
        double M_n1 = 0.0;
        double H_n1 = 0.0;
        double H_d_n1 = 0.0;
        HysteresisState hpState;
    };

    // ==========================| JG_Tape_Model Node Class |==========================
    template <int NV> struct JG_Tape_Model : public data::base
    {
        // ==========================| Metadata Definitions |==========================
        SNEX_NODE(JG_Tape_Model);

        struct MetadataClass
        {
            SN_NODE_ID("JG_Tape_Model");
        };

        static constexpr bool isModNode() { return false; }
        static constexpr bool isPolyphonic() { return NV > 1; }
        static constexpr bool hasTail() { return false; }
        static constexpr bool isSuspendedOnSilence() { return false; }
        static constexpr int getFixChannelAmount() { return 2; }

        static constexpr int NumTables = 0;
        static constexpr int NumSliderPacks = 0;
        static constexpr int NumAudioFiles = 0;
        static constexpr int NumFilters = 0;
        static constexpr int NumDisplayBuffers = 0;

        // ==========================| External Variables |==========================
        float Drive = 0.2f;
        float Width = 0.5f;
        float Saturation = 0.1f;
        bool V1 = false;

        // ==========================| Internal Variables |==========================
        float sampleRate = 48000.0f;
        ModValue modValue;
        HysteresisProcessing hProc;

        // ==========================| Helper Functions |==========================

        // ==========================| Preparation and Reset |==========================
        void prepare(PrepareSpecs prepareSpecs)
        {
            sampleRate = prepareSpecs.sampleRate;
            hProc.setSampleRate(sampleRate);
            hProc.configure(Drive, Width, Saturation, V1); // Example parameters
        }

        void reset()
        {
            hProc.reset();
        }

        // ==========================| Processing |==========================
        // Non Interleaved Sample Processing
        template <typename T>
        void process(T& data)
        {
            auto& fixData = data.template as<ProcessData<getFixChannelAmount()>>();
            auto fd = fixData.toFrameData();

            while (fd.next())
            {
                processFrame(fd.toSpan());
            }
        }

        // Interleaved Sample Processing
        template <typename T> void processFrame(T& data)
        {
            const size_t numChannels = 2; // We only have two channels
            const size_t numSamples = data.size() / numChannels;

            // Process each channel differently
            for (size_t channel = 0; channel < numChannels; ++channel)
            {
                for (size_t i = 0; i < numSamples; ++i)
                {
                    size_t index = i * numChannels + channel;

                    if (channel == 0)
                    {
                        // Process left channel
                        data[index] = hProc.process(data[index]);
                    }
                    else if (channel == 1)
                    {
                        // Process right channel
                        data[index] = hProc.process(data[index]);
                    }
                }
            }
        }

        // ==========================| Modulation and Parameters |==========================
        int handleModulation(double& value)
        {
            return modValue.getChangedValue(value);
        }

        void setExternalData(const ExternalData& data, int index) {}
        void handleHiseEvent(HiseEvent& e) {}

        template <int P>
        void setParameter(double v)
        {
            if (P == 0)
            {
                Drive = static_cast<float>(v);
                hProc.configure(Drive, Width, Saturation, V1); // Update processing
            }
            else if (P == 1)
            {
                Width = static_cast<float>(v);
                hProc.configure(Drive, Width, Saturation, V1); // Update processing
            }
            else if (P == 2)
            {
                Saturation = static_cast<float>(v);
                hProc.configure(Drive, Width, Saturation, V1); // Update processing
            }
            else if (P == 3)
            {
                V1 = static_cast<bool>(v);
                hProc.configure(Drive, Width, Saturation, V1); // Update processing
            }
            else if (P == 4)
            {
                upperLim = static_cast<bool>(v);
                hProc.configure(Drive, Width, Saturation, V1); // Update processing
            }
        }

        void createParameters(ParameterDataList& data)
        {
            {
                parameter::data p("Drive", { 0.01, 1.0 });
                registerCallback<0>(p);
                p.setDefaultValue(0.2);
                data.add(std::move(p));
            }
            {
                parameter::data p("Width", { 0.0, 1.0 });
                registerCallback<1>(p);
                p.setDefaultValue(0.5);
                data.add(std::move(p));
            }
            {
                parameter::data p("Saturation", { 0.0, 1.0 });
                registerCallback<2>(p);
                p.setDefaultValue(0.1);
                data.add(std::move(p));
            }
            {
                parameter::data p("V1", { 0.0, 1.0 });
                registerCallback<3>(p);
                p.setDefaultValue(0.0);
                data.add(std::move(p));
            }
            {
                parameter::data p("Limit", { 0.0, 30.0 });
                registerCallback<4>(p);
                p.setDefaultValue(20.0);
                data.add(std::move(p));
            }
        }
    };
}

Place it in your hise project, under DspNetworks > Thirdparty

Then open the project in hise and select Export > Compile dsp Networks as dll

griffinboy

@griffinboy

I uploaded the wrong file a second ago, so if you imported it, I've fixed it now lol

Morphoice

@griffinboy alright let me walk through this real quick, I guarantee I'll explode the mac in a minute ;))) everything C++ I tough is doomed but I'll try my best

Morphoice

@griffinboy alright, I complied the DSP network, how do I access the node inside my DSP network now? I can only see it as a hardcoded fx

griffinboy

@Morphoice

In scriptnode open a new node, and go under the 'projects' tab category of nodes

Morphoice

@griffinboy brilliant!!! I'll try the clipping thing now

griffinboy

@Morphoice

set the gain node to like 0 to 60 range. Playing audio through the high end of that will cause the output to be silent because the model has collapsed. To reset this you have to do a reset, this can be done using the hise keyboard there is a button with a circle and an exclamation mark

also wrap the whole chain in oversample 2x or 4x.
2x should be fine. The more oversampling the less aliasing, and the more volume you can boost before a collapse.

Morphoice

@griffinboy I used a pretty hard mastered cyberpunk track as test audio, I can't get it to collapse and go silent, but the drive drops out at certain hard points, and it just becomes quieter, which is what I suppose is the collapse... upon turning the drive button it comes back and stays if I limit enough

this is what it looks
Screenshot 2024-11-27 at 23.06.30.jpg

Morphoice

@griffinboy it does become mono btw when drive is below around 0.1, what would be a good range for drive/saturation and width to offer the user? clearly the extreme values are a nono

griffinboy

@Morphoice

hmm it's not done that on my machine.
I wonder if I've given you a bad version, I've done a lot of tinkering with that script over time.

The default settings should be left as is! Ignore all parameters on the actual unit apart from saturation which can be adjusted

griffinboy

@Morphoice

Ah sorry not limiter node.

I meant Math.clip node

griffinboy

@Morphoice

To find the right place to clip, use a hise Math.clip node before the hysteresis, and a simple gain node before that. Now drive the gain node into the hysteresis at somepoint it will collapse. Lower the Clip node until it no longer collapses.

Node graph:

Simple Gain

Math.clip

Oversample 4x - JG_Tape_Model -

send a sine wave into it and visualise using fft

Morphoice

@griffinboy there's also a significant gain reduction even without limiting, which I guess is fine, I wonder if it would make more sense to implement the JG tape model with a dry/wet knob or give the user actual control over the drive/saturation

Retro 80s Tape Wow & Flutter with faust

18

1.9k

10.5k

91.7k