如何创建一个没有舍入误差的简单 iir 低通滤波器? (16 位 pcm 数据)
how to create a simple iir low pass filter with not round errors? (16 bit pcm data)
我有一个 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);
}
我的目标是在波形上有绝对的精度,并且只使用整数直接低通
以失去过滤器的分辨率为代价(我没意见)..我看到了很多例子,但我真的什么都不懂...你们中的某个人会像你一样温柔地解释这是如何完成的能给小宝宝解释一下吗?(用代码或伪代码表示)谢谢
假设函数 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 初始化的,因此可能需要几个样本才能按照可能的步骤到达第一个输入值。根据您的数据,最好根据第一个输入值初始化过滤器。
我有一个 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);
}
我的目标是在波形上有绝对的精度,并且只使用整数直接低通 以失去过滤器的分辨率为代价(我没意见)..我看到了很多例子,但我真的什么都不懂...你们中的某个人会像你一样温柔地解释这是如何完成的能给小宝宝解释一下吗?(用代码或伪代码表示)谢谢
假设函数 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 初始化的,因此可能需要几个样本才能按照可能的步骤到达第一个输入值。根据您的数据,最好根据第一个输入值初始化过滤器。