# Wrong result when realizing fixed-point biquad IIR filter

I'm implementing a fixed-point direct form I biquad filter with Q1.15 precision. I use fraction saving to improve quantization error. However, my function gives me the different results compared by MATLAB's filter. I didn't find anything wrong, please help. My code is as follows

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

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

I have 100 samples for input signal, and two cascaded biquad IIR filters. The block size is 50 so I have 2 blocks to process in order. The main function for test is given by

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

The input signal is random noise generated in MATLAB and the filter coefficients is also designed in MATLAB. The MATLAB validation code is

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

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

Since some of the filter coefficients are out of range of $$[-1, 1)$$, I impose a post shift that shrinks all filter coefficients in the range of $$[-1, 1)$$. Thus the conversion to fixed-point is multiplied by $$2^{14}$$ and gives a Q2.14 precision.

After debugging, I found that the output signal of the first biquad section seems to be correct, but after cascading the second section, the final output is not the same as MATLAB gives. This problem has bothered me for two days. Any advice would be greatly appreciated.

Edit: I plot the signal before and after the second section of biquad filter. The input signal is random noise

and the filtered signal after the first secion of filter is

We can see that the C output is nearly identical with the MATLAB output. However, after the second biquad, they are quite different.

Does it mean that 16-bit word length is not long enough for my application (audio processing)?

• If the first block works fine, then the code probably works. Is your output on the same format as your input? It should be...
– Ben
Nov 3, 2021 at 13:32
• @PeterK. No, I got the same result. The filter states are just two previous samples of input and output, I didn't find anything wrong with the states. Nov 4, 2021 at 1:35
• Instead of b0, b1, b2, a1, and a2 being short and you casting them to int in the inner computation (which is repeated), why not declare those coef variables as int and cast them when loading them from the coef array? Nov 5, 2021 at 4:11
• @robertbristow-johnson and btw, fraction saving is learned from one of your answers : ) I've learned a lot on this site! Nov 5, 2021 at 6:12
• I think it's coefficient quantization. Your input signal is $\pm$0.5 and your output signal is about $\frac{1}{30}$ of the size. Why not scale the b0, b1, b2 coefficients 20 times bigger before quantizing them to the nearest integer and see if you get better agreement between the MATLAB and C filtering? Nov 6, 2021 at 6:21