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