我正在实现一个具有 Q1.15 精度的定点直接形式 I biquad 滤波器。我使用分数保存来改善量化误差。但是，我的函数给了我与 MATLAB 比较的不同结果filter。我没发现有什么问题，请帮忙。我的代码如下

void iir_biquad_df1_16(
    const short* pCoeffs,             /* b0, b1, b2, a1, a2 */
    short* pStates,                   /* x[n-1], x[n-2], y[n-1], y[n-2], fraction_state */
    const short* pSrc,                /* Input buffer */
    short* pDst,                      /* Output buffer */
    unsigned int blockSize,           /* Frame length */
    const unsigned int postShift,     /* Post shift when coefficients out of range [-1,1) */
    const unsigned int nSections      /* Number of cascaded sections */
)
{
    short* pIn = pSrc;                              /* Source pointer */
    short* pOut = pDst;                             /* Destination pointer */
    short* pState = pStates;                        /* State pointer */
    const short* pCoeff = pCoeffs;                  /* Coefficient pointer */
    int acc;                                        /* Accumulator */
    int b0, b1, b2, a1, a2;                         /* Filter coefficients */
    int x1, x2, y1, y2;                             /* Filter state */
    int xIn;                                        /* Temporary input */
    unsigned int shift = 15U - postShift;           /* Post shift for output */
    unsigned int mask = (1 << shift) - 1;           /* Fraction saving */
    int saturation = (0x00008000 << shift) - 1;
    unsigned int samples, section = nSections;      /* Loop counters */
    
    do
    {
        /* Read the coefficients */
        b0 = (int)(*pCoeff++);
        b1 = (int)(*pCoeff++);
        b2 = (int)(*pCoeff++);
        a1 = (int)(*pCoeff++);
        a2 = (int)(*pCoeff++);
        
        /* Read the state values */
        x1 = (int)(*pState++);
        x2 = (int)(*pState++);
        y1 = (int)(*pState++);
        y2 = (int)(*pState++);
        acc = (int)(*pState);
        pState -= 4;
        
        samples = blockSize;
        
        while (samples > 0U)
        {
            xIn = (int)(*pIn++);
            /* acc =  b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2] */
            acc += b0*xIn + b1*x1 + b2*x2 - a1*y1 - a2*y2;
            
            if (acc >  saturation) acc =  saturation;         /* saturate if necessary */
            if (acc < -saturation) acc = -saturation;
            
            x2 = x1;
            x1 = xIn;
            y2 = y1;
            y1 = acc >> shift;
            
            
            *pOut++ = (short)y1;
            acc &= mask;

            samples--;
        }

        /* Store the updated state variables back into the pState array */
        *pState++ = (short)x1;
        *pState++ = (short)x2;
        *pState++ = (short)y1;
        *pState++ = (short)y2;
        *pState++ = (short)acc;

        /* Subsequent sections take previous output buffer as input */
        pIn = pDst;
        pOut = pDst;
    
    } while (--section);
}

我有 100 个输入信号样本和两个级联双二阶 IIR 滤波器。块大小是 50，所以我有 2 个块要按顺序处理。测试的主要功能由下式给出

int main()
{
    short x[100] = { 10313,13297,-12223,13546,4337,-13188,-7258,1536,14992,15233,-11219,15420,14980,-479,9840,-11735,-2564,13623,9575,15057,5103,-15214,11440,14221,5857,8446,7967,-3532,5095,-10775,6752,-15341,-7310,-14871,-13201,10599,6384,-5993,14753,-15255,-2007,-3881,8700,9673,-10261,-335,-1783,4794,6860,8346,-7339,5888,5082,-11056,-12485,-54,15065,-5230,2794,-9050,8234,-8025,195,6523,12809,15050,1547,-11842,-11492,-7946,11165,-8052,10298,-8404,14066,-4916,-9942,-8156,3803,-875,-4861,10841,2794,1629,13671,-7018,8428,8314,-3918,2222,-13898,-14616,1009,9148,14222,-12127,2255,-1003,-15994,-5337 };
    short y[100] = { 0 };
    short coeffs[5 * 2] = { 575, -1150, 575, -31232, 14931, 575, 1150, 575, -32621, 16238 };
    short states[5 * 2] = { 0 };
    unsigned int blockSize = 50;
    unsigned int postShift = 1;
    unsigned int nSections = 2;
    // the first block of 50 samples
    iir_biquad_df1_16(coeffs, states, x, y, blockSize, postShift, nSections);
    // the second block of 50 samples
    iir_biquad_df1_16(coeffs, states, x + blockSize, y + blockSize, blockSize, postShift, nSections);

    return 0;
}

输入信号为 MATLAB 生成的随机噪声，滤波器系数也在 MATLAB 中设计。MATLAB 验证代码是

% clear
%% filter design
fs = 48000;
f1 = 50;
f2 = 600;
w1 = f1 / (fs / 2);
w2 = f2 / (fs / 2);
N = 2;
[b_bp, a_bp] = butter(N, [w1, w2], 'bandpass');

%% convert to cascade biquad
[sos,g] = tf2sos(b_bp, a_bp);
b1 = sos(1,1:3) * sqrt(g); a1 = sos(1, 4:6);
b2 = sos(2,1:3) * sqrt(g); a2 = sos(2, 4:6);

postShift = 1;
b1_q14 = int16(b1*2^14);
b2_q14 = int16(b2*2^14);
a1_q14 = int16(a1*2^14);
a2_q14 = int16(a2*2^14);

%% input signal generation
rng(0)
x = rand(1, 100) - 0.5; % uniform distribution in range (-0.5, 0.5)
x_q15 = int16(x*2^15);
x = x';

%% output signal
y1 = filter(b1, a1, x);
y2 = filter(b2, a2, y1);
y1_q15 = int16(y1*2^15);
y2_q15 = int16(y2*2^15);

由于某些滤波器系数超出范围$[-1, 1)$，我施加了一个后移，它缩小了范围内的所有滤波器系数$[-1, 1)$. 因此，到定点的转换乘以$2^{14}$并给出 Q2.14 精度。

经过调试，发现第一个biquad段的输出信号似乎是正确的，但是级联第二段后，最终的输出和MATLAB给出的不一样。这个问题困扰了我两天。任何建议将不胜感激。

编辑：我在双二阶滤波器的第二部分之前和之后绘制信号。输入信号是随机噪声

第一段滤波后的滤波信号为

我们可以看到 C 输出与 MATLAB 输出几乎相同。但是，在第二个双二阶之后，它们就大不相同了。

这是否意味着 16 位字长对于我的应用程序（音频处理）来说不够长？