Alc/mastering.cpp

   1 #include "config.h"
   2
   3 #include <math.h>
   4
   5 #include "mastering.h"
   6 #include "alu.h"
   7 #include "almalloc.h"
   8 #include "math_defs.h"
   9
  10
  11 /* Early MSVC lacks round/roundf */
  12 #if defined(_MSC_VER) && _MSC_VER < 1800
  13 static double round(double val)
  14 {
  15     if(val < 0.0)
  16         return ceil(val-0.5);
  17     return floor(val+0.5);
  18 }
  19 #define roundf(f) ((float)round((float)(f)))
  20 #endif
  21
  22
  23 /* These structures assume BUFFERSIZE is a power of 2. */
  24 static_assert((BUFFERSIZE & (BUFFERSIZE-1)) == 0, "BUFFERSIZE is not a power of 2");
  25
  26 typedef struct SlidingHold {
  27     ALfloat Values[BUFFERSIZE];
  28     ALsizei Expiries[BUFFERSIZE];
  29     ALsizei LowerIndex;
  30     ALsizei UpperIndex;
  31     ALsizei Length;
  32 } SlidingHold;
  33
  34 /* General topology and basic automation was based on the following paper:
  35  *
  36  *   D. Giannoulis, M. Massberg and J. D. Reiss,
  37  *   "Parameter Automation in a Dynamic Range Compressor,"
  38  *   Journal of the Audio Engineering Society, v61 (10), Oct. 2013
  39  *
  40  * Available (along with supplemental reading) at:
  41  *
  42  *   http://c4dm.eecs.qmul.ac.uk/audioengineering/compressors/
  43  */
  44 typedef struct Compressor {
  45     ALsizei NumChans;
  46     ALuint SampleRate;
  47
  48     struct {
  49         ALuint Knee : 1;
  50         ALuint Attack : 1;
  51         ALuint Release : 1;
  52         ALuint PostGain : 1;
  53         ALuint Declip : 1;
  54     } Auto;
  55
  56     ALsizei LookAhead;
  57
  58     ALfloat PreGain;
  59     ALfloat PostGain;
  60
  61     ALfloat Threshold;
  62     ALfloat Slope;
  63     ALfloat Knee;
  64
  65     ALfloat Attack;
  66     ALfloat Release;
  67
  68     alignas(16) ALfloat SideChain[2*BUFFERSIZE];
  69     alignas(16) ALfloat CrestFactor[BUFFERSIZE];
  70
  71     SlidingHold *Hold;
  72     ALfloat (*Delay)[BUFFERSIZE];
  73     ALsizei DelayIndex;
  74
  75     ALfloat CrestCoeff;
  76     ALfloat GainEstimate;
  77     ALfloat AdaptCoeff;
  78
  79     ALfloat LastPeakSq;
  80     ALfloat LastRmsSq;
  81     ALfloat LastRelease;
  82     ALfloat LastAttack;
  83     ALfloat LastGainDev;
  84 } Compressor;
  85
  86
  87 /* This sliding hold follows the input level with an instant attack and a
  88  * fixed duration hold before an instant release to the next highest level.
  89  * It is a sliding window maximum (descending maxima) implementation based on
  90  * Richard Harter's ascending minima algorithm available at:
  91  *
  92  *   http://www.richardhartersworld.com/cri/2001/slidingmin.html
  93  */
  94 static ALfloat UpdateSlidingHold(SlidingHold *Hold, const ALsizei i, const ALfloat in)
  95 {
  96     const ALsizei mask = BUFFERSIZE - 1;
  97     const ALsizei length = Hold->Length;
  98     ALfloat *RESTRICT values = Hold->Values;
  99     ALsizei *RESTRICT expiries = Hold->Expiries;
 100     ALsizei lowerIndex = Hold->LowerIndex;
 101     ALsizei upperIndex = Hold->UpperIndex;
 102
 103     if(i >= expiries[upperIndex])
 104         upperIndex = (upperIndex + 1) & mask;
 105
 106     if(in >= values[upperIndex])
 107     {
 108         values[upperIndex] = in;
 109         expiries[upperIndex] = i + length;
 110         lowerIndex = upperIndex;
 111     }
 112     else
 113     {
 114         do {
 115             do {
 116                 if(!(in >= values[lowerIndex]))
 117                     goto found_place;
 118             } while(lowerIndex--);
 119             lowerIndex = mask;
 120         } while(1);
 121     found_place:
 122
 123         lowerIndex = (lowerIndex + 1) & mask;
 124         values[lowerIndex] = in;
 125         expiries[lowerIndex] = i + length;
 126     }
 127
 128     Hold->LowerIndex = lowerIndex;
 129     Hold->UpperIndex = upperIndex;
 130
 131     return values[upperIndex];
 132 }
 133
 134 static void ShiftSlidingHold(SlidingHold *Hold, const ALsizei n)
 135 {
 136     const ALsizei lowerIndex = Hold->LowerIndex;
 137     ALsizei *RESTRICT expiries = Hold->Expiries;
 138     ALsizei i = Hold->UpperIndex;
 139
 140     if(lowerIndex < i)
 141     {
 142         for(;i < BUFFERSIZE;i++)
 143             expiries[i] -= n;
 144         i = 0;
 145     }
 146     for(;i < lowerIndex;i++)
 147         expiries[i] -= n;
 148
 149     expiries[i] -= n;
 150 }
 151
 152 /* Multichannel compression is linked via the absolute maximum of all
 153  * channels.
 154  */
 155 static void LinkChannels(Compressor *Comp, const ALsizei SamplesToDo, ALfloat (*RESTRICT OutBuffer)[BUFFERSIZE])
 156 {
 157     const ALsizei index = Comp->LookAhead;
 158     const ALsizei numChans = Comp->NumChans;
 159     ALfloat *RESTRICT sideChain = Comp->SideChain;
 160     ALsizei c, i;
 161
 162     ASSUME(SamplesToDo > 0);
 163     ASSUME(numChans > 0);
 164
 165     for(i = 0;i < SamplesToDo;i++)
 166         sideChain[index + i] = 0.0f;
 167
 168     for(c = 0;c < numChans;c++)
 169     {
 170         ALsizei offset = index;
 171         for(i = 0;i < SamplesToDo;i++)
 172         {
 173             sideChain[offset] = maxf(sideChain[offset], fabsf(OutBuffer[c][i]));
 174             ++offset;
 175         }
 176     }
 177 }
 178
 179 /* This calculates the squared crest factor of the control signal for the
 180  * basic automation of the attack/release times.  As suggested by the paper,
 181  * it uses an instantaneous squared peak detector and a squared RMS detector
 182  * both with 200ms release times.
 183  */
 184 static void CrestDetector(Compressor *Comp, const ALsizei SamplesToDo)
 185 {
 186     const ALfloat a_crest = Comp->CrestCoeff;
 187     const ALsizei index = Comp->LookAhead;
 188     const ALfloat *RESTRICT sideChain = Comp->SideChain;
 189     ALfloat *RESTRICT crestFactor = Comp->CrestFactor;
 190     ALfloat y2_peak = Comp->LastPeakSq;
 191     ALfloat y2_rms = Comp->LastRmsSq;
 192     ALsizei i;
 193
 194     ASSUME(SamplesToDo > 0);
 195
 196     for(i = 0;i < SamplesToDo;i++)
 197     {
 198         ALfloat x_abs = sideChain[index + i];
 199         ALfloat x2 = maxf(0.000001f, x_abs * x_abs);
 200
 201         y2_peak = maxf(x2, lerp(x2, y2_peak, a_crest));
 202         y2_rms = lerp(x2, y2_rms, a_crest);
 203         crestFactor[i] = y2_peak / y2_rms;
 204     }
 205
 206     Comp->LastPeakSq = y2_peak;
 207     Comp->LastRmsSq = y2_rms;
 208 }
 209
 210 /* The side-chain starts with a simple peak detector (based on the absolute
 211  * value of the incoming signal) and performs most of its operations in the
 212  * log domain.
 213  */
 214 static void PeakDetector(Compressor *Comp, const ALsizei SamplesToDo)
 215 {
 216     const ALsizei index = Comp->LookAhead;
 217     ALfloat *RESTRICT sideChain = Comp->SideChain;
 218     ALsizei i;
 219
 220     ASSUME(SamplesToDo > 0);
 221
 222     for(i = 0;i < SamplesToDo;i++)
 223     {
 224         const ALuint offset = index + i;
 225         const ALfloat x_abs = sideChain[offset];
 226
 227         sideChain[offset] = logf(maxf(0.000001f, x_abs));
 228     }
 229 }
 230
 231 /* An optional hold can be used to extend the peak detector so it can more
 232  * solidly detect fast transients.  This is best used when operating as a
 233  * limiter.
 234  */
 235 static void PeakHoldDetector(Compressor *Comp, const ALsizei SamplesToDo)
 236 {
 237     const ALsizei index = Comp->LookAhead;
 238     ALfloat *RESTRICT sideChain = Comp->SideChain;
 239     SlidingHold *hold = Comp->Hold;
 240     ALsizei i;
 241
 242     ASSUME(SamplesToDo > 0);
 243
 244     for(i = 0;i < SamplesToDo;i++)
 245     {
 246         const ALsizei offset = index + i;
 247         const ALfloat x_abs = sideChain[offset];
 248         const ALfloat x_G = logf(maxf(0.000001f, x_abs));
 249
 250         sideChain[offset] = UpdateSlidingHold(hold, i, x_G);
 251     }
 252
 253     ShiftSlidingHold(hold, SamplesToDo);
 254 }
 255
 256 /* This is the heart of the feed-forward compressor.  It operates in the log
 257  * domain (to better match human hearing) and can apply some basic automation
 258  * to knee width, attack/release times, make-up/post gain, and clipping
 259  * reduction.
 260  */
 261 static void GainCompressor(Compressor *Comp, const ALsizei SamplesToDo)
 262 {
 263     const bool autoKnee = Comp->Auto.Knee;
 264     const bool autoAttack = Comp->Auto.Attack;
 265     const bool autoRelease = Comp->Auto.Release;
 266     const bool autoPostGain = Comp->Auto.PostGain;
 267     const bool autoDeclip = Comp->Auto.Declip;
 268     const ALsizei lookAhead = Comp->LookAhead;
 269     const ALfloat threshold = Comp->Threshold;
 270     const ALfloat slope = Comp->Slope;
 271     const ALfloat attack = Comp->Attack;
 272     const ALfloat release = Comp->Release;
 273     const ALfloat c_est = Comp->GainEstimate;
 274     const ALfloat a_adp = Comp->AdaptCoeff;
 275     const ALfloat *RESTRICT crestFactor = Comp->CrestFactor;
 276     ALfloat *RESTRICT sideChain = Comp->SideChain;
 277     ALfloat postGain = Comp->PostGain;
 278     ALfloat knee = Comp->Knee;
 279     ALfloat t_att = attack;
 280     ALfloat t_rel = release - attack;
 281     ALfloat a_att = expf(-1.0f / t_att);
 282     ALfloat a_rel = expf(-1.0f / t_rel);
 283     ALfloat y_1 = Comp->LastRelease;
 284     ALfloat y_L = Comp->LastAttack;
 285     ALfloat c_dev = Comp->LastGainDev;
 286     ALsizei i;
 287
 288     ASSUME(SamplesToDo > 0);
 289
 290     for(i = 0;i < SamplesToDo;i++)
 291     {
 292         const ALfloat y2_crest = crestFactor[i];
 293         const ALfloat x_G = sideChain[lookAhead + i];
 294         const ALfloat x_over = x_G - threshold;
 295         ALfloat knee_h;
 296         ALfloat y_G;
 297         ALfloat x_L;
 298
 299         if(autoKnee)
 300             knee = maxf(0.0f, 2.5f * (c_dev + c_est));
 301         knee_h = 0.5f * knee;
 302
 303         /* This is the gain computer.  It applies a static compression curve
 304          * to the control signal.
 305          */
 306         if(x_over <= -knee_h)
 307             y_G = 0.0f;
 308         else if(fabsf(x_over) < knee_h)
 309             y_G = (x_over + knee_h) * (x_over + knee_h) / (2.0f * knee);
 310         else
 311             y_G = x_over;
 312
 313         x_L = -slope * y_G;
 314
 315         if(autoAttack)
 316         {
 317             t_att = 2.0f * attack / y2_crest;
 318             a_att = expf(-1.0f / t_att);
 319         }
 320
 321         if(autoRelease)
 322         {
 323             t_rel = 2.0f * release / y2_crest - t_att;
 324             a_rel = expf(-1.0f / t_rel);
 325         }
 326
 327         /* Gain smoothing (ballistics) is done via a smooth decoupled peak
 328          * detector.  The attack time is subtracted from the release time
 329          * above to compensate for the chained operating mode.
 330          */
 331         y_1 = maxf(x_L, lerp(x_L, y_1, a_rel));
 332         y_L = lerp(y_1, y_L, a_att);
 333
 334         /* Knee width and make-up gain automation make use of a smoothed
 335          * measurement of deviation between the control signal and estimate.
 336          * The estimate is also used to bias the measurement to hot-start its
 337          * average.
 338          */
 339         c_dev = lerp(-y_L - c_est, c_dev, a_adp);
 340
 341         if(autoPostGain)
 342         {
 343             /* Clipping reduction is only viable when make-up gain is being
 344              * automated.  It modifies the deviation to further attenuate the
 345              * control signal when clipping is detected.  The adaptation
 346              * time is sufficiently long enough to suppress further clipping
 347              * at the same output level.
 348              */
 349             if(autoDeclip)
 350                 c_dev = maxf(c_dev, sideChain[i] - y_L - threshold - c_est);
 351
 352             postGain = -(c_dev + c_est);
 353         }
 354
 355         sideChain[i] = expf(postGain - y_L);
 356     }
 357
 358     Comp->LastRelease = y_1;
 359     Comp->LastAttack = y_L;
 360     Comp->LastGainDev = c_dev;
 361 }
 362
 363 /* Combined with the hold time, a look-ahead delay can improve handling of
 364  * fast transients by allowing the envelope time to converge prior to
 365  * reaching the offending impulse.  This is best used when operating as a
 366  * limiter.
 367  */
 368 static void SignalDelay(Compressor *Comp, const ALsizei SamplesToDo, ALfloat (*RESTRICT OutBuffer)[BUFFERSIZE])
 369 {
 370     const ALsizei mask = BUFFERSIZE - 1;
 371     const ALsizei numChans = Comp->NumChans;
 372     const ALsizei indexIn = Comp->DelayIndex;
 373     const ALsizei indexOut = Comp->DelayIndex - Comp->LookAhead;
 374     ALfloat (*RESTRICT delay)[BUFFERSIZE] = Comp->Delay;
 375     ALsizei c, i;
 376
 377     ASSUME(SamplesToDo > 0);
 378     ASSUME(numChans > 0);
 379
 380     for(c = 0;c < numChans;c++)
 381     {
 382         for(i = 0;i < SamplesToDo;i++)
 383         {
 384             ALfloat sig = OutBuffer[c][i];
 385
 386             OutBuffer[c][i] = delay[c][(indexOut + i) & mask];
 387             delay[c][(indexIn + i) & mask] = sig;
 388         }
 389     }
 390
 391     Comp->DelayIndex = (indexIn + SamplesToDo) & mask;
 392 }
 393
 394 /* The compressor is initialized with the following settings:
 395  *
 396  *   NumChans       - Number of channels to process.
 397  *   SampleRate     - Sample rate to process.
 398  *   AutoKnee       - Whether to automate the knee width parameter.
 399  *   AutoAttack     - Whether to automate the attack time parameter.
 400  *   AutoRelease    - Whether to automate the release time parameter.
 401  *   AutoPostGain   - Whether to automate the make-up (post) gain parameter.
 402  *   AutoDeclip     - Whether to automate clipping reduction.  Ignored when
 403  *                    not automating make-up gain.
 404  *   LookAheadTime  - Look-ahead time (in seconds).
 405  *   HoldTime       - Peak hold-time (in seconds).
 406  *   PreGainDb      - Gain applied before detection (in dB).
 407  *   PostGainDb     - Make-up gain applied after compression (in dB).
 408  *   ThresholdDb    - Triggering threshold (in dB).
 409  *   Ratio          - Compression ratio (x:1).  Set to INFINITY for true
 410  *                    limiting.  Ignored when automating knee width.
 411  *   KneeDb         - Knee width (in dB).  Ignored when automating knee
 412  *                    width.
 413  *   AttackTimeMin  - Attack time (in seconds).  Acts as a maximum when
 414  *                    automating attack time.
 415  *   ReleaseTimeMin - Release time (in seconds).  Acts as a maximum when
 416  *                    automating release time.
 417  */
 418 Compressor* CompressorInit(const ALsizei NumChans, const ALuint SampleRate,
 419                            const ALboolean AutoKnee, const ALboolean AutoAttack,
 420                            const ALboolean AutoRelease, const ALboolean AutoPostGain,
 421                            const ALboolean AutoDeclip, const ALfloat LookAheadTime,
 422                            const ALfloat HoldTime, const ALfloat PreGainDb,
 423                            const ALfloat PostGainDb, const ALfloat ThresholdDb,
 424                            const ALfloat Ratio, const ALfloat KneeDb,
 425                            const ALfloat AttackTime, const ALfloat ReleaseTime)
 426 {
 427     Compressor *Comp;
 428     ALsizei lookAhead;
 429     ALsizei hold;
 430     size_t size;
 431
 432     lookAhead = (ALsizei)clampf(roundf(LookAheadTime*SampleRate), 0.0f, BUFFERSIZE-1);
 433     hold = (ALsizei)clampf(roundf(HoldTime*SampleRate), 0.0f, BUFFERSIZE-1);
 434     /* The sliding hold implementation doesn't handle a length of 1. A 1-sample
 435      * hold is useless anyway, it would only ever give back what was just given
 436      * to it.
 437      */
 438     if(hold == 1)
 439         hold = 0;
 440
 441     size = sizeof(*Comp);
 442     if(lookAhead > 0)
 443     {
 444         size += sizeof(*Comp->Delay) * NumChans;
 445         if(hold > 0)
 446             size += sizeof(*Comp->Hold);
 447     }
 448
 449     Comp = static_cast<Compressor*>(al_calloc(16, size));
 450     Comp->NumChans = NumChans;
 451     Comp->SampleRate = SampleRate;
 452     Comp->Auto.Knee = AutoKnee;
 453     Comp->Auto.Attack = AutoAttack;
 454     Comp->Auto.Release = AutoRelease;
 455     Comp->Auto.PostGain = AutoPostGain;
 456     Comp->Auto.Declip = AutoPostGain && AutoDeclip;
 457     Comp->LookAhead = lookAhead;
 458     Comp->PreGain = powf(10.0f, PreGainDb / 20.0f);
 459     Comp->PostGain = PostGainDb * logf(10.0f) / 20.0f;
 460     Comp->Threshold = ThresholdDb * logf(10.0f) / 20.0f;
 461     Comp->Slope = 1.0f / maxf(1.0f, Ratio) - 1.0f;
 462     Comp->Knee = maxf(0.0f, KneeDb * logf(10.0f) / 20.0f);
 463     Comp->Attack = maxf(1.0f, AttackTime * SampleRate);
 464     Comp->Release = maxf(1.0f, ReleaseTime * SampleRate);
 465
 466     /* Knee width automation actually treats the compressor as a limiter.  By
 467      * varying the knee width, it can effectively be seen as applying
 468      * compression over a wide range of ratios.
 469      */
 470     if(AutoKnee)
 471         Comp->Slope = -1.0f;
 472
 473     if(lookAhead > 0)
 474     {
 475         if(hold > 0)
 476         {
 477             Comp->Hold = (SlidingHold*)(Comp + 1);
 478             Comp->Hold->Values[0] = -HUGE_VALF;
 479             Comp->Hold->Expiries[0] = hold;
 480             Comp->Hold->Length = hold;
 481             Comp->Delay = (ALfloat(*)[BUFFERSIZE])(Comp->Hold + 1);
 482         }
 483         else
 484         {
 485             Comp->Delay = (ALfloat(*)[BUFFERSIZE])(Comp + 1);
 486         }
 487     }
 488
 489     Comp->CrestCoeff = expf(-1.0f / (0.200f * SampleRate)); // 200ms
 490     Comp->GainEstimate = Comp->Threshold * -0.5f * Comp->Slope;
 491     Comp->AdaptCoeff = expf(-1.0f / (2.0f * SampleRate)); // 2s
 492
 493     return Comp;
 494 }
 495
 496 void ApplyCompression(Compressor *Comp, const ALsizei SamplesToDo, ALfloat (*RESTRICT OutBuffer)[BUFFERSIZE])
 497 {
 498     const ALsizei numChans = Comp->NumChans;
 499     const ALfloat preGain = Comp->PreGain;
 500     ALfloat *RESTRICT sideChain;
 501     ALsizei c, i;
 502
 503     ASSUME(SamplesToDo > 0);
 504     ASSUME(numChans > 0);
 505
 506     if(preGain != 1.0f)
 507     {
 508         for(c = 0;c < numChans;c++)
 509         {
 510             for(i = 0;i < SamplesToDo;i++)
 511                 OutBuffer[c][i] *= preGain;
 512         }
 513     }
 514
 515     LinkChannels(Comp, SamplesToDo, OutBuffer);
 516
 517     if(Comp->Auto.Attack || Comp->Auto.Release)
 518         CrestDetector(Comp, SamplesToDo);
 519
 520     if(Comp->Hold)
 521         PeakHoldDetector(Comp, SamplesToDo);
 522     else
 523         PeakDetector(Comp, SamplesToDo);
 524
 525     GainCompressor(Comp, SamplesToDo);
 526
 527     if(Comp->Delay)
 528         SignalDelay(Comp, SamplesToDo, OutBuffer);
 529
 530     sideChain = Comp->SideChain;
 531     for(c = 0;c < numChans;c++)
 532     {
 533         for(i = 0;i < SamplesToDo;i++)
 534             OutBuffer[c][i] *= sideChain[i];
 535     }
 536
 537     memmove(sideChain, sideChain+SamplesToDo, Comp->LookAhead*sizeof(ALfloat));
 538 }
 539
 540
 541 ALsizei GetCompressorLookAhead(const Compressor *Comp)
 542 { return Comp->LookAhead; }