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
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);
}
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');
%% 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);
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)?
b0,b1,b2,a1, anda2beingshortand you casting them tointin the inner computation (which is repeated), why not declare those coef variables asintand cast them when loading them from the coef array? $\endgroup$b0,b1,b2coefficients 20 times bigger before quantizing them to the nearest integer and see if you get better agreement between the MATLAB and C filtering? $\endgroup$