git.sesse.net Git - nageru/blob - nageru/stereocompressor.cpp

   1 #include "stereocompressor.h"
   2
   3 #include <assert.h>
   4 #include <algorithm>
   5 #include <math.h>
   6 #include <stddef.h>
   7
   8 using namespace std;
   9
  10 namespace {
  11
  12 // Implement a less accurate but faster pow(x, y). We use the standard identity
  13 //
  14 //    x^y = exp(y * ln(x))
  15 //
  16 // with the ranges:
  17 //
  18 //    x in 1..(1/threshold)
  19 //    y in -1..0
  20 //
  21 // Assume threshold goes from 0 to -40 dB. That means 1/threshold = 100,
  22 // so input to ln(x) can be 1..100. Worst case for end accuracy is y=-1.
  23 // To get a good minimax approximation (not the least wrt. continuity
  24 // at x=1), I had to make a piecewise linear function for the two ranges:
  25 //
  26 //   with(numapprox):
  27 //   f1 := minimax(ln, 1..6, [3, 3], x -> 1/x, 'maxerror');
  28 //   f2 := minimax(ln, 6..100, [3, 3], x -> 1/x, 'maxerror');
  29 //   f := x -> piecewise(x < 6, f1(x), f2(x));
  30 //
  31 // (Continuity: Error is down to the 1e-6 range for x=1, difference between
  32 // f1 and f2 range at the crossover point is in the 1e-5 range. The cutoff
  33 // point at x=6 is chosen to get maxerror pretty close between f1 and f2.)
  34 //
  35 // Maximum output of ln(x) here is of course ln(100) ~= 4.605. So we can find
  36 // an approximation for exp over the range -4.605..0, where we care mostly
  37 // about the relative error:
  38 //
  39 //   g := minimax(exp, -ln(100)..0, [3, 3], x -> 1/exp(x), 'maxerror');
  40 //
  41 // We can find the worst-case error in dB from this through a simple plot:
  42 //
  43 //   dbdiff := (x, y) -> abs(20 * log10(x / y));
  44 //   plot(dbdiff(g(-f(x)), 1/x), x=1..100);
  45 //
  46 // which readily shows the error never to be above ~0.001 dB or so
  47 // (actually 0.00119 dB, for the case of x=100). y=-1 remains the worst case,
  48 // it would seem.
  49 //
  50 // If we cared even more about speed, we could probably fuse y into
  51 // the coefficients for ln_nom and postgain into the coefficients for ln_den.
  52 // But if so, we should probably rather just SIMD the entire thing instead.
  53 inline float fastpow(float x, float y)
  54 {
  55         float ln_nom, ln_den;
  56         if (x < 6.0f) {
  57                 ln_nom = -0.059237648f + (-0.0165117771f + (0.06818859075f + 0.007560968243f * x) * x) * x;
  58                 ln_den = 0.0202509098f + (0.08419174188f + (0.03647189417f + 0.001642577975f * x) * x) * x;
  59         } else {
  60                 ln_nom = -0.005430534f + (0.00633589178f + (0.0006319155549f + 0.4789541675e-5f * x) * x) * x;
  61                 ln_den = 0.0064785099f + (0.003219629109f + (0.0001531823694f + 0.6884656640e-6f * x) * x) * x;
  62         }
  63         float v = y * ln_nom / ln_den;
  64         float exp_nom = 0.2195097621f + (0.08546059868f + (0.01208501759f + 0.0006173448113f * v) * v) * v;
  65         float exp_den = 0.2194980791f + (-0.1343051968f + (0.03556072737f - 0.006174398513f * v) * v) * v;
  66         return exp_nom / exp_den;
  67 }
  68
  69 inline float compressor_knee(float x, float threshold, float inv_threshold, float inv_ratio_minus_one, float postgain)
  70 {
  71         assert(inv_ratio_minus_one <= 0.0f);
  72         if (x > threshold) {
  73                 return postgain * fastpow(x * inv_threshold, inv_ratio_minus_one);
  74         } else {
  75                 return postgain;
  76         }
  77 }
  78
  79 }  // namespace
  80
  81 void StereoCompressor::process(float *buf, size_t num_samples, float threshold, float ratio,
  82             float attack_time, float release_time, float makeup_gain)
  83 {
  84         float attack_increment = float(pow(2.0f, 1.0f / (attack_time * sample_rate + 1)));
  85         if (attack_time == 0.0f) attack_increment = 100000;  // For instant attack reaction.
  86
  87         const float release_increment = float(pow(2.0f, -1.0f / (release_time * sample_rate + 1)));
  88         const float peak_increment = float(pow(2.0f, -1.0f / (0.003f * sample_rate + 1)));
  89
  90         float inv_ratio_minus_one = 1.0f / ratio - 1.0f;
  91         if (ratio > 63) inv_ratio_minus_one = -1.0f;  // Infinite ratio.
  92         float inv_threshold = 1.0f / threshold;
  93
  94         float *left_ptr = buf;
  95         float *right_ptr = buf + 1;
  96
  97         if (inv_ratio_minus_one >= 0.0) {
  98                 for (size_t i = 0; i < num_samples; ++i) {
  99                         *left_ptr *= makeup_gain;
 100                         left_ptr += 2;
 101
 102                         *right_ptr *= makeup_gain;
 103                         right_ptr += 2;
 104                 }
 105                 return;
 106         }
 107
 108         float peak_level = this->peak_level;
 109         float compr_level = this->compr_level;
 110
 111         for (size_t i = 0; i < num_samples; ++i) {
 112                 if (fabs(*left_ptr) > peak_level) peak_level = float(fabs(*left_ptr));
 113                 if (fabs(*right_ptr) > peak_level) peak_level = float(fabs(*right_ptr));
 114
 115                 if (peak_level > compr_level) {
 116                         compr_level = min(compr_level * attack_increment, peak_level);
 117                 } else {
 118                         compr_level = max(compr_level * release_increment, 0.0001f);
 119                 }
 120
 121                 float scalefactor_with_gain = compressor_knee(compr_level, threshold, inv_threshold, inv_ratio_minus_one, makeup_gain);
 122
 123                 *left_ptr *= scalefactor_with_gain;
 124                 left_ptr += 2;
 125
 126                 *right_ptr *= scalefactor_with_gain;
 127                 right_ptr += 2;
 128
 129                 peak_level = max(peak_level * peak_increment, 0.0001f);
 130         }
 131
 132         // Store attenuation level for debug/visualization.
 133         scalefactor = compressor_knee(compr_level, threshold, inv_threshold, inv_ratio_minus_one, 1.0f);
 134
 135         this->peak_level = peak_level;
 136         this->compr_level = compr_level;
 137 }
 138