how to create a simple iir low pass filter with not round errors? (16 bit pcm data)

Question

i have an array of n length fullfilled by 16 bit (int16) pcm raw data,the data is in 44100 sample_rate and stereo,so i have in my array first 2 bytes left channel then right channel etc...i tried to implement a simple low pass converting my array into floating points -1 1,the low pass works but there are round errors that cause little pops in the sound now i do simply this :

    INT32  left_id  = 0;
    INT32  right_id = 1;
    DOUBLE  filtered_l_db = 0.0;
    DOUBLE  filtered_r_db = 0.0;
    DOUBLE  last_filtered_left  = 0;
    DOUBLE  last_filtered_right = 0;  
    DOUBLE  l_db = 0.0;
    DOUBLE  r_db = 0.0; 
    DOUBLE  low_filter  =  filter_freq(core->audio->low_pass_cut);  
    for(UINT32 a = 0; a < (buffer_size/2);++a)  
    {    

        
  
        l_db = ((DOUBLE)input_buffer[left_id])  / (DOUBLE)32768;
        r_db = ((DOUBLE)input_buffer[right_id]) / (DOUBLE)32768;      
        ///////////////LOW PASS
        filtered_l_db = last_filtered_left  + 
        (low_filter * (l_db -last_filtered_left ));
        filtered_r_db = last_filtered_right + 
        (low_filter * (r_db - last_filtered_right));   
        last_filtered_left  = filtered_l_db;
        last_filtered_right = filtered_r_db;
  
        INT16 l = (INT16)(filtered_l_db * (DOUBLE)32768);
        INT16 r = (INT16)(filtered_r_db * (DOUBLE)32768);
        output_buffer[left_id]  =  (output_buffer[left_id]  + l);
        output_buffer[right_id] =  (output_buffer[right_id] + r);     
          
      left_id +=2;
      right_id +=2;
    }

PS: the input buffer is an int16 array with the pcm data from -32767 to 32767;

i found this function here Low Pass filter in C

and was the only one that i could understand xd

   DOUBLE filter_freq(DOUBLE cut_freq)
   {
     DOUBLE a = 1.0/(cut_freq * 2 * PI);
     DOUBLE b = 1.0/SAMPLE_RATE;
     return b/(a+b);  
   }

my aim is instead to have absolute precision on the wave,and to directly low pass using only integers with the cost to lose resolution on the filter(and i'm ok with it)..i saw a lot of examples but i really didnt understand anything...someone of you would be so gentle to explain how this is done like you would explain to a little baby?(in code or pseudo code rapresentation) thank you

Can you please include all variable declarations. Are they all `float` or just some of them? — Lundin, Feb 19 '21 at 14:32
Please show/link in your question the examples you saw. There are still several variable definitions missing. You can do the same calculation using integer calculation. — Bodo, Feb 19 '21 at 14:48
How your buffers are organised is very unclear to me. Here `input_buffer[b][left_id]`, you are incrementing both `b`and `left_id` — Damien, Feb 19 '21 at 14:48
65536 /2 32767 if unsigned(0 65535)...btw i don't wanna use floating points at all so it's not relevant ...to damien : yes sorry cause i have multiple samples,now i edit it to be more clear i forgot and how whould you do it with integers? — 1bytecharhaterxxx, Feb 19 '21 at 14:48

score 1 · Accepted Answer · answered Feb 19 '21 at 15:43

Assuming the result of function filter_freq can be written as a fraction m/n your filter calculation basically is

y_new = y_old + (m/n) * (x - y_old);

which can be transformed to

y_new = ((n * y_old) + m * (x - y_old)) / n;

The integer division / n truncates the result towards 0. If you want rounding instead of truncation you can implement it as

y_tmp = ((n * y_old) + m * (x - y_old));
if(y_tmp < 0) y_tmp -= (n / 2);
else y_tmp += (n / 2);
y_new = y_tmp / n

In order to avoid losing precision from dividing the result by n in one step and multiplying it by n in the next step you can save the value y_tmp before the division and use it in the next cycle.

y_tmp = (y_tmp + m * (x - y_old));
if(y_tmp < 0) y_new = y_tmp - (n / 2);
else y_new = y_tmp + (n / 2);
y_new /= n;

If your input data is int16_t I suggest to implement the calculation using int32_t to avoid overflows.

I tried to convert the filter in your code without checking other parts for possible problems.


    INT32  left_id  = 0;
    INT32  right_id = 1;
    int32_t filtered_l_out = 0; // output value after division
    int32_t filtered_r_out = 0;
    int32_t  filtered_l_tmp = 0; // used to keep the output value before division
    int32_t  filtered_r_tmp = 0;  
    int32_t  l_in = 0; // input value
    int32_t  r_in = 0; 

    DOUBLE  low_filter  =  filter_freq(core->audio->low_pass_cut);
    // define denominator and calculate numerator
    // use power of 2 to allow bit-shift instead of division
    const uint32_t filter_shift = 16U;
    const int32_t filter_n = 1U << filter_shift;
    int32_t filter_m = (int32_t)(low_filter * filter_n)

    for(UINT32 a = 0; a < (buffer_size/2);++a)  
    {    
  
        l_in = input_buffer[left_id]);
        r_in = input_buffer[right_id];

        ///////////////LOW PASS
        filtered_l_tmp = filtered_l_tmp + filter_m * (l_in - filtered_l_out);
        if(last_filtered_left < 0) {
            filtered_l_out = last_filtered_left - filter_n/2;
        } else {
            filtered_l_out = last_filtered_left + filter_n/2;
        }
        //filtered_l_out /= filter_n;
        filtered_l_out >>= filter_shift;

        /* same calculation for right */
  
        INT16 l = (INT16)(filtered_l_out);
        INT16 r = (INT16)(filtered_r_out);

        output_buffer[left_id]  =  (output_buffer[left_id]  + l);
        output_buffer[right_id] =  (output_buffer[right_id] + r);     

        left_id +=2;
        right_id +=2;
    }

As your filter is initialized with 0 it may need several samples to follow a possible step to the first input value. Depending on your data it might be better to initialize the filter based on the first input value.

Hi thank you ,you were so gentle and your code did the job,i thought that the problem was in the conversion,but simply is that the new amplitude near the cuttoff produce the pops,well sorry for wasting your time,i feel really embarassed, now i will try to implement this: https://www.dsprelated.com/freebooks/filters/Simplest_Lowpass_Filter_I.html#eq:simplp in figure 1.1 — 1bytecharhaterxxx, Feb 19 '21 at 17:01

how to create a simple iir low pass filter with not round errors? (16 bit pcm data)

1 Answers1