# Strange Results from Biquad Peaking EQ Filter

I have a signal that I am filtering on an MCU (STM32F405) that I get from a codec (CS42L52). It is being sampled at a sampling frequency of 48kHz (actually between 48400 and 48600). I have a ping pong buffer setup with DMA and am using I2S to send the signal between the MCU and codec.

I seem to have a glitch somewhere when filtering. It is somehow related to the buffer size, it happens at (buffer_size / 4) samples, but only when I am filtering. Below is a picture of both channels of my output, the top signal is the left channel unfiltered and the bottom is the right channel filtered.

I can switch the channel that is filtered and have the same results on the opposite channel. It is a fixed point filter with coefficients in Q29 format and look like

531850828 (b0), -952252845 (b1), 507513243 (b2), 952252845 (a1), -502493159 (a2)


for a center frequency of 3150Hz, Q of 7.2, Gain of -3, and sampling frequency of 48kHz. My filter uses the equation

y = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]


I am able to process all of the data in my buffer in less than 5ms depending on buffer size.

While the glitch appears at all frequencies, it tends to be most prominent around the frequency that I am filtering. Does this mean my coefficients are somehow incorrect?

Any ideas on what I am doing wrong?

## Edit:

I am using the CMSIS library from ARM to do my filtering. I have followed this example from the CMSIS documentation. According to the documentation, I need to scale my coefficients to be between [-1, +1) which is why I use Q29 format, and then shift them by 2. There is a postshift parameter coded into the Q31 filter. Here is the CMSIS code for the filter.

From the codec I get 24-bits of data which I pass through the filter and then send back to the codec. The accumulator is 64 bit and the code is in C. I guess I have made an assumption that 24-bits is small enough to not have overflow from the filter, but I am new to signal processing. How do I correct for overflow?

## Edit 2

I have been playing with the coefficients and when I set the gain to 0, the signal comes out clean. If I leave all the parameters except gain the same and change the gain (I tried +1.0 and -1.0), then I have distorted output. This is making me think more and more that it is arithmetic overflow. I tried multiplying the input of the filter by 0.125 to scale it down by 1/8, but that did not seem to help. What is the proper way to deal with this?

## Edit 3

Using the code RBJ provided, I get the same sort of output as before. This time I have converted the input into floating point format and then back to q31 for output. If I change the coefficients to have a 0 gain, I get a clean output, but as soon as I change the gain to anything but 0, it creates noise.

In my code, I am filtering one value at a time through filterBuffer. When I receive the data from DMA, I have to rotate the data right by 16 places to correct for the DMA's FIFO packing the data in the wrong order. For example, when the codec sends 24-bits 0xABCD EF00, the DMA will pack it as 0xEF00 ABCD. This is fixed by using the asm ror #16 instruction.

// f0 = 3150, Q = 4.3, Gain = 3.0, Fs = 48000
sectionData *coeffs = malloc(sizeof(sectionData));
coeffs->filterCoefficients[0] = 1.0155644882603565; // b0
coeffs->filterCoefficients[1] = -1.7632428285856931; // b1
coeffs->filterCoefficients[2] = 0.90897819415372805; // b2
coeffs->filterCoefficients[3] = 1.7632428285856931; // -a1
coeffs->filterCoefficients[4] = -0.92454268241408455; // -a2

static float32_t states[4] = {0.0,0.0,0.0,0.0}

/* loop in main after receiving data in rxbuf */
for (i = 0; i < BUFFERSIZE; i++) {
*ptrIn = *rxbuf++;
asm("mov %[result], %[value], ror #16" : [result] "=r" (*ptrIn) : [value] "r" (*ptrIn));

fIn = (float32_t) *ptrIn / 2147483648u;
filterBuffer(1, 1, &fIn, &fOut, &states[0], &coeffs);
filtered[i] = fOut * 2147483648u;

asm("mov %[result], %[value], ror #16" : [result] "=r" (*filtered[i]) : [value] "r" (filtered[i]));
asm("mov %[result], %[value], ror #16" : [result] "=r" (*ptrIn) : [value] "r" (*ptrIn));

// separate left and right channel, only filter left channel.
if (i%2 == 0)
*txbuf++ = (uint32_t)filtered[i];
else
*txbuf++ = *ptrIn++;
}

• how, exactly, are you doing your fixed-point arithmetic? by "Q29", do you mean that there are 3 bits left of the binary point? and how are you calculating coefficients. it appears to me that you got the sign of $a_2$ correct (if you're adding in the difference equation). but since your coefficients are fixed point, somewhere in your code, you have to be shifting right by 29 bits. is your accumulator 64 bits? are you doing this in C? are you using long and long long types or explicit int32 and int64? you probably have an arithmetic overflow and you don't know it yet. – robert bristow-johnson Oct 16 '14 at 19:42
• consider posting a snippet of code. – robert bristow-johnson Oct 16 '14 at 19:43
• also, if you're coding this equation explicitly > y = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] then you will have circular address wrap-around issues similarly to when folks encode an FIR filter or long delay line. "While the glitch appears at all frequencies, it tends to be most prominent around the frequency that I am filtering." -- i'll bet it's arithmetic overflow in your fixed-point registers. you might need to go to Q4.28 or Q5.27 or Q6.26 . you have plenty of bits to the right, may as well pad a couple more on the left. – robert bristow-johnson Oct 16 '14 at 19:46
• Thanks for the reply! I am using your equations to calculate the coefficients from here. And then changing the sign of the 'a' coefficients to fit my filter equation. In order to put them in to Q29 format, I am taking the floating point value I get from your equations and then multiplying them by 2^29, then truncating them as integers. I edited the post to answer your other questions. – Rick Oct 16 '14 at 20:45
• It looks like the state of your filter is reset at the beginning of every block. Are you sure it is persistent? Check the scope of all variables which might store the state of the filter. – pichenettes Oct 16 '14 at 21:04

Rick, just because i'm lazy, i don't wanna rewrite your code, but i'm pasting a snippet of old code of mine that does cascaded biquad filtering. it's in single-precision float, but could be modified to be Q3.29 (or i would recommend Q4.28) with use of the q64_t and q32_t types and a little right shifting. it swaps the nesting of the sections and the block of samples, it does the sections (or "stages") on the inside loop and reuses those states as i have described. but your nesting (with the section loop on the outside) might be better (faster).

just an idea about how to structure the code. otherwise, it's gonna be about how to do the fixed-point arithmetic, and i surely don't see why you're doing this bit-field manipulation. just shift the 64-bit acc to the right by 28 or 29 bits (depending on your scaling factor).

typedef struct { float filterCoefficients[5]; /* b0, b1, b2, a1, a2 (a0 is normalized to 1) */ } sectionData;

int filterBuffer(long Nsamples, int Nsections, float *input, float *output, float *filterStates, sectionData **sections) { register long n; register int i; register sectionData **section_ptr; register float sampleValue, state1, state2, *state_ptr, *coef_ptr;

for (n=0; n<Nsamples; n++)
{
section_ptr = sections;                         /* reset section pointer */
state_ptr = filterStates;                       /* reset state pointer */

sampleValue = input[n];                         /* get input sample */

state2 = *state_ptr++;                      /*   x[n-2]            */
state1 = *state_ptr--;                      /*   x[n-1]            */
*state_ptr++ = state1;                      /*   x[n-1] -> x[n-2]  */
*state_ptr++ = sampleValue;                 /*   x[n]   -> x[n-1]  */

i = Nsections;
while (--i >= 0)
{
coef_ptr = &((*section_ptr++)->filterCoefficients[0]);  /* point to section filter coefficients */

sampleValue =  *coef_ptr++ * sampleValue;   /*   b0*x[n]           */
sampleValue += *coef_ptr++ * state1;        /*   b1*x[n-1]         */
sampleValue += *coef_ptr++ * state2;        /*   b2*x[n-2]         */

state2 = *state_ptr++;                      /*   y[n-2]            */
state1 = *state_ptr--;                      /*   y[n-1]            */
*state_ptr++ = state1;                      /*   y[n-1] -> y[n-2]  */

sampleValue += *coef_ptr++ * state1;        /*  -a1*y[n-1]         */
sampleValue += *coef_ptr++ * state2;        /*  -a2*y[n-2]         */

*state_ptr++ = sampleValue;                 /*   y[n]   -> y[n-1]  */


/* * at this point sampleValue is the output of the section just completed * and will be the input to the section about to be computed unless this * is the last section, in which case sampleValue is the input. * state1 and state2 are the output states of the section just completed * and will be the input states of of the section about to be computed. */ }

    output[n] = sampleValue;                        /* store output sample */
}

return 0;
}

• I appreciate you taking the time to do this! Unfortunately, I am still having the same problem. I have updated the post with "Edit 3" with the details, and added the code that I am using to call the filter. – Rick Oct 17 '14 at 22:51
• well, i dunno what the problem is. at first it looked like it was synchronized with the sinusoid so i thought it was distortion due to wrap-around overflow somewhere. shouldn't happen with float and it doesn't appear synchronous with the waveform. do you know the number of samples spacing between the glitches? maybe it has to do with your input buffer or how you're moving your samples in and out. – robert bristow-johnson Oct 18 '14 at 1:45
• @Rick, can you modify your program so that it's just a "wire" (so that the output samples are a copy of the input samples) and see if this glitch ends up in the output? – robert bristow-johnson Oct 18 '14 at 4:13
• In my original code, the glitch seemed to happen every (buffer_size / 4) samples, but it does not seem to be the same using your code. I have generated a sampled sine wave and placed it in my code in order to test the filter and not use the input from the codec and there are no glitches in the output, everything sounds clean. I have also noticed that if I try to boost or cut the input from the codec, there is no change in amplitude, but there is when I am filtering generated signal. – Rick Oct 19 '14 at 1:49
• I have also commented out the call to the filter in my code and run everything as is -- input going to ouput-- and the output is clean. No glitches whatsoever. Could this be a sample rate problem? I measure my sample rate at around 48500 Hz, when I'd really like it at 48000. – Rick Oct 19 '14 at 1:52