5 #include "stereocompressor.h"
9 // Implement a less accurate but faster pow(x, y). We use the standard identity
11 // x^y = exp(y * ln(x))
15 // x in 1..(1/threshold)
18 // Assume threshold goes from 0 to -40 dB. That means 1/threshold = 100,
19 // so input to ln(x) can be 1..100. Worst case for end accuracy is y=-1.
20 // To get a good minimax approximation (not the least wrt. continuity
21 // at x=1), I had to make a piecewise linear function for the two ranges:
24 // f1 := minimax(ln, 1..6, [3, 3], x -> 1/x, 'maxerror');
25 // f2 := minimax(ln, 6..100, [3, 3], x -> 1/x, 'maxerror');
26 // f := x -> piecewise(x < 6, f1(x), f2(x));
28 // (Continuity: Error is down to the 1e-6 range for x=1, difference between
29 // f1 and f2 range at the crossover point is in the 1e-5 range. The cutoff
30 // point at x=6 is chosen to get maxerror pretty close between f1 and f2.)
32 // Maximum output of ln(x) here is of course ln(100) ~= 4.605. So we can find
33 // an approximation for exp over the range -4.605..0, where we care mostly
34 // about the relative error:
36 // g := minimax(exp, -ln(100)..0, [3, 3], x -> 1/exp(x), 'maxerror');
38 // We can find the worst-case error in dB from this through a simple plot:
40 // dbdiff := (x, y) -> abs(20 * log10(x / y));
41 // plot(dbdiff(g(-f(x)), 1/x), x=1..100);
43 // which readily shows the error never to be above ~0.001 dB or so
44 // (actually 0.00119 dB, for the case of x=100). y=-1 remains the worst case,
47 // If we cared even more about speed, we could probably fuse y into
48 // the coefficients for ln_nom and postgain into the coefficients for ln_den.
49 // But if so, we should probably rather just SIMD the entire thing instead.
50 inline float fastpow(float x, float y)
54 ln_nom = -0.059237648 + (-0.0165117771 + (0.06818859075 + 0.007560968243 * x) * x) * x;
55 ln_den = 0.0202509098 + (0.08419174188 + (0.03647189417 + 0.001642577975 * x) * x) * x;
57 ln_nom = -0.005430534 + (0.00633589178 + (0.0006319155549 + 0.4789541675e-5 * x) * x) * x;
58 ln_den = 0.0064785099 + (0.003219629109 + (0.0001531823694 + 0.6884656640e-6 * x) * x) * x;
60 float v = y * ln_nom / ln_den;
61 float exp_nom = 0.2195097621 + (0.08546059868 + (0.01208501759 + 0.0006173448113 * v) * v) * v;
62 float exp_den = 0.2194980791 + (-0.1343051968 + (0.03556072737 - 0.006174398513 * v) * v) * v;
63 return exp_nom / exp_den;
66 inline float compressor_knee(float x, float threshold, float inv_threshold, float inv_ratio_minus_one, float postgain)
68 assert(inv_ratio_minus_one <= 0.0f);
69 if (x > threshold && inv_ratio_minus_one < 0.0f) {
70 return postgain * fastpow(x * inv_threshold, inv_ratio_minus_one);
78 void StereoCompressor::process(float *buf, size_t num_samples, float threshold, float ratio,
79 float attack_time, float release_time, float makeup_gain)
81 float attack_increment = float(pow(2.0f, 1.0f / (attack_time * sample_rate + 1)));
82 if (attack_time == 0.0f) attack_increment = 100000; // For instant attack reaction.
84 const float release_increment = float(pow(2.0f, -1.0f / (release_time * sample_rate + 1)));
85 const float peak_increment = float(pow(2.0f, -1.0f / (0.003f * sample_rate + 1)));
87 float inv_ratio_minus_one = 1.0f / ratio - 1.0f;
88 if (ratio > 63) inv_ratio_minus_one = -1.0f; // Infinite ratio.
89 float inv_threshold = 1.0f / threshold;
91 float *left_ptr = buf;
92 float *right_ptr = buf + 1;
94 float peak_level = this->peak_level;
95 float compr_level = this->compr_level;
97 for (size_t i = 0; i < num_samples; ++i) {
98 if (fabs(*left_ptr) > peak_level) peak_level = float(fabs(*left_ptr));
99 if (fabs(*right_ptr) > peak_level) peak_level = float(fabs(*right_ptr));
101 if (peak_level > compr_level) {
102 compr_level = std::min(compr_level * attack_increment, peak_level);
104 compr_level = std::max(compr_level * release_increment, 0.0001f);
107 float scalefactor_with_gain = compressor_knee(compr_level, threshold, inv_threshold, inv_ratio_minus_one, makeup_gain);
109 *left_ptr *= scalefactor_with_gain;
112 *right_ptr *= scalefactor_with_gain;
115 peak_level = std::max(peak_level * peak_increment, 0.0001f);
118 // Store attenuation level for debug/visualization.
119 scalefactor = compressor_knee(compr_level, threshold, inv_threshold, inv_ratio_minus_one, 1.0f);
121 this->peak_level = peak_level;
122 this->compr_level = compr_level;