AnalogTapeModel

Physical modelling signal processing for analog tape recording
Log | Files | Refs | Submodules | README | LICENSE

HysteresisProcessor.cpp (14397B)


      1 #include "HysteresisProcessor.h"
      2 
      3 namespace
      4 {
      5 enum
      6 {
      7     numSteps = 500,
      8 };
      9 
     10 constexpr double v1Norm = 1.414 / 10000.0;
     11 
     12 template <typename T>
     13 static void interleaveSamples (const T** source, T* dest, int numSamples, int numChannels)
     14 {
     15     for (int chan = 0; chan < numChannels; ++chan)
     16     {
     17         auto i = chan;
     18         const auto* src = source[chan];
     19 
     20         for (int j = 0; j < numSamples; ++j)
     21         {
     22             dest[i] = src[j];
     23             i += numChannels;
     24         }
     25     }
     26 }
     27 
     28 template <typename T>
     29 static void deinterleaveSamples (const T* source, T** dest, int numSamples, int numChannels)
     30 {
     31     for (int chan = 0; chan < numChannels; ++chan)
     32     {
     33         auto i = chan;
     34         auto* dst = dest[chan];
     35 
     36         for (int j = 0; j < numSamples; ++j)
     37         {
     38             dst[j] = source[i];
     39             i += numChannels;
     40         }
     41     }
     42 }
     43 } // namespace
     44 
     45 HysteresisProcessor::HysteresisProcessor (AudioProcessorValueTreeState& vts) : osManager (vts)
     46 {
     47     using namespace chowdsp::ParamUtils;
     48     loadParameterPointer (driveParam, vts, "drive");
     49     loadParameterPointer (satParam, vts, "sat");
     50     loadParameterPointer (widthParam, vts, "width");
     51     modeParam = vts.getRawParameterValue ("mode");
     52     onOffParam = vts.getRawParameterValue ("hyst_onoff");
     53 }
     54 
     55 void HysteresisProcessor::createParameterLayout (chowdsp::Parameters& params)
     56 {
     57     using namespace chowdsp::ParamUtils;
     58     emplace_param<chowdsp::BoolParameter> (params, "hyst_onoff", "Tape On/Off", true);
     59     emplace_param<chowdsp::FloatParameter> (params, "drive", "Tape Drive", NormalisableRange { 0.0f, 1.0f }, 0.5f, &floatValToString, &stringToFloatVal);
     60     emplace_param<chowdsp::FloatParameter> (params, "sat", "Tape Saturation", NormalisableRange { 0.0f, 1.0f }, 0.5f, &floatValToString, &stringToFloatVal);
     61     emplace_param<chowdsp::FloatParameter> (params, "width", "Tape Bias", NormalisableRange { 0.0f, 1.0f }, 0.5f, &floatValToString, &stringToFloatVal);
     62 
     63     emplace_param<chowdsp::ChoiceParameter> (params, "mode", "Tape Mode", StringArray ({ "RK2", "RK4", "NR4", "NR8", "STN", "V1" }), 0);
     64 
     65     using OSManager = decltype (osManager);
     66     OSManager::createParameterLayout (params, OSManager::OSFactor::TwoX, OSManager::OSMode::MinPhase);
     67 }
     68 
     69 void HysteresisProcessor::setSolver (int newSolver)
     70 {
     71     // Hack for V1 solver mode
     72     useV1 = newSolver == SolverType::NUM_SOLVERS;
     73     solver = useV1 ? RK4 : static_cast<SolverType> (newSolver);
     74 
     75     // set clip level for solver
     76     switch (solver)
     77     {
     78         case RK2:
     79             clipLevel = 8.0f;
     80             return;
     81 
     82         case RK4:
     83             clipLevel = 10.0f;
     84             return;
     85 
     86         case NR4:
     87         case NR8:
     88             clipLevel = 12.5f;
     89             return;
     90 
     91         default:
     92             clipLevel = 20.0f;
     93     };
     94 }
     95 
     96 double HysteresisProcessor::calcMakeup()
     97 {
     98     return (1.0 + 0.6 * width[0].getTargetValue()) / (0.5 + 1.5 * (1.0 - sat[0].getTargetValue()));
     99 }
    100 
    101 void HysteresisProcessor::setDrive (float newDrive)
    102 {
    103     for (auto& driveVal : drive)
    104         driveVal.setTargetValue ((double) newDrive);
    105 }
    106 
    107 void HysteresisProcessor::setWidth (float newWidth)
    108 {
    109     for (auto& widthVal : width)
    110         widthVal.setTargetValue ((double) newWidth);
    111 }
    112 
    113 void HysteresisProcessor::setSaturation (float newSaturation)
    114 {
    115     for (auto& satVal : sat)
    116         satVal.setTargetValue ((double) newSaturation);
    117 }
    118 
    119 void HysteresisProcessor::setOversampling()
    120 {
    121     if (osManager.updateOSFactor())
    122     {
    123         for (size_t ch = 0; ch < hProcs.size(); ++ch)
    124         {
    125             hProcs[ch].setSampleRate (fs * osManager.getOSFactor());
    126             hProcs[ch].cook (drive[ch].getCurrentValue(), width[ch].getCurrentValue(), sat[ch].getCurrentValue(), wasV1);
    127             hProcs[ch].reset();
    128         }
    129 
    130         calcBiasFreq();
    131     }
    132 }
    133 
    134 void HysteresisProcessor::calcBiasFreq()
    135 {
    136     biasFreq = fs * osManager.getOSFactor() / 2.0;
    137 }
    138 
    139 void HysteresisProcessor::prepareToPlay (double sampleRate, int samplesPerBlock, int numChannels)
    140 {
    141     fs = sampleRate;
    142     wasV1 = useV1;
    143 
    144     osManager.prepareToPlay (sampleRate, samplesPerBlock, numChannels);
    145     calcBiasFreq();
    146 
    147     drive.resize ((size_t) numChannels);
    148     for (auto& val : drive)
    149         val.reset (numSteps);
    150 
    151     width.resize ((size_t) numChannels);
    152     for (auto& val : width)
    153         val.reset (numSteps);
    154 
    155     sat.resize ((size_t) numChannels);
    156     for (auto& val : sat)
    157         val.reset (numSteps);
    158 
    159     hProcs.resize ((size_t) numChannels);
    160     for (size_t ch = 0; ch < (size_t) numChannels; ++ch)
    161     {
    162         hProcs[ch].setSampleRate (sampleRate * osManager.getOSFactor());
    163         hProcs[ch].cook (drive[ch].getCurrentValue(), width[ch].getCurrentValue(), sat[ch].getCurrentValue(), wasV1);
    164         hProcs[ch].reset();
    165     }
    166 
    167     biasAngle.resize ((size_t) numChannels, 0.0);
    168     makeup.reset (numSteps);
    169 
    170     dcBlocker.resize ((size_t) numChannels);
    171     for (auto& filt : dcBlocker)
    172         filt.prepare (sampleRate, dcFreq);
    173 
    174     doubleBuffer.setSize (numChannels, samplesPerBlock);
    175     bypass.prepare (samplesPerBlock, numChannels, bypass.toBool (onOffParam));
    176 
    177 #if HYSTERESIS_USE_SIMD
    178     const auto maxOSBlockSize = (uint32) samplesPerBlock * 16;
    179     const auto numVecChannels = chowdsp::Math::ceiling_divide ((size_t) numChannels, Vec2::size);
    180     interleavedBlock = chowdsp::AudioBlock<Vec2> (interleavedBlockData, numVecChannels, maxOSBlockSize);
    181     zeroBlock = chowdsp::AudioBlock<double> (zeroData, Vec2::size, maxOSBlockSize);
    182     zeroBlock.clear();
    183 
    184     channelPointers.resize (numVecChannels * Vec2::size);
    185 #endif
    186 }
    187 
    188 void HysteresisProcessor::releaseResources()
    189 {
    190     osManager.reset();
    191 }
    192 
    193 float HysteresisProcessor::getLatencySamples() const noexcept
    194 {
    195     // latency of oversampling + fudge factor for hysteresis
    196     return onOffParam->load() == 1.0f ? osManager.getLatencySamples() + 1.4f // on
    197                                       : 0.0f; // off
    198 }
    199 
    200 void HysteresisProcessor::processBlock (AudioBuffer<float>& buffer)
    201 {
    202     const auto numChannels = buffer.getNumChannels();
    203 
    204     if (! bypass.processBlockIn (buffer, bypass.toBool (onOffParam)))
    205         return;
    206 
    207     setSolver ((int) *modeParam);
    208     setDrive (*driveParam);
    209     setSaturation (*satParam);
    210     setWidth (1.0f - *widthParam);
    211     makeup.setTargetValue (calcMakeup());
    212     setOversampling();
    213 
    214     bool needsSmoothing = drive[0].isSmoothing() || width[0].isSmoothing() || sat[0].isSmoothing() || wasV1 != useV1;
    215 
    216     if (useV1 != wasV1)
    217     {
    218         for (auto& hProc : hProcs)
    219             hProc.reset();
    220     }
    221 
    222     wasV1 = useV1;
    223 
    224     // clip input to avoid unstable hysteresis
    225     for (int ch = 0; ch < numChannels; ++ch)
    226     {
    227         auto* bufferPtr = buffer.getWritePointer (ch);
    228         FloatVectorOperations::clip (bufferPtr,
    229                                      bufferPtr,
    230                                      -clipLevel,
    231                                      clipLevel,
    232                                      buffer.getNumSamples());
    233     }
    234 
    235     doubleBuffer.makeCopyOf (buffer, true);
    236 
    237     dsp::AudioBlock<double> block (doubleBuffer);
    238     dsp::AudioBlock<double> osBlock = osManager.processSamplesUp (block);
    239 
    240 #if HYSTERESIS_USE_SIMD
    241     const auto n = osBlock.getNumSamples();
    242     auto* inout = channelPointers.data();
    243     const auto numChannelsPadded = channelPointers.size();
    244     for (size_t ch = 0; ch < numChannelsPadded; ++ch)
    245         inout[ch] = (ch < osBlock.getNumChannels() ? const_cast<double*> (osBlock.getChannelPointer (ch)) : zeroBlock.getChannelPointer (ch % Vec2::size));
    246 
    247     // interleave channels
    248     for (size_t ch = 0; ch < numChannelsPadded; ch += Vec2::size)
    249     {
    250         auto* simdBlockData = reinterpret_cast<double*> (interleavedBlock.getChannelPointer (ch / Vec2::size));
    251         interleaveSamples (&inout[ch], simdBlockData, static_cast<int> (n), static_cast<int> (Vec2::size));
    252     }
    253 
    254     auto&& processBlock = interleavedBlock.getSubBlock (0, n);
    255 
    256     using ProcessType = Vec2;
    257 #else
    258     auto&& processBlock = osBlock;
    259 
    260     using ProcessType = double;
    261 #endif
    262 
    263     if (useV1)
    264     {
    265         if (needsSmoothing)
    266             processSmoothV1<ProcessType> (processBlock);
    267         else
    268             processV1<ProcessType> (processBlock);
    269     }
    270     else
    271     {
    272         switch (solver)
    273         {
    274             case RK2:
    275                 if (needsSmoothing)
    276                     processSmooth<RK2, ProcessType> (processBlock);
    277                 else
    278                     process<RK2, ProcessType> (processBlock);
    279                 break;
    280             case RK4:
    281                 if (needsSmoothing)
    282                     processSmooth<RK4, ProcessType> (processBlock);
    283                 else
    284                     process<RK4, ProcessType> (processBlock);
    285                 break;
    286             case NR4:
    287                 if (needsSmoothing)
    288                     processSmooth<NR4, ProcessType> (processBlock);
    289                 else
    290                     process<NR4, ProcessType> (processBlock);
    291                 break;
    292             case NR8:
    293                 if (needsSmoothing)
    294                     processSmooth<NR8, ProcessType> (processBlock);
    295                 else
    296                     process<NR8, ProcessType> (processBlock);
    297                 break;
    298             case STN:
    299                 if (needsSmoothing)
    300                     processSmooth<STN, ProcessType> (processBlock);
    301                 else
    302                     process<STN, ProcessType> (processBlock);
    303                 break;
    304             default:
    305                 jassertfalse; // unknown solver!
    306         };
    307     }
    308 
    309 #if HYSTERESIS_USE_SIMD
    310     // de-interleave channels
    311     for (size_t ch = 0; ch < numChannelsPadded; ch += Vec2::size)
    312     {
    313         auto* simdBlockData = reinterpret_cast<double*> (interleavedBlock.getChannelPointer (ch / Vec2::size));
    314         deinterleaveSamples (simdBlockData,
    315                              const_cast<double**> (&inout[ch]),
    316                              static_cast<int> (n),
    317                              static_cast<int> (Vec2::size));
    318     }
    319 #endif
    320 
    321     osManager.processSamplesDown (block);
    322 
    323     buffer.makeCopyOf (doubleBuffer, true);
    324     applyDCBlockers (buffer);
    325 
    326     bypass.processBlockOut (buffer, bypass.toBool (onOffParam));
    327 }
    328 
    329 template <typename T, typename SmoothType>
    330 void applyMakeup (chowdsp::AudioBlock<T>& block, SmoothType& makeup)
    331 {
    332 #if HYSTERESIS_USE_SIMD
    333     const auto numSamples = block.getNumSamples();
    334     const auto numChannels = block.getNumChannels();
    335 
    336     if (makeup.isSmoothing())
    337     {
    338         for (size_t i = 0; i < numSamples; ++i)
    339         {
    340             const auto scaler = makeup.getNextValue();
    341 
    342             for (size_t ch = 0; ch < numChannels; ++ch)
    343                 block.getChannelPointer (ch)[i] *= scaler;
    344         }
    345     }
    346     else
    347     {
    348         const auto scaler = makeup.getTargetValue();
    349         for (size_t ch = 0; ch < numChannels; ++ch)
    350         {
    351             auto* x = block.getChannelPointer (ch);
    352             for (size_t i = 0; i < numSamples; ++i)
    353                 x[i] *= scaler;
    354         }
    355     }
    356 #else
    357     block *= makeup;
    358 #endif
    359 }
    360 
    361 template <SolverType solverType, typename T>
    362 void HysteresisProcessor::process (chowdsp::AudioBlock<T>& block)
    363 {
    364     const auto numChannels = block.getNumChannels();
    365     const auto numSamples = block.getNumSamples();
    366 
    367     for (size_t channel = 0; channel < numChannels; ++channel)
    368     {
    369         auto* x = block.getChannelPointer (channel);
    370         auto& hProc = hProcs[channel];
    371         for (size_t samp = 0; samp < numSamples; samp++)
    372             x[samp] = hProc.process<solverType> (x[samp]);
    373     }
    374 
    375     applyMakeup<T> (block, makeup);
    376 }
    377 
    378 template <SolverType solverType, typename T>
    379 void HysteresisProcessor::processSmooth (chowdsp::AudioBlock<T>& block)
    380 {
    381     const auto numChannels = block.getNumChannels();
    382     const auto numSamples = block.getNumSamples();
    383 
    384     for (size_t channel = 0; channel < numChannels; ++channel)
    385     {
    386         auto* x = block.getChannelPointer (channel);
    387         auto& hProc = hProcs[channel];
    388         for (size_t samp = 0; samp < numSamples; samp++)
    389         {
    390             hProc.cook (drive[channel].getNextValue(), width[channel].getNextValue(), sat[channel].getNextValue(), false);
    391             x[samp] = hProc.process<solverType> (x[samp]);
    392         }
    393     }
    394 
    395     applyMakeup<T> (block, makeup);
    396 }
    397 
    398 template <typename T>
    399 void HysteresisProcessor::processV1 (chowdsp::AudioBlock<T>& block)
    400 {
    401     const auto numChannels = block.getNumChannels();
    402     const auto numSamples = block.getNumSamples();
    403     const auto angleDelta = MathConstants<double>::twoPi * biasFreq / (fs * osManager.getOSFactor());
    404 
    405     for (size_t channel = 0; channel < numChannels; ++channel)
    406     {
    407         auto* x = block.getChannelPointer (channel);
    408         auto& hProc = hProcs[channel];
    409         auto& bAngle = biasAngle[channel];
    410         const auto bAngleMult = biasGain * (1.0 - width[channel].getCurrentValue());
    411         for (size_t samp = 0; samp < numSamples; samp++)
    412         {
    413             auto bias = bAngleMult * std::sin (bAngle);
    414             bAngle += angleDelta;
    415             bAngle -= MathConstants<double>::twoPi * (bAngle >= MathConstants<double>::twoPi);
    416 
    417             x[samp] = hProc.process<RK4> ((x[samp] + bias) * 10000.0) * v1Norm;
    418         }
    419     }
    420 }
    421 
    422 template <typename T>
    423 void HysteresisProcessor::processSmoothV1 (chowdsp::AudioBlock<T>& block)
    424 {
    425     const auto numChannels = block.getNumChannels();
    426     const auto numSamples = block.getNumSamples();
    427     const auto angleDelta = MathConstants<double>::twoPi * biasFreq / (fs * osManager.getOSFactor());
    428 
    429     for (size_t channel = 0; channel < numChannels; ++channel)
    430     {
    431         auto* x = block.getChannelPointer (channel);
    432         auto& hProc = hProcs[channel];
    433         auto& bAngle = biasAngle[channel];
    434         for (size_t samp = 0; samp < numSamples; samp++)
    435         {
    436             hProc.cook (drive[channel].getNextValue(), width[channel].getNextValue(), sat[channel].getNextValue(), true);
    437 
    438             auto bias = biasGain * (1.0 - width[channel].getCurrentValue()) * std::sin (bAngle);
    439             bAngle += angleDelta;
    440             bAngle -= MathConstants<double>::twoPi * (bAngle >= MathConstants<double>::twoPi);
    441 
    442             x[samp] = hProc.process<RK4> ((x[samp] + bias) * 10000.0) * v1Norm;
    443         }
    444     }
    445 }
    446 
    447 void HysteresisProcessor::applyDCBlockers (AudioBuffer<float>& buffer)
    448 {
    449     for (int ch = 0; ch < buffer.getNumChannels(); ++ch)
    450         dcBlocker[(size_t) ch].processBlock (buffer.getWritePointer (ch), buffer.getNumSamples());
    451 }