Retro 80s Tape Wow & Flutter with faust

[Morphoice]

@HISEnberg brilliant thanks. I had already considered to add some randomness something like

rate = ma.SR/10000.0; // new random value every 100 samples (ma.SR from maths.lib)
process = no.lfnoiseN(3,rate)*0.5 + no.lfnoiseN(3,400)*0.1;

instead of a fixed wow frequency. (this just being a test to create a low frequency random noise)
also trying to figure out some ways to "degrade" the signal. I've indeed studied Jatin's work on the tape model but it's giving me a hard time to adapt to faust

[HISEnberg]

@Morphoice Same here, it is insightful but I am still making my way through it.

I suppose for degrading, you could consider frequency-dependent saturation? I think tape saturation ought to effect low frequencies more than high ones.
Other than that, a slight compression at the end.
I also considered something like an "Tape Age" slider which would represent newer to older tape. This could probably be implemented with another noise source and LPF, but some type of saturation would be nice too. Also it would cause greater wow and flutter, but this get's a little unwieldy so I think the values need to be readjusted.

Another thing I did was introduce a bit of stereo processing, so the modulations on the left and right were somewhat offset. I find it helps create a better stereo field but this isn't necessarily something tape does.

import ("stdfaust.lib");
import ("analyzers.lib");

// Tape Age Control
tape_age = hslider("Tape Age", 0.0, 0.0, 1.0, 0.01);  // 0 = New, 1 = Old

// Wow and Flutter
random_mod = hslider("random_mod", 0.002, 0, 0.01, 0.0001);

wow = (os.lf_trianglepos(hslider("wow", 0.5, 0, 2, 0.01)) 
      + no.noise * random_mod 
      + os.oscsin(0.3) * 0.0005) 
      : si.smooth(0.9999) * 1000.;

flutter = ((os.lf_trianglepos(hslider("flutter", 8, 2, 20, 0.01)) 
          + no.noise * random_mod) * (0.8 + os.oscsin(0.1) * 0.2) 
          + os.oscsin(15) * 0.0003) 
          : si.smooth(0.9999) * 1000.;

wow_intensity = hslider("Wow Intensity", 0.3, 0, 1, 0.01) * (1 + tape_age * 0.5) : si.smooth(0.9999);  // Scaled by tape age
flutter_intensity = hslider("Flutter Intensity", 0.1, 0, 1, 0.01) * (1 + tape_age * 0.5) : si.smooth(0.9999);  // Scaled by tape age

// Stereo Spread
stereo = hslider("Stereo Spread", 0.1, 0, 1, 0.01);

wow_left = wow;
wow_right = wow + no.noise * (stereo * 0.001);  // Subtle difference
flutter_left = flutter;
flutter_right = flutter + no.noise * (stereo);  // Subtle difference

// Noise
hiss = hslider("Hiss", 0.1, 0, 1, 0.01);

tapenoise = (no.pink_noise * (0.6 + os.lf_triangle(0.1) * 0.05) 
             + no.sparse_noise(10) * (0.2 + os.lf_triangle(0.05) * 0.05) 
             + os.oscsin(50) * 0.01) 
             : fi.highpass(1, 50) 
             : fi.lowpass(1, 15000) 
             / 5 
             + os.triangle(60) * 0.05 
             + os.oscsin(528) * 0.003 
             + os.oscsin(110) * 0.002;

// Apply Tape Age to Noise
aged_tapenoise = tapenoise : fi.lowpass(1, 15000 - tape_age * 8000) * (1 + tape_age * 0.5);  // Aged noise

// Saturator
drive = hslider("Drive", 0, 0, 1, 0.001) + tape_age * 0.2;  // Increased drive with tape age
tanh(x) = x * (27 + x * x) / (27 + 9 * x * x);
saturator = ef.dryWetMixerConstantPower(drive, tanh * drive * 20 : co.limiter_1176_R4_mono : tanh);

// Combine Wow, Flutter, and Noise
wowflutterLeft = de.fdelay1(ma.SR, 100. + wow_left * wow_intensity + flutter_left * (flutter_intensity / 5)) 
                  <: saturator 
                  : fi.lowpass(2, 8000) 
                  + aged_tapenoise * hiss;

wowflutterRight = de.fdelay1(ma.SR, 100. + wow_right * wow_intensity + flutter_right * (flutter_intensity / 5)) 
                   <: saturator 
                   : fi.lowpass(2, 8000) 
                   + aged_tapenoise * hiss;

// Output
process = wowflutterLeft, wowflutterRight;

[Morphoice]

@HISEnberg I've thought about stereo processing, but it would imply having two mono tape machines... a single stereo tape would always run the same speed and flutter/vibrate the same as it's a single stereo head picking up the signal...
again though this would offer a bit of variance and introduce an abstract option to go more into experimental sounds

as for the degrading it's really the hysteresis I want to tackle, as all the normal saturation circuits we use are just a "cheat"... also having some transformer things like the Kush Omega Transformer

[HISEnberg]

@Morphoice yes I thought of that too regarding the stereo processing. I suppose it would be useful if you want to retain the original, stereo signal input. Otherwise it's just for experimentation.

[Lindon]

@HISEnberg @Morphoice - as I mentioned before if your seriously thinking of emulating tape machines you will need to account for head bump - and probably make it variable as different machines from different manufacturers had different amounts in different freq...

[Morphoice]

@Lindon good point. I could add some eq to get close to an average tape response curve. Not really looking to model specific tape machines, just want to have that wobbly 80s retro sound much like baby audio's VHS

[griffinboy]

@Morphoice

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.