core/mastering.cpp

   1
   2 #include "config.h"
   3
   4 #include "mastering.h"
   5
   6 #include <algorithm>
   7 #include <cmath>
   8 #include <cstddef>
   9 #include <functional>
  10 #include <iterator>
  11 #include <limits>
  12 #include <new>
  13
  14 #include "almalloc.h"
  15 #include "alnumeric.h"
  16 #include "alspan.h"
  17 #include "opthelpers.h"
  18
  19
  20 /* These structures assume BufferLineSize is a power of 2. */
  21 static_assert((BufferLineSize & (BufferLineSize-1)) == 0, "BufferLineSize is not a power of 2");
  22
  23 struct SlidingHold {
  24     alignas(16) float mValues[BufferLineSize];
  25     uint mExpiries[BufferLineSize];
  26     uint mLowerIndex;
  27     uint mUpperIndex;
  28     uint mLength;
  29 };
  30
  31
  32 namespace {
  33
  34 using namespace std::placeholders;
  35
  36 /* This sliding hold follows the input level with an instant attack and a
  37  * fixed duration hold before an instant release to the next highest level.
  38  * It is a sliding window maximum (descending maxima) implementation based on
  39  * Richard Harter's ascending minima algorithm available at:
  40  *
  41  *   http://www.richardhartersworld.com/cri/2001/slidingmin.html
  42  */
  43 float UpdateSlidingHold(SlidingHold *Hold, const uint i, const float in)
  44 {
  45     static constexpr uint mask{BufferLineSize - 1};
  46     const uint length{Hold->mLength};
  47     float (&values)[BufferLineSize] = Hold->mValues;
  48     uint (&expiries)[BufferLineSize] = Hold->mExpiries;
  49     uint lowerIndex{Hold->mLowerIndex};
  50     uint upperIndex{Hold->mUpperIndex};
  51
  52     if(i >= expiries[upperIndex])
  53         upperIndex = (upperIndex + 1) & mask;
  54
  55     if(in >= values[upperIndex])
  56     {
  57         values[upperIndex] = in;
  58         expiries[upperIndex] = i + length;
  59         lowerIndex = upperIndex;
  60     }
  61     else
  62     {
  63         do {
  64             do {
  65                 if(!(in >= values[lowerIndex]))
  66                     goto found_place;
  67             } while(lowerIndex--);
  68             lowerIndex = mask;
  69         } while(1);
  70     found_place:
  71
  72         lowerIndex = (lowerIndex + 1) & mask;
  73         values[lowerIndex] = in;
  74         expiries[lowerIndex] = i + length;
  75     }
  76
  77     Hold->mLowerIndex = lowerIndex;
  78     Hold->mUpperIndex = upperIndex;
  79
  80     return values[upperIndex];
  81 }
  82
  83 void ShiftSlidingHold(SlidingHold *Hold, const uint n)
  84 {
  85     auto exp_begin = std::begin(Hold->mExpiries) + Hold->mUpperIndex;
  86     auto exp_last = std::begin(Hold->mExpiries) + Hold->mLowerIndex;
  87     if(exp_last-exp_begin < 0)
  88     {
  89         std::transform(exp_begin, std::end(Hold->mExpiries), exp_begin,
  90             std::bind(std::minus<>{}, _1, n));
  91         exp_begin = std::begin(Hold->mExpiries);
  92     }
  93     std::transform(exp_begin, exp_last+1, exp_begin, std::bind(std::minus<>{}, _1, n));
  94 }
  95
  96
  97 /* Multichannel compression is linked via the absolute maximum of all
  98  * channels.
  99  */
 100 void LinkChannels(Compressor *Comp, const uint SamplesToDo, const FloatBufferLine *OutBuffer)
 101 {
 102     const size_t numChans{Comp->mNumChans};
 103
 104     ASSUME(SamplesToDo > 0);
 105     ASSUME(numChans > 0);
 106
 107     auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
 108     std::fill(side_begin, side_begin+SamplesToDo, 0.0f);
 109
 110     auto fill_max = [SamplesToDo,side_begin](const FloatBufferLine &input) -> void
 111     {
 112         const float *RESTRICT buffer{al::assume_aligned<16>(input.data())};
 113         auto max_abs = std::bind(maxf, _1, std::bind(static_cast<float(&)(float)>(std::fabs), _2));
 114         std::transform(side_begin, side_begin+SamplesToDo, buffer, side_begin, max_abs);
 115     };
 116     std::for_each(OutBuffer, OutBuffer+numChans, fill_max);
 117 }
 118
 119 /* This calculates the squared crest factor of the control signal for the
 120  * basic automation of the attack/release times.  As suggested by the paper,
 121  * it uses an instantaneous squared peak detector and a squared RMS detector
 122  * both with 200ms release times.
 123  */
 124 static void CrestDetector(Compressor *Comp, const uint SamplesToDo)
 125 {
 126     const float a_crest{Comp->mCrestCoeff};
 127     float y2_peak{Comp->mLastPeakSq};
 128     float y2_rms{Comp->mLastRmsSq};
 129
 130     ASSUME(SamplesToDo > 0);
 131
 132     auto calc_crest = [&y2_rms,&y2_peak,a_crest](const float x_abs) noexcept -> float
 133     {
 134         const float x2{clampf(x_abs * x_abs, 0.000001f, 1000000.0f)};
 135
 136         y2_peak = maxf(x2, lerpf(x2, y2_peak, a_crest));
 137         y2_rms = lerpf(x2, y2_rms, a_crest);
 138         return y2_peak / y2_rms;
 139     };
 140     auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
 141     std::transform(side_begin, side_begin+SamplesToDo, std::begin(Comp->mCrestFactor), calc_crest);
 142
 143     Comp->mLastPeakSq = y2_peak;
 144     Comp->mLastRmsSq = y2_rms;
 145 }
 146
 147 /* The side-chain starts with a simple peak detector (based on the absolute
 148  * value of the incoming signal) and performs most of its operations in the
 149  * log domain.
 150  */
 151 void PeakDetector(Compressor *Comp, const uint SamplesToDo)
 152 {
 153     ASSUME(SamplesToDo > 0);
 154
 155     /* Clamp the minimum amplitude to near-zero and convert to logarithm. */
 156     auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
 157     std::transform(side_begin, side_begin+SamplesToDo, side_begin,
 158         [](const float s) -> float { return std::log(maxf(0.000001f, s)); });
 159 }
 160
 161 /* An optional hold can be used to extend the peak detector so it can more
 162  * solidly detect fast transients.  This is best used when operating as a
 163  * limiter.
 164  */
 165 void PeakHoldDetector(Compressor *Comp, const uint SamplesToDo)
 166 {
 167     ASSUME(SamplesToDo > 0);
 168
 169     SlidingHold *hold{Comp->mHold};
 170     uint i{0};
 171     auto detect_peak = [&i,hold](const float x_abs) -> float
 172     {
 173         const float x_G{std::log(maxf(0.000001f, x_abs))};
 174         return UpdateSlidingHold(hold, i++, x_G);
 175     };
 176     auto side_begin = std::begin(Comp->mSideChain) + Comp->mLookAhead;
 177     std::transform(side_begin, side_begin+SamplesToDo, side_begin, detect_peak);
 178
 179     ShiftSlidingHold(hold, SamplesToDo);
 180 }
 181
 182 /* This is the heart of the feed-forward compressor.  It operates in the log
 183  * domain (to better match human hearing) and can apply some basic automation
 184  * to knee width, attack/release times, make-up/post gain, and clipping
 185  * reduction.
 186  */
 187 void GainCompressor(Compressor *Comp, const uint SamplesToDo)
 188 {
 189     const bool autoKnee{Comp->mAuto.Knee};
 190     const bool autoAttack{Comp->mAuto.Attack};
 191     const bool autoRelease{Comp->mAuto.Release};
 192     const bool autoPostGain{Comp->mAuto.PostGain};
 193     const bool autoDeclip{Comp->mAuto.Declip};
 194     const uint lookAhead{Comp->mLookAhead};
 195     const float threshold{Comp->mThreshold};
 196     const float slope{Comp->mSlope};
 197     const float attack{Comp->mAttack};
 198     const float release{Comp->mRelease};
 199     const float c_est{Comp->mGainEstimate};
 200     const float a_adp{Comp->mAdaptCoeff};
 201     const float *crestFactor{Comp->mCrestFactor};
 202     float postGain{Comp->mPostGain};
 203     float knee{Comp->mKnee};
 204     float t_att{attack};
 205     float t_rel{release - attack};
 206     float a_att{std::exp(-1.0f / t_att)};
 207     float a_rel{std::exp(-1.0f / t_rel)};
 208     float y_1{Comp->mLastRelease};
 209     float y_L{Comp->mLastAttack};
 210     float c_dev{Comp->mLastGainDev};
 211
 212     ASSUME(SamplesToDo > 0);
 213
 214     for(float &sideChain : al::span<float>{Comp->mSideChain, SamplesToDo})
 215     {
 216         if(autoKnee)
 217             knee = maxf(0.0f, 2.5f * (c_dev + c_est));
 218         const float knee_h{0.5f * knee};
 219
 220         /* This is the gain computer.  It applies a static compression curve
 221          * to the control signal.
 222          */
 223         const float x_over{std::addressof(sideChain)[lookAhead] - threshold};
 224         const float y_G{
 225             (x_over <= -knee_h) ? 0.0f :
 226             (std::fabs(x_over) < knee_h) ? (x_over + knee_h) * (x_over + knee_h) / (2.0f * knee) :
 227             x_over};
 228
 229         const float y2_crest{*(crestFactor++)};
 230         if(autoAttack)
 231         {
 232             t_att = 2.0f*attack/y2_crest;
 233             a_att = std::exp(-1.0f / t_att);
 234         }
 235         if(autoRelease)
 236         {
 237             t_rel = 2.0f*release/y2_crest - t_att;
 238             a_rel = std::exp(-1.0f / t_rel);
 239         }
 240
 241         /* Gain smoothing (ballistics) is done via a smooth decoupled peak
 242          * detector.  The attack time is subtracted from the release time
 243          * above to compensate for the chained operating mode.
 244          */
 245         const float x_L{-slope * y_G};
 246         y_1 = maxf(x_L, lerpf(x_L, y_1, a_rel));
 247         y_L = lerpf(y_1, y_L, a_att);
 248
 249         /* Knee width and make-up gain automation make use of a smoothed
 250          * measurement of deviation between the control signal and estimate.
 251          * The estimate is also used to bias the measurement to hot-start its
 252          * average.
 253          */
 254         c_dev = lerpf(-(y_L+c_est), c_dev, a_adp);
 255
 256         if(autoPostGain)
 257         {
 258             /* Clipping reduction is only viable when make-up gain is being
 259              * automated. It modifies the deviation to further attenuate the
 260              * control signal when clipping is detected. The adaptation time
 261              * is sufficiently long enough to suppress further clipping at the
 262              * same output level.
 263              */
 264             if(autoDeclip)
 265                 c_dev = maxf(c_dev, sideChain - y_L - threshold - c_est);
 266
 267             postGain = -(c_dev + c_est);
 268         }
 269
 270         sideChain = std::exp(postGain - y_L);
 271     }
 272
 273     Comp->mLastRelease = y_1;
 274     Comp->mLastAttack = y_L;
 275     Comp->mLastGainDev = c_dev;
 276 }
 277
 278 /* Combined with the hold time, a look-ahead delay can improve handling of
 279  * fast transients by allowing the envelope time to converge prior to
 280  * reaching the offending impulse.  This is best used when operating as a
 281  * limiter.
 282  */
 283 void SignalDelay(Compressor *Comp, const uint SamplesToDo, FloatBufferLine *OutBuffer)
 284 {
 285     const size_t numChans{Comp->mNumChans};
 286     const uint lookAhead{Comp->mLookAhead};
 287
 288     ASSUME(SamplesToDo > 0);
 289     ASSUME(numChans > 0);
 290     ASSUME(lookAhead > 0);
 291
 292     for(size_t c{0};c < numChans;c++)
 293     {
 294         float *inout{al::assume_aligned<16>(OutBuffer[c].data())};
 295         float *delaybuf{al::assume_aligned<16>(Comp->mDelay[c].data())};
 296
 297         auto inout_end = inout + SamplesToDo;
 298         if(SamplesToDo >= lookAhead) [[likely]]
 299         {
 300             auto delay_end = std::rotate(inout, inout_end - lookAhead, inout_end);
 301             std::swap_ranges(inout, delay_end, delaybuf);
 302         }
 303         else
 304         {
 305             auto delay_start = std::swap_ranges(inout, inout_end, delaybuf);
 306             std::rotate(delaybuf, delay_start, delaybuf + lookAhead);
 307         }
 308     }
 309 }
 310
 311 } // namespace
 312
 313
 314 std::unique_ptr<Compressor> Compressor::Create(const size_t NumChans, const float SampleRate,
 315     const bool AutoKnee, const bool AutoAttack, const bool AutoRelease, const bool AutoPostGain,
 316     const bool AutoDeclip, const float LookAheadTime, const float HoldTime, const float PreGainDb,
 317     const float PostGainDb, const float ThresholdDb, const float Ratio, const float KneeDb,
 318     const float AttackTime, const float ReleaseTime)
 319 {
 320     const auto lookAhead = static_cast<uint>(
 321         clampf(std::round(LookAheadTime*SampleRate), 0.0f, BufferLineSize-1));
 322     const auto hold = static_cast<uint>(
 323         clampf(std::round(HoldTime*SampleRate), 0.0f, BufferLineSize-1));
 324
 325     size_t size{sizeof(Compressor)};
 326     if(lookAhead > 0)
 327     {
 328         size += sizeof(*Compressor::mDelay) * NumChans;
 329         /* The sliding hold implementation doesn't handle a length of 1. A 1-
 330          * sample hold is useless anyway, it would only ever give back what was
 331          * just given to it.
 332          */
 333         if(hold > 1)
 334             size += sizeof(*Compressor::mHold);
 335     }
 336
 337     auto Comp = CompressorPtr{al::construct_at(static_cast<Compressor*>(al_calloc(16, size)))};
 338     Comp->mNumChans = NumChans;
 339     Comp->mAuto.Knee = AutoKnee;
 340     Comp->mAuto.Attack = AutoAttack;
 341     Comp->mAuto.Release = AutoRelease;
 342     Comp->mAuto.PostGain = AutoPostGain;
 343     Comp->mAuto.Declip = AutoPostGain && AutoDeclip;
 344     Comp->mLookAhead = lookAhead;
 345     Comp->mPreGain = std::pow(10.0f, PreGainDb / 20.0f);
 346     Comp->mPostGain = PostGainDb * std::log(10.0f) / 20.0f;
 347     Comp->mThreshold = ThresholdDb * std::log(10.0f) / 20.0f;
 348     Comp->mSlope = 1.0f / maxf(1.0f, Ratio) - 1.0f;
 349     Comp->mKnee = maxf(0.0f, KneeDb * std::log(10.0f) / 20.0f);
 350     Comp->mAttack = maxf(1.0f, AttackTime * SampleRate);
 351     Comp->mRelease = maxf(1.0f, ReleaseTime * SampleRate);
 352
 353     /* Knee width automation actually treats the compressor as a limiter. By
 354      * varying the knee width, it can effectively be seen as applying
 355      * compression over a wide range of ratios.
 356      */
 357     if(AutoKnee)
 358         Comp->mSlope = -1.0f;
 359
 360     if(lookAhead > 0)
 361     {
 362         if(hold > 1)
 363         {
 364             Comp->mHold = al::construct_at(reinterpret_cast<SlidingHold*>(Comp.get() + 1));
 365             Comp->mHold->mValues[0] = -std::numeric_limits<float>::infinity();
 366             Comp->mHold->mExpiries[0] = hold;
 367             Comp->mHold->mLength = hold;
 368             Comp->mDelay = reinterpret_cast<FloatBufferLine*>(Comp->mHold + 1);
 369         }
 370         else
 371             Comp->mDelay = reinterpret_cast<FloatBufferLine*>(Comp.get() + 1);
 372         std::uninitialized_fill_n(Comp->mDelay, NumChans, FloatBufferLine{});
 373     }
 374
 375     Comp->mCrestCoeff = std::exp(-1.0f / (0.200f * SampleRate)); // 200ms
 376     Comp->mGainEstimate = Comp->mThreshold * -0.5f * Comp->mSlope;
 377     Comp->mAdaptCoeff = std::exp(-1.0f / (2.0f * SampleRate)); // 2s
 378
 379     return Comp;
 380 }
 381
 382 Compressor::~Compressor()
 383 {
 384     if(mHold)
 385         al::destroy_at(mHold);
 386     mHold = nullptr;
 387     if(mDelay)
 388         al::destroy_n(mDelay, mNumChans);
 389     mDelay = nullptr;
 390 }
 391
 392
 393 void Compressor::process(const uint SamplesToDo, FloatBufferLine *OutBuffer)
 394 {
 395     const size_t numChans{mNumChans};
 396
 397     ASSUME(SamplesToDo > 0);
 398     ASSUME(numChans > 0);
 399
 400     const float preGain{mPreGain};
 401     if(preGain != 1.0f)
 402     {
 403         auto apply_gain = [SamplesToDo,preGain](FloatBufferLine &input) noexcept -> void
 404         {
 405             float *buffer{al::assume_aligned<16>(input.data())};
 406             std::transform(buffer, buffer+SamplesToDo, buffer,
 407                 std::bind(std::multiplies<float>{}, _1, preGain));
 408         };
 409         std::for_each(OutBuffer, OutBuffer+numChans, apply_gain);
 410     }
 411
 412     LinkChannels(this, SamplesToDo, OutBuffer);
 413
 414     if(mAuto.Attack || mAuto.Release)
 415         CrestDetector(this, SamplesToDo);
 416
 417     if(mHold)
 418         PeakHoldDetector(this, SamplesToDo);
 419     else
 420         PeakDetector(this, SamplesToDo);
 421
 422     GainCompressor(this, SamplesToDo);
 423
 424     if(mDelay)
 425         SignalDelay(this, SamplesToDo, OutBuffer);
 426
 427     const float (&sideChain)[BufferLineSize*2] = mSideChain;
 428     auto apply_comp = [SamplesToDo,&sideChain](FloatBufferLine &input) noexcept -> void
 429     {
 430         float *buffer{al::assume_aligned<16>(input.data())};
 431         const float *gains{al::assume_aligned<16>(&sideChain[0])};
 432         std::transform(gains, gains+SamplesToDo, buffer, buffer,
 433             std::bind(std::multiplies<float>{}, _1, _2));
 434     };
 435     std::for_each(OutBuffer, OutBuffer+numChans, apply_comp);
 436
 437     auto side_begin = std::begin(mSideChain) + SamplesToDo;
 438     std::copy(side_begin, side_begin+mLookAhead, std::begin(mSideChain));
 439 }