Electronic – IIR filter on a 16bit MCU (PIC24F)

digital-logicfilteriir

Im trying to implement a simple first order IIR filter on a MCU (PIC24FJ32GA002), without success until now.
The filter is a DC tracking filter (low pass filter) whose purpose is to track DC component of a 1.5Hz signal. The difference equation was taken from a TI application note:

y(n)=Kx(n)+y(n-1)(1-K)
with K = 1/2^8

I made a MATLAB script to test it and it works well in the simulation.
Code used:

K=1/2^8
b = K
a = [1 -(1-K)]

Fs=200; // sampling frequency
Ts=1/Fs;
Nx=5000; // number of samples
nT=Ts*(0:Nx-1);
fin=1.5; // signal frequency

randn('state',sum(100*clock));
noise=randn(1,Nx);
noise=noise-mean(noise);

xin=200+9*(cos(2*pi*fin*nT));
xin=xin+noise;

out = filter(b,a,xin);

Simulation
Frequence response

However I can't implement it on a PIC24F microcontroller. i'm representing the coefficients in Q15 (1.15) format, storing them in short variables and using a long one for multiplications. Here it is the code:

short xn;
short y;
short b0 = 128, a1 = 32640; // Q15
long aux1, aux2;

// (...)

while(1){

    xn = readADC(adc_ch);

    aux1 = ((long)b0*xn) << 1;
    aux2 = ((long)a1*y) << 1;
    y = ((aux1 + aux2) >> 16);

    delay_ms(5);
}

Long cast is used to extend the signal so the multiplying operation is done correctly. After each multiplication I shift left one bit to remove the extended signal bit. When summing I shift right 16 bits to get y in Q15 format.

Im debugging the MCU with Pickit2 and "View->Watch" window (MPLAB IDE 8.53) and testing the filter with a DC signal (I change the DC signal with a potenciometer to test different values). The ADC has 10bit resolution and the MCU is supplied with 3.3V.
Some results:

1V –> xn = 312 (correct), yn = 226 (incorrect)
1.5V –> xn = 470 (correct), yn = 228 (completely wrong)

What am I doing wrong?
Any suggestions on how to implement this IIR filter on a 16bit MCU?

Many thanks in advance 🙂

Best Answer

I didn't dive super far into your filter design, but just looking at the source code brings a couple of things up. For example, these lines:

aux1 = ((long)b0*xn) << 1;
aux2 = ((long)a1*y) << 1;
y = ((aux1 + aux2) >> 16);

The first issue I see is the ((long)b0*xn). I have ran across compilers that would compile this incorrectly as ((long)(b0*xn)), which is entirely wrong. Just to be on the safe side, I would write this as (((long)b0)*((long)xn)). To be sure, this is paranoid programming, but...

Next, when you do the "<<1", this is NOT the same as "*2". For most things, it's close, but not for DSP. It has to do with how the MCU/DSP handles overflow conditions and sign extensions, etc. But even if it did work as a *2, you are removing one bit of resolution that you don't need to remove. If you really have to do a *2, then do a *2 and let the compiler figure out if it could substitute a <<1 instead.

The >>16 is also problematic. Off the top of my head, I don't know if it's going to do a logical or arithmetic shift. You want an arithmetic shift. Arithmetic shifts will handle the sign bit correctly where a logical shift will insert zeros for the new bits. Besides, you can save bits of resolution by getting rid of the <<1 and changing the >>16 to >>15. Well, and changing all of these to normal multiplies and divides.

So here's the code I would use:

aux1 = ((long)b0) * ((long)xn);
aux2 = ((long)a1) * ((long)y);
y = (short)((aux1+aux2) / 32768);

Now, I don't claim that this will solve your problem. It may or may not, but it does improve your code.

Related Topic