Trouble with fixed-point arithmetic 2nd order IIR implementation in c

Question

DSP newcomer here!

I am tinkering around with a TI DSP and am trying to implement a second order IIR filter in C. Input is 16 bit 2's complement, as is the output, the accumulator is 32 bit wide. I have tried several structures: Direct form I, Transposed direct form II, Gold & Rader structure. However, I can't seem to get anything to work and I am not sure whether it is a typo in the code of all of my filter implementations or whether it is because of me misunderstanding the principle of fixed-point arithmetic. I assume - or better: hope - that it's the same mistake that is affecting all implementations.

Here's my code for the transposed df II implementation:

// state variables
int w[3] = {0};

short filter_dfii(short value) {
    w[0] = ((value * sos_filter_coeffs_b[0] + w[1]) >> 15 ) & 0xffff;
    w[1] = value * sos_filter_coeffs_b[1]
       - w[0] * sos_filter_coeffs_a[1]
       + w[2];
    w[2] = value * sos_filter_coeffs_b[2]
       - w[0] * sos_filter_coeffs_a[2];

    return (short) w[0];
}

My df I implementation looks like this:

int last_value[2] = {0};

short filter_dfi(short value) {
    w[0] =
        ((value * sos_filter_coeffs_b[0]
        + last_value[0] * sos_filter_coeffs_b[1]
        + last_value[1] * sos_filter_coeffs_b[2]
        - w[1] * sos_filter_coeffs_a[1]
        - w[2] * sos_filter_coeffs_a[2])
        >> 15) & 0xffff;
    last_value[1] = last_value[0];
    last_value[0] = value;
    w[2] = w[1];
    w[1] = w[0];
    return (short) w[0];
}

The filter coefficients are:

int sos_filter_coefficients_a = {32767, -51150, 21015};
int sos_filter_coefficients_b = {658, 1316, 658};

How did I get those coefficients? I utilized the Matlab 'butter' function and specified a 2nd order Butterworth lowpass with a cutoff frequency of 0.05*fs or 2400 kHz in my specific case. Then I multiplied each coefficient with 32767 and rounded it to an integer value. If I call Matlab freqz function with those quantized coefficients, I get a bode plot with the expected result, so I believe quantization error is not the issue here.

However, the result is not a lowpass but just garbage. Unfortunately, I lack the proper testing equipment to describe the result scientifically, but it sounds like digital oscillation that is slightly modulated by the actual input signal and is present even when there's no input at all.

So, in my opinion, this can't be an unstable filter since the poles are well within the unit circle. It can't be overflow limit cycle either because then it shouldn't be happening when there's only low input level or no input at all.. right?

If I alter the coefficients to

int sos_filter_coefficients_a = {32767, 0, 0};
int sos_filter_coefficients_b = {32767, 0, 0};

input gets returned basically unprocessed as expected, so my issues most likely aren't caused by any external factors - there must be something wrong with my coefficients, with the function itself, or both.

Any help is much appreciated!

Answer 1

The " & 0xffff" operation is questionable. It works for positive numbers but for negative numbers it will mask out the sign bit as well, so you will turn everything into positive numbers.

For the DF1, it would be better to delay the shift and mask until you return the actual output and do the conversion to short. This way you can keep the state variable in full double precision.

Answer 2

For starters, what size is "int" ? 16 bits or 32 bits? One of your coefficient is "-51150" which requires 17 bits to represent, hence the need for a 32-bit integer for coefficients.

Make sure "int" is indeed a 32-bit integer, otherwise it would not work. Or, you can maybe multiply coefficients by 16384 instead of 32768 in order to make them fit in 16 bits.