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