Dynamic range compression - Need pseudo algorithm for normalizing a signal

Question

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.

Answer 1

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.

Answer 2

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;
   }