Strange Results from Biquad Peaking EQ Filter

Question

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

Answer 1

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