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

32位计算得到这些:

-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 ^个顺序部分.请注意,虽然您提供的代码可以直接用于实现过滤器部分,但您可能会意识到,作为额外的奖励,可以更有效地优化固定顺序构建块.

在导出滤波器部分的系数之前,我们退一步以获得滤波器设计参数.既然你提到滤波器是一个5阶低通巴特沃斯滤波器,我假设你选择省略a[0],b[0]哪两个都是1.从你提供的信息回来看,看起来滤波器是从一个生成的模拟滤波器的截止频率规格为45Hz,并使用双线性变换映射到以44100Hz的采样率工作的数字滤波器.作为完整性检查,可以从MATLAB或此applet确认过滤系数:

 a   b
  1, +1
  5, -4.9792522401964
 10, +9.91722403267282
 10, -9.87615728025693
  5, +4.91765142871949
  1, -0.97946594092826

获得等效滤波器所需的第一步是获得极点和零点.该过滤器的极点和零点可以通过以下任一方式获得:

• 用上面给出的a系数(对于零)和b系数(对于极点)分解多项式,或
• 通过对术语(ss _k)应用双线性变换,其中s _k是s平面中的极点,可以从例如维基百科(取代w _c = 2pi*(45/44100)弧度获得,给出您的滤波器设计规格).

得到的零在z = -1(所有5个).类似地,z平面中产生的极点是:

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

的2 ^次阶滤波器部分可以然后通过匹配磁极的复共轭对,并用零的等效数结合来获得.因此,将z = x + yj处的极点与其复共轭和2个零(z = -1)组合,滤波器部分的系数为:

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

类似地,1个^ST可以为真实的极点在Z = X的待获得的阶滤波器部,具有零一起(在z = -1):

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);

请注意,某些临时存储可以重复使用,但为了清楚起见,这是遗漏的.