C# 中的 IIR 低通滤波器在 x86 模式下中断

Question

我正在尝试在 C# 中使用 IIR LP 滤波器。它是一个五阶巴特沃斯滤波器。 该代码在 64 位模式下工作，但在 32 位模式下中断。调试显示，参数略有不同，输出提高到无穷大/NAN。 我使用双打进行计算和存储。 正确的参数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”以避免模数/比较...

任何想法为什么这不起作用以及如何使其稳定？

Answer 1

IIR 滤波器因对数值精度敏感而臭名昭著，尤其是当递归方程中的项数增加时。幸运的是，在这种情况下，通过将滤波器实现为较小的 1^st 和 2^nd 阶部分的级联，可以获得不太敏感的不同滤波器拓扑.请注意，虽然您提供的代码可以直接用于实现过滤器部分，但您可能会意识到，作为额外的好处，可以更有效地优化那些固定顺序的构建块。

在推导滤波器部分的系数之前，我们后退一步来获取滤波器设计参数。由于您提到该滤波器是一个 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-s_k) 应用双线性变换，其中 s_k 是 s 平面中的极点，例如可以从 @987654322 获得@（根据您的滤波器设计规格，替换为 w_c=2pi*(45/44100) 弧度）。

生成的零位于 z=-1（全部 5 个）。类似地，在 z 平面中得到的极点是：

0.998002182897612 +/- 0.00608551812433764 j
0.994819411183486 +/- 0.00374906249111718 j
0.993609052034203

2^nd 阶滤波器部分然后可以通过匹配复共轭极点对，并与等效数量的零点组合来获得。 因此，将 z = x + y j 处的极点与其复共轭和 2 个零点（在 z=-1 处）相结合，滤波器部分的系数为：

1, 1
2, -2x
1, x*x+y*y

同样，对于 z = x 处的实极点以及零点（在 z=-1 处），可以获得 1^st 阶滤波器部分：

1, 1
1, -x

提供的 5^th 阶过滤器的 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);

请注意，一些临时存储可以重复使用，但为了清楚起见，将其省略。