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