|
| 1 | +/* From https://github.com/chipaudette/OpenAudio_ArduinoLibrary */ |
| 2 | + |
| 3 | +/* |
| 4 | + AudioEffectCompressor |
| 5 | +
|
| 6 | + Created: Chip Audette, Dec 2016 - Jan 2017 |
| 7 | + Purpose; Apply dynamic range compression to the audio stream. |
| 8 | + Assumes floating-point data. |
| 9 | +
|
| 10 | + This processes a single stream fo audio data (ie, it is mono) |
| 11 | +
|
| 12 | + MIT License. use at your own risk. |
| 13 | +*/ |
| 14 | + |
| 15 | +#include <circle/logger.h> |
| 16 | +#include <cstdlib> |
| 17 | +#include "effect_compressor.h" |
| 18 | + |
| 19 | +LOGMODULE ("compressor"); |
| 20 | + |
| 21 | +Compressor::Compressor(const float32_t sample_rate_Hz) { |
| 22 | + //setDefaultValues(AUDIO_SAMPLE_RATE); resetStates(); |
| 23 | + setDefaultValues(sample_rate_Hz); |
| 24 | + resetStates(); |
| 25 | +} |
| 26 | + |
| 27 | +void Compressor::setDefaultValues(const float32_t sample_rate_Hz) { |
| 28 | + setThresh_dBFS(-20.0f); //set the default value for the threshold for compression |
| 29 | + setCompressionRatio(5.0f); //set the default copression ratio |
| 30 | + setAttack_sec(0.005f, sample_rate_Hz); //default to this value |
| 31 | + setRelease_sec(0.200f, sample_rate_Hz); //default to this value |
| 32 | + setHPFilterCoeff(); enableHPFilter(true); //enable the HP filter to remove any DC offset from the audio |
| 33 | +} |
| 34 | + |
| 35 | +//Compute the instantaneous desired gain, including the compression ratio and |
| 36 | +//threshold for where the comrpession kicks in |
| 37 | +void Compressor::calcInstantaneousTargetGain(float32_t *audio_level_dB_block, float32_t *inst_targ_gain_dB_block, uint16_t len) |
| 38 | +{ |
| 39 | + // how much are we above the compression threshold? |
| 40 | + float32_t above_thresh_dB_block[len]; |
| 41 | + |
| 42 | + //arm_copy_f32(zeroblock_f32,above_thresh_dB_block,len); |
| 43 | + |
| 44 | + arm_offset_f32(audio_level_dB_block, //CMSIS DSP for "add a constant value to all elements" |
| 45 | + -thresh_dBFS, //this is the value to be added |
| 46 | + above_thresh_dB_block, //this is the output |
| 47 | + len); |
| 48 | + |
| 49 | + // scale by the compression ratio...this is what the output level should be (this is our target level) |
| 50 | + arm_scale_f32(above_thresh_dB_block, //CMSIS DSP for "multiply all elements by a constant value" |
| 51 | + 1.0f / comp_ratio, //this is the value to be multiplied |
| 52 | + inst_targ_gain_dB_block, //this is the output |
| 53 | + len); |
| 54 | + |
| 55 | + // compute the instantaneous gain...which is the difference between the target level and the original level |
| 56 | + arm_sub_f32(inst_targ_gain_dB_block, //CMSIS DSP for "subtract two vectors element-by-element" |
| 57 | + above_thresh_dB_block, //this is the vector to be subtracted |
| 58 | + inst_targ_gain_dB_block, //this is the output |
| 59 | + len); |
| 60 | + |
| 61 | + // limit the target gain to attenuation only (this part of the compressor should not make things louder!) |
| 62 | + for (uint16_t i=0; i < len; i++) { |
| 63 | + if (inst_targ_gain_dB_block[i] > 0.0f) inst_targ_gain_dB_block[i] = 0.0f; |
| 64 | + } |
| 65 | + |
| 66 | + return; //output is passed through inst_targ_gain_dB_block |
| 67 | +} |
| 68 | + |
| 69 | +//this method applies the "attack" and "release" constants to smooth the |
| 70 | +//target gain level through time. |
| 71 | +void Compressor::calcSmoothedGain_dB(float32_t *inst_targ_gain_dB_block, float32_t *gain_dB_block, uint16_t len) |
| 72 | +{ |
| 73 | + float32_t gain_dB; |
| 74 | + float32_t one_minus_attack_const = 1.0f - attack_const; |
| 75 | + float32_t one_minus_release_const = 1.0f - release_const; |
| 76 | + for (uint16_t i = 0; i < len; i++) { |
| 77 | + gain_dB = inst_targ_gain_dB_block[i]; |
| 78 | + |
| 79 | + //smooth the gain using the attack or release constants |
| 80 | + if (gain_dB < prev_gain_dB) { //are we in the attack phase? |
| 81 | + gain_dB_block[i] = attack_const*prev_gain_dB + one_minus_attack_const*gain_dB; |
| 82 | + } else { //or, we're in the release phase |
| 83 | + gain_dB_block[i] = release_const*prev_gain_dB + one_minus_release_const*gain_dB; |
| 84 | + } |
| 85 | + |
| 86 | + //save value for the next time through this loop |
| 87 | + prev_gain_dB = gain_dB_block[i]; |
| 88 | + } |
| 89 | + |
| 90 | + return; //the output here is gain_block |
| 91 | +} |
| 92 | + |
| 93 | +// Here's the method that estimates the level of the audio (in dB) |
| 94 | +// It squares the signal and low-pass filters to get a time-averaged |
| 95 | +// signal power. It then |
| 96 | +void Compressor::calcAudioLevel_dB(float32_t *wav_block, float32_t *level_dB_block, uint16_t len) { |
| 97 | + |
| 98 | + // calculate the instantaneous signal power (square the signal) |
| 99 | + float32_t wav_pow_block[len]; |
| 100 | + |
| 101 | + //arm_copy_f32(zeroblock_f32,wav_pow_block,len); |
| 102 | + |
| 103 | + arm_mult_f32(wav_block, wav_block, wav_pow_block, len); |
| 104 | + |
| 105 | + // low-pass filter and convert to dB |
| 106 | + float32_t c1 = level_lp_const, c2 = 1.0f - c1; //prepare constants |
| 107 | + for (uint16_t i = 0; i < len; i++) { |
| 108 | + // first-order low-pass filter to get a running estimate of the average power |
| 109 | + wav_pow_block[i] = c1*prev_level_lp_pow + c2*wav_pow_block[i]; |
| 110 | + |
| 111 | + // save the state of the first-order low-pass filter |
| 112 | + prev_level_lp_pow = wav_pow_block[i]; |
| 113 | + |
| 114 | + //now convert the signal power to dB (but not yet multiplied by 10.0) |
| 115 | + level_dB_block[i] = log10f_approx(wav_pow_block[i]); |
| 116 | + } |
| 117 | + |
| 118 | + //limit the amount that the state of the smoothing filter can go toward negative infinity |
| 119 | + if (prev_level_lp_pow < (1.0E-13)) prev_level_lp_pow = 1.0E-13; //never go less than -130 dBFS |
| 120 | + |
| 121 | + //scale the wav_pow_block by 10.0 to complete the conversion to dB |
| 122 | + arm_scale_f32(level_dB_block, 10.0f, level_dB_block, len); //use ARM DSP for speed! |
| 123 | + |
| 124 | + return; //output is passed through level_dB_block |
| 125 | + } |
| 126 | + |
| 127 | + //This method computes the desired gain from the compressor, given an estimate |
| 128 | + //of the signal level (in dB) |
| 129 | +void Compressor::calcGain(float32_t *audio_level_dB_block, float32_t *gain_block,uint16_t len) |
| 130 | +{ |
| 131 | + //first, calculate the instantaneous target gain based on the compression ratio |
| 132 | + float32_t inst_targ_gain_dB_block[len]; |
| 133 | + //arm_copy_f32(zeroblock_f32,inst_targ_gain_dB_block,len); |
| 134 | + |
| 135 | + calcInstantaneousTargetGain(audio_level_dB_block, inst_targ_gain_dB_block,len); |
| 136 | + |
| 137 | + //second, smooth in time (attack and release) by stepping through each sample |
| 138 | + float32_t gain_dB_block[len]; |
| 139 | + //arm_copy_f32(zeroblock_f32,gain_dB_block,len); |
| 140 | + |
| 141 | + calcSmoothedGain_dB(inst_targ_gain_dB_block,gain_dB_block, len); |
| 142 | + |
| 143 | + //finally, convert from dB to linear gain: gain = 10^(gain_dB/20); (ie this takes care of the sqrt, too!) |
| 144 | + arm_scale_f32(gain_dB_block, 1.0f/20.0f, gain_dB_block, len); //divide by 20 |
| 145 | + for (uint16_t i = 0; i < len; i++) gain_block[i] = pow10f(gain_dB_block[i]); //do the 10^(x) |
| 146 | + |
| 147 | + return; //output is passed through gain_block |
| 148 | +} |
| 149 | + |
| 150 | +//here's the method that does all the work |
| 151 | +void Compressor::doCompression(float32_t *audio_block, uint16_t len) { |
| 152 | + //Serial.println("AudioEffectGain_F32: updating."); //for debugging. |
| 153 | + if (!audio_block) { |
| 154 | + LOGERR("No audio_block available for Compressor!"); |
| 155 | + return; |
| 156 | + } |
| 157 | + |
| 158 | + //apply a high-pass filter to get rid of the DC offset |
| 159 | + if (use_HP_prefilter) |
| 160 | + arm_biquad_cascade_df1_f32(&hp_filt_struct, audio_block, audio_block, len); |
| 161 | + |
| 162 | + //apply the pre-gain...a negative gain value will disable |
| 163 | + if (pre_gain > 0.0f) |
| 164 | + arm_scale_f32(audio_block, pre_gain, audio_block, len); //use ARM DSP for speed! |
| 165 | + |
| 166 | + //calculate the level of the audio (ie, calculate a smoothed version of the signal power) |
| 167 | + float32_t audio_level_dB_block[len]; |
| 168 | + |
| 169 | + //arm_copy_f32(zeroblock_f32,audio_level_dB_block,len); |
| 170 | + |
| 171 | + calcAudioLevel_dB(audio_block, audio_level_dB_block, len); //returns through audio_level_dB_block |
| 172 | + |
| 173 | + //compute the desired gain based on the observed audio level |
| 174 | + float32_t gain_block[len]; |
| 175 | + |
| 176 | + //arm_copy_f32(zeroblock_f32,gain_block,len); |
| 177 | + |
| 178 | + calcGain(audio_level_dB_block, gain_block, len); //returns through gain_block |
| 179 | + |
| 180 | + //apply the desired gain...store the processed audio back into audio_block |
| 181 | + arm_mult_f32(audio_block, gain_block, audio_block, len); |
| 182 | +} |
| 183 | + |
| 184 | +//methods to set parameters of this module |
| 185 | +void Compressor::resetStates(void) |
| 186 | +{ |
| 187 | + prev_level_lp_pow = 1.0f; |
| 188 | + prev_gain_dB = 0.0f; |
| 189 | + |
| 190 | + //initialize the HP filter. (This also resets the filter states,) |
| 191 | + arm_biquad_cascade_df1_init_f32(&hp_filt_struct, hp_nstages, hp_coeff, hp_state); |
| 192 | +} |
| 193 | + |
| 194 | +void Compressor::setPreGain(float32_t g) |
| 195 | +{ |
| 196 | + pre_gain = g; |
| 197 | +} |
| 198 | + |
| 199 | +void Compressor::setPreGain_dB(float32_t gain_dB) |
| 200 | +{ |
| 201 | + setPreGain(pow(10.0, gain_dB / 20.0)); |
| 202 | +} |
| 203 | + |
| 204 | +void Compressor::setCompressionRatio(float32_t cr) |
| 205 | +{ |
| 206 | + comp_ratio = max(0.001f, cr); //limit to positive values |
| 207 | + updateThresholdAndCompRatioConstants(); |
| 208 | +} |
| 209 | + |
| 210 | +void Compressor::setAttack_sec(float32_t a, float32_t fs_Hz) |
| 211 | +{ |
| 212 | + attack_sec = a; |
| 213 | + attack_const = expf(-1.0f / (attack_sec * fs_Hz)); //expf() is much faster than exp() |
| 214 | + |
| 215 | + //also update the time constant for the envelope extraction |
| 216 | + setLevelTimeConst_sec(min(attack_sec,release_sec) / 5.0, fs_Hz); //make the level time-constant one-fifth the gain time constants |
| 217 | +} |
| 218 | + |
| 219 | +void Compressor::setRelease_sec(float32_t r, float32_t fs_Hz) |
| 220 | +{ |
| 221 | + release_sec = r; |
| 222 | + release_const = expf(-1.0f / (release_sec * fs_Hz)); //expf() is much faster than exp() |
| 223 | + |
| 224 | + //also update the time constant for the envelope extraction |
| 225 | + setLevelTimeConst_sec(min(attack_sec,release_sec) / 5.0, fs_Hz); //make the level time-constant one-fifth the gain time constants |
| 226 | +} |
| 227 | + |
| 228 | +void Compressor::setLevelTimeConst_sec(float32_t t_sec, float32_t fs_Hz) |
| 229 | +{ |
| 230 | + const float32_t min_t_sec = 0.002f; //this is the minimum allowed value |
| 231 | + level_lp_sec = max(min_t_sec,t_sec); |
| 232 | + level_lp_const = expf(-1.0f / (level_lp_sec * fs_Hz)); //expf() is much faster than exp() |
| 233 | +} |
| 234 | + |
| 235 | +void Compressor::setThresh_dBFS(float32_t val) |
| 236 | +{ |
| 237 | + thresh_dBFS = val; |
| 238 | + setThreshPow(pow(10.0, thresh_dBFS / 10.0)); |
| 239 | +} |
| 240 | + |
| 241 | +void Compressor::enableHPFilter(boolean flag) |
| 242 | +{ |
| 243 | + use_HP_prefilter = flag; |
| 244 | +} |
| 245 | + |
| 246 | +void Compressor::setHPFilterCoeff_N2IIR_Matlab(float32_t b[], float32_t a[]) |
| 247 | +{ |
| 248 | + //https://www.keil.com/pack/doc/CMSIS/DSP/html/group__BiquadCascadeDF1.html#ga8e73b69a788e681a61bccc8959d823c5 |
| 249 | + //Use matlab to compute the coeff for HP at 20Hz: [b,a]=butter(2,20/(44100/2),'high'); %assumes fs_Hz = 44100 |
| 250 | + hp_coeff[0] = b[0]; hp_coeff[1] = b[1]; hp_coeff[2] = b[2]; //here are the matlab "b" coefficients |
| 251 | + hp_coeff[3] = -a[1]; hp_coeff[4] = -a[2]; //the DSP needs the "a" terms to have opposite sign vs Matlab |
| 252 | +} |
| 253 | + |
| 254 | +void Compressor::setHPFilterCoeff(void) |
| 255 | +{ |
| 256 | + //https://www.keil.com/pack/doc/CMSIS/DSP/html/group__BiquadCascadeDF1.html#ga8e73b69a788e681a61bccc8959d823c5 |
| 257 | + //Use matlab to compute the coeff for HP at 20Hz: [b,a]=butter(2,20/(44100/2),'high'); %assumes fs_Hz = 44100 |
| 258 | + float32_t b[] = {9.979871156751189e-01, -1.995974231350238e+00, 9.979871156751189e-01}; //from Matlab |
| 259 | + float32_t a[] = { 1.000000000000000e+00, -1.995970179642828e+00, 9.959782830576472e-01}; //from Matlab |
| 260 | + setHPFilterCoeff_N2IIR_Matlab(b, a); |
| 261 | + //hp_coeff[0] = b[0]; hp_coeff[1] = b[1]; hp_coeff[2] = b[2]; //here are the matlab "b" coefficients |
| 262 | + //hp_coeff[3] = -a[1]; hp_coeff[4] = -a[2]; //the DSP needs the "a" terms to have opposite sign vs Matlab |
| 263 | +} |
| 264 | + |
| 265 | +void Compressor::updateThresholdAndCompRatioConstants(void) |
| 266 | +{ |
| 267 | + comp_ratio_const = 1.0f-(1.0f / comp_ratio); |
| 268 | + thresh_pow_FS_wCR = powf(thresh_pow_FS, comp_ratio_const); |
| 269 | +} |
| 270 | + |
| 271 | +void Compressor::setThreshPow(float32_t t_pow) |
| 272 | +{ |
| 273 | + thresh_pow_FS = t_pow; |
| 274 | + updateThresholdAndCompRatioConstants(); |
| 275 | +} |
| 276 | + |
| 277 | +// Accelerate the powf(10.0,x) function |
| 278 | +static float32_t pow10f(float32_t x) |
| 279 | +{ |
| 280 | + //return powf(10.0f,x) //standard, but slower |
| 281 | + return expf(2.302585092994f*x); //faster: exp(log(10.0f)*x) |
| 282 | +} |
| 283 | + |
| 284 | +// Accelerate the log10f(x) function? |
| 285 | +static float32_t log10f_approx(float32_t x) |
| 286 | +{ |
| 287 | + //return log10f(x); //standard, but slower |
| 288 | + return log2f_approx(x)*0.3010299956639812f; //faster: log2(x)/log2(10) |
| 289 | +} |
| 290 | + |
| 291 | +/* ---------------------------------------------------------------------- |
| 292 | +** Fast approximation to the log2() function. It uses a two step |
| 293 | +** process. First, it decomposes the floating-point number into |
| 294 | +** a fractional component F and an exponent E. The fraction component |
| 295 | +** is used in a polynomial approximation and then the exponent added |
| 296 | +** to the result. A 3rd order polynomial is used and the result |
| 297 | +** when computing db20() is accurate to 7.984884e-003 dB. |
| 298 | +** ------------------------------------------------------------------- */ |
| 299 | +//https://community.arm.com/tools/f/discussions/4292/cmsis-dsp-new-functionality-proposal/22621#22621 |
| 300 | +//float32_t log2f_approx_coeff[4] = {1.23149591368684f, -4.11852516267426f, 6.02197014179219f, -3.13396450166353f}; |
| 301 | +static float32_t log2f_approx(float32_t X) |
| 302 | +{ |
| 303 | + //float32_t *C = &log2f_approx_coeff[0]; |
| 304 | + float32_t Y; |
| 305 | + float32_t F; |
| 306 | + int E; |
| 307 | + |
| 308 | + // This is the approximation to log2() |
| 309 | + F = frexpf(fabsf(X), &E); |
| 310 | + // Y = C[0]*F*F*F + C[1]*F*F + C[2]*F + C[3] + E; |
| 311 | + //Y = *C++; |
| 312 | + Y = 1.23149591368684f; |
| 313 | + Y *= F; |
| 314 | + //Y += (*C++); |
| 315 | + Y += -4.11852516267426f; |
| 316 | + Y *= F; |
| 317 | + //Y += (*C++); |
| 318 | + Y += 6.02197014179219f; |
| 319 | + Y *= F; |
| 320 | + //Y += (*C++); |
| 321 | + Y += -3.13396450166353f; |
| 322 | + Y += E; |
| 323 | + |
| 324 | + return(Y); |
| 325 | +} |
0 commit comments