0
$\begingroup$

I have a signal that is NOT an audio signal (but looks like an audio signal - probably you can even "play" it). I would like to amplify the signal in regions where its amplitude is low. In very rare occasions the signal amplitude may be too high so I want to decrease it a bit. At the beginning and at the end of the 'track' I have some noise.

I think what what I want is called "Dynamic range compression". Is this correct? How I do this programmatically? I need some pseudo-code.

enter image description here

$\endgroup$
  • $\begingroup$ Did you check dsp.stackexchange.com/questions/10536/… $\endgroup$ – Laurent Duval Jul 21 '16 at 22:35
  • 2
    $\begingroup$ "Dynamc range compression" is a good label. so also is "Level compression" and "Automatic Gain Control" (AGC) are terms you might wanna look up. don't confuse "Level compression" with "Data compression". two different things. sometimes in the audio world it is ambiguous when someone refers to "compression of audio". to deal with the noise, you might need a "gate", which is also implemented with some audio level compressors. $\endgroup$ – robert bristow-johnson Jul 21 '16 at 22:35
  • $\begingroup$ @LaurentDuval-Yes. The question is a bit different. The answers: One of the answers is too generic the other one just provide some equations (and some possible drawback for using those equations as they are optimized for voice). $\endgroup$ – WeGoToMars Jul 21 '16 at 22:38
  • $\begingroup$ @robertbristow-johnson-THANKS. Good to know. I hope I can find more info now by searching for "Automatic Gain Control" $\endgroup$ – WeGoToMars Jul 21 '16 at 22:41
  • 1
    $\begingroup$ @MBaz, i don't think that mu-law or A-law (or arcsinh-law) is what this AGC is about. $\endgroup$ – robert bristow-johnson Jul 22 '16 at 0:53
3
$\begingroup$

from your posted waveform, i am assuming that this is a unipolar signal. that is

$$ x[n] \ge 0 $$

in audio, it would be the same, except that we would be working on $|x[n]|$ instead.

so first you want a sliding maximum of your signal, where the window length is $L$.

$$ x_1[n] = \max_{0 \le i < L} \Big| x[n-i] \Big| $$

since your input signal $x[n]$ appears to be unipolar, you can leave off the absolute value operation.

then you want to apply a low-pass filter (LPF) to that sliding max. you want the gain of the LPF at DC (0 Hz) to be 1 (or 0 dB gain). a simple first-order LPF is

$$ \begin{align} x_2[n] &= (1-p) \cdot x_1[n] + p \cdot x_2[n-1] \\ \\ &= x_1[n] + p \cdot (x_2[n-1] - x_1[n]) \\ \\ &= x_2[n-1] + (1-p) \cdot (x_1[n] - x_2[n-1]) \\ \end{align} $$

$p$ is the pole value of the LPF and

$$ 0 < 1-p \ll 1 $$

so

$x_2[n]$ will be an envelope for your input signal $x[n]$ and will be delayed by about half of the max window length plus about 4 times the "time constant" of the LPF:

$$d = \frac{L}{2} - \frac{4}{\log(p)} \quad \quad \text{samples}$$.

so, to normalize, you want to invert (compute the reciprocal of) the envelope $x_2[n]$ and multiply your input signal by that inverted envelope. the inverted envelope is

$$ x_3[n] = \frac{A}{x_2[n] + \epsilon} $$

$A$ is the normalized amplitude you want (it can be $A=1$ or $A=$ any other positive number that you like). $\epsilon$ is a tiny number, much smaller than most of your non-zero $x[n]$ that you need to add to the denominator to keep from dividing by zero (which is a bad thing).

but you should line up the signal and the delayed inverted envelope, so your normalized output is also delayed by $d$ samples (defined above):

$$ y[n] = x_3[n] \cdot x[n-d] $$

the most computationally expensive operation is the sliding max. we were just talking about this sliding max in the music-dsp mailing list. i will dig up C code for it and post that in a following answer.

$\endgroup$
  • $\begingroup$ Thanks Robert. I will try to implement this in Pascal. One question: As you can see, my signal has some areas where it is flat (no peaks). Will this affect the algorithm in any way? $\endgroup$ – WeGoToMars Aug 1 '16 at 9:55
  • $\begingroup$ about the flat parts, it will still look for the sliding max. but when it computes the inverse gain, it will be limited by the setting of that $\epsilon$ parameter. this prevents division by zero and prevents those flat areas from being amplified by a gain of $\infty$. $\endgroup$ – robert bristow-johnson Aug 3 '16 at 11:22
2
$\begingroup$

here's an efficient sliding maximum algorithm that has cost that is $O(\log_2(L))$. below window_width is $L$.

comes from

Brookes: "Algorithms for Max and Min Filters with Improved Worst-Case Performance" IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 47, NO. 9, SEPTEMBER 2000

#define A_REALLY_LARGE_NUMBER 3.40e38

typedef struct
   {
   unsigned long window_width;          // array_size/2 < window_width <= array_size
   unsigned long array_size;            // must be power of 2 for this simple implementation
   unsigned long input_index;           // the actual sample placement is at (array_size + input_index);
   float* big_array_base;               // the big array is malloc() separately and is actually twice array_size;
   } search_tree_array_data;


void initSearchArray(unsigned long window_width, search_tree_array_data* array_data)
   {
   array_data->window_width = window_width;

   array_data->array_size = 1;
   window_width--;
   while (window_width > 0)
      {
      array_data->array_size <<= 1;
      window_width >>= 1;
      }
   // array_size is a power of 2 such that
   // window_width <= array_size < 2*window_width
   // array_size = 2^ceil(log2(window_width)) = 2^(1+floor(log2(window_width-1)))

   array_data->input_index = 0;

   array_data->big_array_base = (float*)malloc(sizeof(float)*2*array_data->array_size);        // dunno what to do if malloc() fails.

   for (unsigned long n=0; n<2*array_data->array_size; n++)
      {
      array_data->big_array_base[n] = -A_REALLY_LARGE_NUMBER;        // init array.
      }                                                              // array_base[0] is never used.
    }



/*
 *   findMaxSample(value, &array_data) will place "value" into the circular
 *   buffer in the latter half of the array pointed to by array_data->big_array_base .
 *   it will then compare the value in "value" to its "sibling" value, takes the
 *   greater of the two and then pops up one generation to the parent node where 
 *   this parent also has a sibling and repeats the process.  since the other parent  
 *   nodes already have the max value of the two child nodes, when getting to the
 *   top-level parent node, this node will have the maximum value of all the samples
 *   in the big_array.  the number of iterations of this loop is ceil(log2(window_width)).
 */

float findMaxSample(float value, search_tree_array_data* array_data)
   {
   register float* big_array = array_data->big_array_base;

   register unsigned long index = array_data->array_size + array_data->input_index;        // our main buffer is in the latter half of the big array.

   while (index > 1UL)
      {
      big_array[index] = value;

      register float sibling_value = big_array[index ^ 1UL];        // toggle LSB, the upper bits of the sibling address are the same.

      if (value < sibling_value)
         {
         value = sibling_value;                        // use maximum of the two values
         }

      index >>= 1;                                     // parent address is index/2 (drop remainder or "sibling bit")
      }

   array_data->input_index++;
   if (array_data->input_index >= array_data->window_width)
      {
      array_data->input_index = 0;
      }

   return value;
   }
$\endgroup$
  • 1
    $\begingroup$ Interesting code. Unfortunately I cannot see what exactly is this pair of functions doing; (the findmax I mean) i.e., how do we use them in our main code ? And how does it compress a given buffer of samples ? (and in O($\log_2(L)$)) time?) $\endgroup$ – Fat32 Jul 22 '16 at 22:53
  • $\begingroup$ it's for envelope following, @Fat32. it's there to find the maximum value of the most current $L$ samples. so the output of this filter will often be a piecewise-constant function and will always "cover" the previous $L$. the compression is done in the other answer. i thought it best to put the $$ x_1[n] = \max_{0 \le i < L} \big| x[n] \big|$$ in another answer. $\endgroup$ – robert bristow-johnson Jul 22 '16 at 23:06
  • $\begingroup$ ok... I see...thanks! btw you mean $$ x_1[n] = \max_{0 \le i < L} \Big| x[n-i] \Big| $$ (excerpted from your other answer). To Sum Up: The code finds the max of the most recent $L$ samples of the audio buffer, and then compression is perfomed according to the other answer. What is the recommended min/max value of L? (assuming 44,100 Hz audio). I still have doubt of how can we find the max in logarithmic time? Also the sliding window is jumping by $L$ samples (into a complete new set of samples) or just moving one sample. I think I should post this to the other aswer $\endgroup$ – Fat32 Jul 22 '16 at 23:30
  • 1
    $\begingroup$ yes. thanks for the correction. i dunno, for an audio compressor/limiter, maybe $L=$ window_width can be about 4000 for 100 ms. that is a perceptual issue. it also depends on how fast you want your compressor/limiter to "pump". the algorithm for a compressor/limiter/gate is not exactly what is described above. the attack time and decay time on the envelope would be different for an audio compressor/limiter/gate. (and i didn't put in the gate alg. anyway.) $\endgroup$ – robert bristow-johnson Jul 23 '16 at 0:05
  • 1
    $\begingroup$ well, it's an AGC implemented without feedback and, in my opinion, equivalent to an audio level compressor with a sliding max on the input, identical attack and decay times (regarding the LPF), compression ratio set to $\infty$ to 1, and the "knee" set to $20 \log_{10}(A)$ in dB relative to whatever dB the level $1$ is. and delay compensating for the delay in the sliding max and the simple LPF. $\endgroup$ – robert bristow-johnson Jul 23 '16 at 0:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.