如何创建一个没有舍入误差的简单 iir 低通滤波器？ （16 位 pcm 数据）

Posted 2023-02-25

【中文标题】如何创建一个没有舍入误差的简单 iir 低通滤波器？ （16 位 pcm 数据）

【英文标题】：how to create a simple iir low pass filter with not round errors? (16 bit pcm data)

【发布时间】：2021-02-19 14:29:09

【问题描述】：

我有一个长度为 n 的数组，由 16 位 (int16) pcm 原始数据填充，数据为 44100 sample_rate 和立体声，所以我的数组中有前 2 个字节的左声道，然后是右声道等...我尝试实现一个简单的低通，将我的数组转换为浮点数 -1 1，低通有效，但存在导致圆形错误声音中的小爆裂声 现在我只做这个：

    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：输入缓冲区是一个int16数组，pcm数据从-32767到32767；

我在这里找到了这个功能 Low Pass filter in C

是我唯一能理解的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);  
   

我的目标是在波形上具有绝对精度，并仅使用整数直接低通 以失去过滤器分辨率为代价（我没问题）..我看到了很多例子，但我真的什么都不懂......你们中的某个人会像你一样温柔地解释这是如何完成的给小宝宝解释一下？（用代码或伪代码表示）谢谢

【问题讨论】：

能否请您包括所有变量声明。他们都是float 还是只是其中的一部分？

好的，抱歉：D

嗯，我真的帮不了你，但为什么是 32768 而不是 32767？

请在您的问题中显示/链接您看到的示例。仍然缺少几个变量定义。您可以使用整数计算进行相同的计算。

缓冲区的组织方式对我来说非常不清楚。这里input_buffer[b][left_id]，你正在增加b和left_id

【参考方案1】：

假设函数filter_freq的结果可以写成分数m/n你的过滤器计算基本上是

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

可以转化为

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

整数除法/ n 将结果截断为 0。如果您想要舍入而不是截断，您可以将其实现为

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

为了避免在一步中将结果除以n 并在下一步将其乘以n 时丢失精度，您可以在除法之前保存值y_tmp 并在下一个循环中使用它。

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;

如果您的输入数据是int16_t，我建议使用int32_t 来实现计算以避免溢出。

我尝试在您的代码中转换过滤器，而不检查其他部分是否存在问题。

    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;
    

由于您的过滤器初始化为 0，因此可能需要多个样本才能按照可能的步骤到达第一个输入值。根据您的数据，根据第一个输入值初始化过滤器可能会更好。

【讨论】：

您好，谢谢您，您非常温柔，您的代码完成了这项工作，我认为问题出在转换中，但只是截止点附近的新幅度产生了爆裂声，抱歉浪费你的时间，我真的很尴尬，现在我将尝试实现这一点：dsprelated.com/freebooks/filters/… 图 1.1