我正在尝试在 C# 中使用 IIR LP 滤波器。它是一个 5 阶巴特沃斯滤波器。 该代码在 64 位模式下工作,但在 32 位模式下中断。调试显示,参数略有不同,输出增加到无穷大/NAN。 我正在使用 double 进行计算和存储。 正确的参数 a[i],b[i] 是:
-5 -4,9792522401964
10 9,91722403267282
-10 -9,87615728025693
5 4,91765142871949
-1 -0,979465940928259
32Bit计算得到这些:
-5 -4,97925281524658
10 9,91722583770752
-10 -9,87615966796875
5 4,91765308380127
-1 -0,979466378688812
过滤代码:
public void FilterBuffer(float[] srcBuf, long srcPos, float[] dstBuf, long dstPos, long nLen)
{
const double kDenormal = 0.000000000000001;
double denormal = m_invertDenormal ? -kDenormal : kDenormal;
m_invertDenormal = !m_invertDenormal;
for (int sampleIdx = 0; sampleIdx < nLen; sampleIdx++)
{
double sum = 0.0f;
m_inHistory[m_histIdx] = srcBuf[srcPos + sampleIdx] + denormal;
for (int idx = 0; idx < m_aCoeff.Length; idx++)
sum += m_aCoeff[idx] * m_inHistory[(m_histIdx - idx) & kHistMask];
for (int idx = 1; idx < m_bCoeff.Length; idx++)
sum -= m_bCoeff[idx] * m_outHistory[(m_histIdx - idx) & kHistMask];
m_outHistory[m_histIdx] = sum;
m_histIdx = (m_histIdx + 1) & kHistMask;
dstBuf[dstPos + sampleIdx] = (float)sum;
}
}
历史是 32 个条目,因此 histMask 为“31”以避免取模/比较...
有什么想法为什么这不起作用以及如何使其稳定吗?
最佳答案
IIR 滤波器因对数值精度敏感而臭名昭著,尤其是当递归方程中的项数增加时。幸运的是,在这种情况下,通过将滤波器实现为较小的 1st 和 2nd 阶部分的级联,可以获得不太敏感的不同滤波器拓扑结构.请注意,虽然您提供的代码可以直接用于实现过滤器部分,但您可能会意识到,作为额外的好处,可以更有效地优化那些固定顺序的构建 block 。
在导出滤波器部分的系数之前,我们退后一步以获得滤波器设计参数。由于您提到该滤波器是 5 阶低通巴特沃斯滤波器,我假设您选择省略 a[0]
和 b[0]
,它们是两者 1. 从您提供的信息回溯,看起来该滤波器是从截止频率规范为 45Hz 的模拟滤波器生成的,并使用双线性变换映射到以 44100Hz 采样率运行的数字滤波器。作为完整性检查,可以从 MATLAB 或 this applet 确认滤波器系数:
a b
1, +1
5, -4.9792522401964
10, +9.91722403267282
10, -9.87615728025693
5, +4.91765142871949
1, -0.97946594092826
获得等效滤波器所需的第一步是获得极点和零点。该滤波器的极点和零点可以通过以下任一方式获得:
- 使用上面给出的
a
系数(对于零点)和b
系数(对于极点)分解多项式,或者 - 通过对项 (s-sk) 应用双线性变换,其中 sk 是 s 平面中的极点,可以从例如 wikipedia (替换为 wc=2pi*(45/44100) 弧度,给定您的滤波器设计规范)。
生成的零点位于 z=-1(全部 5 个)。类似地,在 z 平面中产生的极点是:
0.998002182897612 +/- 0.00608551812433764 j
0.994819411183486 +/- 0.00374906249111718 j
0.993609052034203
然后可以通过匹配复数共轭极点对并与等效数量的零点组合来获得 2nd 阶滤波器部分。 因此,将 z = x + y j 处的极点与其复共轭和 2 个零点(z=-1 处)相结合,滤波器部分的系数为:
1, 1
2, -2x
1, x*x+y*y
类似地,对于 z = x 处的实极点以及零(z=-1 处),可以获得 1st 阶滤波器部分:
1, 1
1, -x
因此提供的第 5 阶滤波器的 3 个滤波器部分是:
// 1st section
// using poles at 0.998002182897612 +/- 0.00608551812433764j, and 2 zeros at -1
1, 1
2, -1.996004365795224
1, 0.996045390599240
// 2nd section
// using poles at 0.994819411183486 +/- 0.00374906249111718j, and 2 zeros at -1
1, 1
2, -1.989638822366971
1, 0.989679716337019
// 3rd section
// using pole at 0.993609052034203 and a zero at -1
1, 1
1, -0.993609052034203
如前所述,输出是通过级联这些部分生成的。从而获得诸如以下序列的东西:
f1.FilterBuffer(srcBuf, srcPos, tmpBuf1, 0, nLen);
f2.FilterBuffer(tmpBuf1, 0, tmpBuf2, 0, nLen);
f3.FilterBuffer(tmpBuf2, 0, dstBuf, dstPos, nLen);
请注意,一些临时存储可以重复使用,但为了清楚起见,将其省略。
关于c# - C# 中的 IIR 低通滤波器在 x86 模式下中断,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/23122901/