Home | History | Annotate | Download | only in audio
      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 "core/platform/audio/DynamicsCompressorKernel.h"
     34 
     35 #include <algorithm>
     36 #include "core/platform/audio/AudioUtilities.h"
     37 #include "core/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(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)
    480