effect_dynamics.cpp

/* Audio Library for Teensy
   Dynamics Processor (Gate, Compressor & Limiter)
   Copyright (c) 2017, Marc Paquette (marc@dacsystemes.com)
   Based on analyse_rms & mixer objects by Paul Stoffregen

   Development of this audio library was funded by PJRC.COM, LLC by sales of
   Teensy and Audio Adaptor boards.  Please support PJRC's efforts to develop
   open source software by purchasing Teensy or other PJRC products.

   Permission is hereby granted, free of charge, to any person obtaining a copy
   of this software and associated documentation files (the "Software"), to deal
   in the Software without restriction, including without limitation the rights
   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   copies of the Software, and to permit persons to whom the Software is
   furnished to do so, subject to the following conditions:

   The above copyright notice, development funding notice, and this permission
   notice shall be included in all copies or substantial portions of the Software.

   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
   THE SOFTWARE.
*/

#include "effect_dynamics.h"
#include "utility/dspinst.h"
#include "utility/sqrt_integer.h"

static float analyse_rms(int16_t *data) {
  //Serial.println(__FUNCTION__);
  uint32_t *p = (uint32_t *)data;
  const uint32_t *end = p + AUDIO_BLOCK_SAMPLES / 2;
  int64_t sum = 0;
  do {
    uint32_t n1 = *p++;
    uint32_t n2 = *p++;
    uint32_t n3 = *p++;
    uint32_t n4 = *p++;
    sum = multiply_accumulate_16tx16t_add_16bx16b(sum, n1, n1);
    sum = multiply_accumulate_16tx16t_add_16bx16b(sum, n2, n2);
    sum = multiply_accumulate_16tx16t_add_16bx16b(sum, n3, n3);
    sum = multiply_accumulate_16tx16t_add_16bx16b(sum, n4, n4);

  } while (p < end);
  if (sum == 0) return 0;
  int32_t meansq = sum / AUDIO_BLOCK_SAMPLES;
  return sqrt_uint32(meansq) / 32767.0f;
}

static void applyGain(int16_t *data, int32_t mult1, int32_t mult2) {
  //Serial.println(__FUNCTION__);
  uint32_t *p = (uint32_t *)data;
  const uint32_t *end = p + AUDIO_BLOCK_SAMPLES / 2;
  int32_t inc = (mult2 - mult1) / (AUDIO_BLOCK_SAMPLES / 2);

  do {
    uint32_t tmp32 = *p; // read 2 samples from *data
    int32_t val1 = signed_multiply_32x16b(mult1, tmp32);
    mult1 += inc;
    int32_t val2 = signed_multiply_32x16t(mult1, tmp32);
    mult1 += inc;
    val1 = signed_saturate_rshift(val1, 16, 0);
    val2 = signed_saturate_rshift(val2, 16, 0);
    *p++ = pack_16b_16b(val2, val1);
  } while (p < end);
}

/* ----------------------------------------------------------------------
  https://community.arm.com/tools/f/discussions/4292/cmsis-dsp-new-functionality-proposal/22621#22621
  Fast approximation to the log2() function.  It uses a two step
  process.  First, it decomposes the floating-point number into
  a fractional component F and an exponent E.  The fraction component
  is used in a polynomial approximation and then the exponent added
  to the result.  A 3rd order polynomial is used and the result
  when computing db20() is accurate to 7.984884e-003 dB.
** ------------------------------------------------------------------- */

float log2f_approx_coeff[4] = {1.23149591368684f, -4.11852516267426f, 6.02197014179219f, -3.13396450166353f};

float log2f_approx(float X)
{
  float *C = &log2f_approx_coeff[0];
  float Y;
  float F;
  int E;

  // This is the approximation to log2()
  F = frexpf(fabsf(X), &E);

  //  Y = C[0]*F*F*F + C[1]*F*F + C[2]*F + C[3] + E;
  Y = *C++;
  Y *= F;
  Y += (*C++);
  Y *= F;
  Y += (*C++);
  Y *= F;
  Y += (*C++);
  Y += E;
  return (Y);
}

// https://codingforspeed.com/using-faster-exponential-approximation/
inline float expf_approx(float x) {
  x = 1.0f + x / 1024;
  x *= x; x *= x; x *= x; x *= x;
  x *= x; x *= x; x *= x; x *= x;
  x *= x; x *= x;
  return x;
}

inline float unitToDb(float unit) {
  return 6.02f * log2f_approx(unit);
}

inline float dbToUnit(float db) {
  return expf_approx(db * 2.302585092994046f * 0.05f);
}

void AudioEffectDynamics::update(void) {

  audio_block_t *block;

  block = receiveWritable(0);

  if (!block) return;

  if (!gateEnabled && !compEnabled && !limiterEnabled) {

    //Transmit & release
    transmit(block);
    release(block);
    return;
  }

  for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) {

    unsigned int sampleIndexPlus1 = (sampleIndex + 1) % sampleBufferSize;

    uint32_t sampleToRemove = samplesSquared[sampleIndexPlus1];
    sumOfSamplesSquared -= (sampleToRemove * sampleToRemove);

    int16_t sample = block->data[i];
    samplesSquared[sampleIndex] = abs(sample);
    uint32_t sampleSquared = sample * sample;
    sumOfSamplesSquared += sampleSquared;

    sampleIndex = (sampleIndex + 1) % sampleBufferSize;

    float rms = sqrt(sumOfSamplesSquared / float(sampleBufferSize)) / 32768.0;

    //Compute block RMS level in Db
    float inputdb = MIN_DB;
    if (rms > 0) inputdb = unitToDb(rms);

    //Gate
    if (gateEnabled) {
      if (inputdb >= gateThresholdOpen) gatedb = (aGateAttack * gatedb) + (aOneMinusGateAttack * MAX_DB);
      else if (inputdb < gateThresholdClose) gatedb = (aGateRelease * gatedb) + (aOneMinusGateRelease * MIN_DB);
    } else gatedb = MAX_DB;

    //Compressor
    if (compEnabled) {
      float attdb = MAX_DB; //Below knee
      if (inputdb >= aLowKnee) {
        if (inputdb <= aHighKnee) {
          //Knee transition
          float knee = inputdb - aLowKnee;
          attdb = aKneeRatio * knee * knee * aTwoKneeWidth;
        } else {
          //Above knee
          attdb = compThreshold + ((inputdb - compThreshold) * compRatio) - inputdb;
        }
      }
      if (attdb <= compdb) compdb = (aCompAttack * compdb) + (aOneMinusCompAttack * attdb);
      else compdb = (aCompRelease * compdb) + (aOneMinusCompRelease * attdb);
    } else compdb = MAX_DB;

    //Brickwall Limiter
    if (limiterEnabled) {
      float outdb = inputdb + compdb + makeupdb;
      if (outdb >= limitThreshold) limitdb = (aLimitAttack * limitdb) +
                                               (aOneMinusLimitAttack * (limitThreshold - outdb));
      else limitdb *= aLimitRelease;
    } else limitdb = MAX_DB;

    //Compute linear gain
    float totalGain = gatedb + compdb + makeupdb + limitdb;

    float multiplier = dbToUnit(totalGain);
    int16_t result = sample * multiplier;
    block->data[i] = result;
    //Apply gain to block
  }

  //Transmit & release
  transmit(block);
  release(block);
}