src/third_party/WebKit/Source/platform/audio/DynamicsCompressorKernel.cpp

   1 /*
   2  * Copyright (C) 2011 Google Inc. All rights reserved.
   3  *
   4  * Redistribution and use in source and binary forms, with or without
   5  * modification, are permitted provided that the following conditions
   6  * are met:
   7  *
   8  * 1.  Redistributions of source code must retain the above copyright
   9  *     notice, this list of conditions and the following disclaimer.
  10  * 2.  Redistributions in binary form must reproduce the above copyright
  11  *     notice, this list of conditions and the following disclaimer in the
  12  *     documentation and/or other materials provided with the distribution.
  13  * 3.  Neither the name of Apple Computer, Inc. ("Apple") nor the names of
  14  *     its contributors may be used to endorse or promote products derived
  15  *     from this software without specific prior written permission.
  16  *
  17  * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
  18  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  19  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  20  * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
  21  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  22  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  23  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  24  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  26  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  27  */
  28
  29 #include "config.h"
  30
  31 #if ENABLE(WEB_AUDIO)
  32
  33 #include "platform/audio/DynamicsCompressorKernel.h"
  34
  35 #include <algorithm>
  36 #include "platform/audio/AudioUtilities.h"
  37 #include "platform/audio/DenormalDisabler.h"
  38 #include "wtf/MathExtras.h"
  39
  40 using namespace std;
  41
  42 namespace WebCore {
  43
  44 using namespace AudioUtilities;
  45
  46 // Metering hits peaks instantly, but releases this fast (in seconds).
  47 const float meteringReleaseTimeConstant = 0.325f;
  48
  49 const float uninitializedValue = -1;
  50
  51 DynamicsCompressorKernel::DynamicsCompressorKernel(float sampleRate, unsigned numberOfChannels)
  52     : m_sampleRate(sampleRate)
  53     , m_lastPreDelayFrames(DefaultPreDelayFrames)
  54     , m_preDelayReadIndex(0)
  55     , m_preDelayWriteIndex(DefaultPreDelayFrames)
  56     , m_ratio(uninitializedValue)
  57     , m_slope(uninitializedValue)
  58     , m_linearThreshold(uninitializedValue)
  59     , m_dbThreshold(uninitializedValue)
  60     , m_dbKnee(uninitializedValue)
  61     , m_kneeThreshold(uninitializedValue)
  62     , m_kneeThresholdDb(uninitializedValue)
  63     , m_ykneeThresholdDb(uninitializedValue)
  64     , m_K(uninitializedValue)
  65 {
  66     setNumberOfChannels(numberOfChannels);
  67
  68     // Initializes most member variables
  69     reset();
  70
  71     m_meteringReleaseK = static_cast<float>(discreteTimeConstantForSampleRate(meteringReleaseTimeConstant, sampleRate));
  72 }
  73
  74 void DynamicsCompressorKernel::setNumberOfChannels(unsigned numberOfChannels)
  75 {
  76     if (m_preDelayBuffers.size() == numberOfChannels)
  77         return;
  78
  79     m_preDelayBuffers.clear();
  80     for (unsigned i = 0; i < numberOfChannels; ++i)
  81         m_preDelayBuffers.append(adoptPtr(new AudioFloatArray(MaxPreDelayFrames)));
  82 }
  83
  84 void DynamicsCompressorKernel::setPreDelayTime(float preDelayTime)
  85 {
  86     // Re-configure look-ahead section pre-delay if delay time has changed.
  87     unsigned preDelayFrames = preDelayTime * sampleRate();
  88     if (preDelayFrames > MaxPreDelayFrames - 1)
  89         preDelayFrames = MaxPreDelayFrames - 1;
  90
  91     if (m_lastPreDelayFrames != preDelayFrames) {
  92         m_lastPreDelayFrames = preDelayFrames;
  93         for (unsigned i = 0; i < m_preDelayBuffers.size(); ++i)
  94             m_preDelayBuffers[i]->zero();
  95
  96         m_preDelayReadIndex = 0;
  97         m_preDelayWriteIndex = preDelayFrames;
  98     }
  99 }
 100
 101 // Exponential curve for the knee.
 102 // It is 1st derivative matched at m_linearThreshold and asymptotically approaches the value m_linearThreshold + 1 / k.
 103 float DynamicsCompressorKernel::kneeCurve(float x, float k)
 104 {
 105     // Linear up to threshold.
 106     if (x < m_linearThreshold)
 107         return x;
 108
 109     return m_linearThreshold + (1 - expf(-k * (x - m_linearThreshold))) / k;
 110 }
 111
 112 // Full compression curve with constant ratio after knee.
 113 float DynamicsCompressorKernel::saturate(float x, float k)
 114 {
 115     float y;
 116
 117     if (x < m_kneeThreshold)
 118         y = kneeCurve(x, k);
 119     else {
 120         // Constant ratio after knee.
 121         float xDb = linearToDecibels(x);
 122         float yDb = m_ykneeThresholdDb + m_slope * (xDb - m_kneeThresholdDb);
 123
 124         y = decibelsToLinear(yDb);
 125     }
 126
 127     return y;
 128 }
 129
 130 // Approximate 1st derivative with input and output expressed in dB.
 131 // This slope is equal to the inverse of the compression "ratio".
 132 // In other words, a compression ratio of 20 would be a slope of 1/20.
 133 float DynamicsCompressorKernel::slopeAt(float x, float k)
 134 {
 135     if (x < m_linearThreshold)
 136         return 1;
 137
 138     float x2 = x * 1.001;
 139
 140     float xDb = linearToDecibels(x);
 141     float x2Db = linearToDecibels(x2);
 142
 143     float yDb = linearToDecibels(kneeCurve(x, k));
 144     float y2Db = linearToDecibels(kneeCurve(x2, k));
 145
 146     float m = (y2Db - yDb) / (x2Db - xDb);
 147
 148     return m;
 149 }
 150
 151 float DynamicsCompressorKernel::kAtSlope(float desiredSlope)
 152 {
 153     float xDb = m_dbThreshold + m_dbKnee;
 154     float x = decibelsToLinear(xDb);
 155
 156     // Approximate k given initial values.
 157     float minK = 0.1;
 158     float maxK = 10000;
 159     float k = 5;
 160
 161     for (int i = 0; i < 15; ++i) {
 162         // A high value for k will more quickly asymptotically approach a slope of 0.
 163         float slope = slopeAt(x, k);
 164
 165         if (slope < desiredSlope) {
 166             // k is too high.
 167             maxK = k;
 168         } else {
 169             // k is too low.
 170             minK = k;
 171         }
 172
 173         // Re-calculate based on geometric mean.
 174         k = sqrtf(minK * maxK);
 175     }
 176
 177     return k;
 178 }
 179
 180 float DynamicsCompressorKernel::updateStaticCurveParameters(float dbThreshold, float dbKnee, float ratio)
 181 {
 182     if (dbThreshold != m_dbThreshold || dbKnee != m_dbKnee || ratio != m_ratio) {
 183         // Threshold and knee.
 184         m_dbThreshold = dbThreshold;
 185         m_linearThreshold = decibelsToLinear(dbThreshold);
 186         m_dbKnee = dbKnee;
 187
 188         // Compute knee parameters.
 189         m_ratio = ratio;
 190         m_slope = 1 / m_ratio;
 191
 192         float k = kAtSlope(1 / m_ratio);
 193
 194         m_kneeThresholdDb = dbThreshold + dbKnee;
 195         m_kneeThreshold = decibelsToLinear(m_kneeThresholdDb);
 196
 197         m_ykneeThresholdDb = linearToDecibels(kneeCurve(m_kneeThreshold, k));
 198
 199         m_K = k;
 200     }
 201     return m_K;
 202 }
 203
 204 void DynamicsCompressorKernel::process(const float* sourceChannels[],
 205                                        float* destinationChannels[],
 206                                        unsigned numberOfChannels,
 207                                        unsigned framesToProcess,
 208
 209                                        float dbThreshold,
 210                                        float dbKnee,
 211                                        float ratio,
 212                                        float attackTime,
 213                                        float releaseTime,
 214                                        float preDelayTime,
 215                                        float dbPostGain,
 216                                        float effectBlend, /* equal power crossfade */
 217
 218                                        float releaseZone1,
 219                                        float releaseZone2,
 220                                        float releaseZone3,
 221                                        float releaseZone4
 222                                        )
 223 {
 224     ASSERT(m_preDelayBuffers.size() == numberOfChannels);
 225
 226     float sampleRate = this->sampleRate();
 227
 228     float dryMix = 1 - effectBlend;
 229     float wetMix = effectBlend;
 230
 231     float k = updateStaticCurveParameters(dbThreshold, dbKnee, ratio);
 232
 233     // Makeup gain.
 234     float fullRangeGain = saturate(1, k);
 235     float fullRangeMakeupGain = 1 / fullRangeGain;
 236
 237     // Empirical/perceptual tuning.
 238     fullRangeMakeupGain = powf(fullRangeMakeupGain, 0.6f);
 239
 240     float masterLinearGain = decibelsToLinear(dbPostGain) * fullRangeMakeupGain;
 241
 242     // Attack parameters.
 243     attackTime = max(0.001f, attackTime);
 244     float attackFrames = attackTime * sampleRate;
 245
 246     // Release parameters.
 247     float releaseFrames = sampleRate * releaseTime;
 248
 249     // Detector release time.
 250     float satReleaseTime = 0.0025f;
 251     float satReleaseFrames = satReleaseTime * sampleRate;
 252
 253     // Create a smooth function which passes through four points.
 254
 255     // Polynomial of the form
 256     // y = a + b*x + c*x^2 + d*x^3 + e*x^4;
 257
 258     float y1 = releaseFrames * releaseZone1;
 259     float y2 = releaseFrames * releaseZone2;
 260     float y3 = releaseFrames * releaseZone3;
 261     float y4 = releaseFrames * releaseZone4;
 262
 263     // All of these coefficients were derived for 4th order polynomial curve fitting where the y values
 264     // match the evenly spaced x values as follows: (y1 : x == 0, y2 : x == 1, y3 : x == 2, y4 : x == 3)
 265     float kA = 0.9999999999999998f*y1 + 1.8432219684323923e-16f*y2 - 1.9373394351676423e-16f*y3 + 8.824516011816245e-18f*y4;
 266     float kB = -1.5788320352845888f*y1 + 2.3305837032074286f*y2 - 0.9141194204840429f*y3 + 0.1623677525612032f*y4;
 267     float kC = 0.5334142869106424f*y1 - 1.272736789213631f*y2 + 0.9258856042207512f*y3 - 0.18656310191776226f*y4;
 268     float kD = 0.08783463138207234f*y1 - 0.1694162967925622f*y2 + 0.08588057951595272f*y3 - 0.00429891410546283f*y4;
 269     float kE = -0.042416883008123074f*y1 + 0.1115693827987602f*y2 - 0.09764676325265872f*y3 + 0.028494263462021576f*y4;
 270
 271     // x ranges from 0 -> 3       0    1    2   3
 272     //                           -15  -10  -5   0db
 273
 274     // y calculates adaptive release frames depending on the amount of compression.
 275
 276     setPreDelayTime(preDelayTime);
 277
 278     const int nDivisionFrames = 32;
 279
 280     const int nDivisions = framesToProcess / nDivisionFrames;
 281
 282     unsigned frameIndex = 0;
 283     for (int i = 0; i < nDivisions; ++i) {
 284         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 285         // Calculate desired gain
 286         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 287
 288         // Fix gremlins.
 289         if (std::isnan(m_detectorAverage))
 290             m_detectorAverage = 1;
 291         if (std::isinf(m_detectorAverage))
 292             m_detectorAverage = 1;
 293
 294         float desiredGain = m_detectorAverage;
 295
 296         // Pre-warp so we get desiredGain after sin() warp below.
 297         float scaledDesiredGain = asinf(desiredGain) / (0.5f * piFloat);
 298
 299         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 300         // Deal with envelopes
 301         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 302
 303         // envelopeRate is the rate we slew from current compressor level to the desired level.
 304         // The exact rate depends on if we're attacking or releasing and by how much.
 305         float envelopeRate;
 306
 307         bool isReleasing = scaledDesiredGain > m_compressorGain;
 308
 309         // compressionDiffDb is the difference between current compression level and the desired level.
 310         float compressionDiffDb = linearToDecibels(m_compressorGain / scaledDesiredGain);
 311
 312         if (isReleasing) {
 313             // Release mode - compressionDiffDb should be negative dB
 314             m_maxAttackCompressionDiffDb = -1;
 315
 316             // Fix gremlins.
 317             if (std::isnan(compressionDiffDb))
 318                 compressionDiffDb = -1;
 319             if (std::isinf(compressionDiffDb))
 320                 compressionDiffDb = -1;
 321
 322             // Adaptive release - higher compression (lower compressionDiffDb)  releases faster.
 323
 324             // Contain within range: -12 -> 0 then scale to go from 0 -> 3
 325             float x = compressionDiffDb;
 326             x = max(-12.0f, x);
 327             x = min(0.0f, x);
 328             x = 0.25f * (x + 12);
 329
 330             // Compute adaptive release curve using 4th order polynomial.
 331             // Normal values for the polynomial coefficients would create a monotonically increasing function.
 332             float x2 = x * x;
 333             float x3 = x2 * x;
 334             float x4 = x2 * x2;
 335             float releaseFrames = kA + kB * x + kC * x2 + kD * x3 + kE * x4;
 336
 337 #define kSpacingDb 5
 338             float dbPerFrame = kSpacingDb / releaseFrames;
 339
 340             envelopeRate = decibelsToLinear(dbPerFrame);
 341         } else {
 342             // Attack mode - compressionDiffDb should be positive dB
 343
 344             // Fix gremlins.
 345             if (std::isnan(compressionDiffDb))
 346                 compressionDiffDb = 1;
 347             if (std::isinf(compressionDiffDb))
 348                 compressionDiffDb = 1;
 349
 350             // As long as we're still in attack mode, use a rate based off
 351             // the largest compressionDiffDb we've encountered so far.
 352             if (m_maxAttackCompressionDiffDb == -1 || m_maxAttackCompressionDiffDb < compressionDiffDb)
 353                 m_maxAttackCompressionDiffDb = compressionDiffDb;
 354
 355             float effAttenDiffDb = max(0.5f, m_maxAttackCompressionDiffDb);
 356
 357             float x = 0.25f / effAttenDiffDb;
 358             envelopeRate = 1 - powf(x, 1 / attackFrames);
 359         }
 360
 361         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 362         // Inner loop - calculate shaped power average - apply compression.
 363         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 364
 365         {
 366             int preDelayReadIndex = m_preDelayReadIndex;
 367             int preDelayWriteIndex = m_preDelayWriteIndex;
 368             float detectorAverage = m_detectorAverage;
 369             float compressorGain = m_compressorGain;
 370
 371             int loopFrames = nDivisionFrames;
 372             while (loopFrames--) {
 373                 float compressorInput = 0;
 374
 375                 // Predelay signal, computing compression amount from un-delayed version.
 376                 for (unsigned i = 0; i < numberOfChannels; ++i) {
 377                     float* delayBuffer = m_preDelayBuffers[i]->data();
 378                     float undelayedSource = sourceChannels[i][frameIndex];
 379                     delayBuffer[preDelayWriteIndex] = undelayedSource;
 380
 381                     float absUndelayedSource = undelayedSource > 0 ? undelayedSource : -undelayedSource;
 382                     if (compressorInput < absUndelayedSource)
 383                         compressorInput = absUndelayedSource;
 384                 }
 385
 386                 // Calculate shaped power on undelayed input.
 387
 388                 float scaledInput = compressorInput;
 389                 float absInput = scaledInput > 0 ? scaledInput : -scaledInput;
 390
 391                 // Put through shaping curve.
 392                 // This is linear up to the threshold, then enters a "knee" portion followed by the "ratio" portion.
 393                 // The transition from the threshold to the knee is smooth (1st derivative matched).
 394                 // The transition from the knee to the ratio portion is smooth (1st derivative matched).
 395                 float shapedInput = saturate(absInput, k);
 396
 397                 float attenuation = absInput <= 0.0001f ? 1 : shapedInput / absInput;
 398
 399                 float attenuationDb = -linearToDecibels(attenuation);
 400                 attenuationDb = max(2.0f, attenuationDb);
 401
 402                 float dbPerFrame = attenuationDb / satReleaseFrames;
 403
 404                 float satReleaseRate = decibelsToLinear(dbPerFrame) - 1;
 405
 406                 bool isRelease = (attenuation > detectorAverage);
 407                 float rate = isRelease ? satReleaseRate : 1;
 408
 409                 detectorAverage += (attenuation - detectorAverage) * rate;
 410                 detectorAverage = min(1.0f, detectorAverage);
 411
 412                 // Fix gremlins.
 413                 if (std::isnan(detectorAverage))
 414                     detectorAverage = 1;
 415                 if (std::isinf(detectorAverage))
 416                     detectorAverage = 1;
 417
 418                 // Exponential approach to desired gain.
 419                 if (envelopeRate < 1) {
 420                     // Attack - reduce gain to desired.
 421                     compressorGain += (scaledDesiredGain - compressorGain) * envelopeRate;
 422                 } else {
 423                     // Release - exponentially increase gain to 1.0
 424                     compressorGain *= envelopeRate;
 425                     compressorGain = min(1.0f, compressorGain);
 426                 }
 427
 428                 // Warp pre-compression gain to smooth out sharp exponential transition points.
 429                 float postWarpCompressorGain = sinf(0.5f * piFloat * compressorGain);
 430
 431                 // Calculate total gain using master gain and effect blend.
 432                 float totalGain = dryMix + wetMix * masterLinearGain * postWarpCompressorGain;
 433
 434                 // Calculate metering.
 435                 float dbRealGain = 20 * log10(postWarpCompressorGain);
 436                 if (dbRealGain < m_meteringGain)
 437                     m_meteringGain = dbRealGain;
 438                 else
 439                     m_meteringGain += (dbRealGain - m_meteringGain) * m_meteringReleaseK;
 440
 441                 // Apply final gain.
 442                 for (unsigned i = 0; i < numberOfChannels; ++i) {
 443                     float* delayBuffer = m_preDelayBuffers[i]->data();
 444                     destinationChannels[i][frameIndex] = delayBuffer[preDelayReadIndex] * totalGain;
 445                 }
 446
 447                 frameIndex++;
 448                 preDelayReadIndex = (preDelayReadIndex + 1) & MaxPreDelayFramesMask;
 449                 preDelayWriteIndex = (preDelayWriteIndex + 1) & MaxPreDelayFramesMask;
 450             }
 451
 452             // Locals back to member variables.
 453             m_preDelayReadIndex = preDelayReadIndex;
 454             m_preDelayWriteIndex = preDelayWriteIndex;
 455             m_detectorAverage = DenormalDisabler::flushDenormalFloatToZero(detectorAverage);
 456             m_compressorGain = DenormalDisabler::flushDenormalFloatToZero(compressorGain);
 457         }
 458     }
 459 }
 460
 461 void DynamicsCompressorKernel::reset()
 462 {
 463     m_detectorAverage = 0;
 464     m_compressorGain = 1;
 465     m_meteringGain = 1;
 466
 467     // Predelay section.
 468     for (unsigned i = 0; i < m_preDelayBuffers.size(); ++i)
 469         m_preDelayBuffers[i]->zero();
 470
 471     m_preDelayReadIndex = 0;
 472     m_preDelayWriteIndex = DefaultPreDelayFrames;
 473
 474     m_maxAttackCompressionDiffDb = -1; // uninitialized state
 475 }
 476
 477 } // namespace WebCore
 478
 479 #endif // ENABLE(WEB_AUDIO)