0
\$\begingroup\$

I have a cascaded IIR filter that has a sampling frequency of 1kHz and consist of a highpass filter with a cutoff frequency of 0.8Hz and a lowpass filter with a cutoff frequency of 5Hz. enter image description here

I have generated a C header and implemented it in some arduino code on a UNO that takes an analog value every millisecond and filters it.

//FilterData
#define M 6
#define N 6
double b[M+1] =  {3.738001453478e-06,1.694065894509e-21,-1.121400436044e-05,8.470329472543e-22,1.121400436044e-05,-1.694065894509e-21,-3.738001453478e-06};
double a[N+1] = {1,   -5.927117642542,    14.63822857612,   -19.28168484046,14.28685989898,   -5.645991366596,   0.9297053744989};

double x[M+1] = {0,0,0,0,0,0,0};
double y[N+1] = {0,0,0,0,0,0,0};

ISR(TIMER1_COMPA_vect){
    int val = analogRead(sensorPin);
    Serial.print("val: ");Serial.println(val);
    val = signalFilter(val);
    Serial.print("Filtered val: ");Serial.println(val);
}

int signalFilter(int invalue){
    Serial.print("invalue=");Serial.println(invalue);
    //FIR
    Serial.print("x=[");
    for(uint32_t k = M; k > 0; k--){
        Serial.print(x[k-1]);Serial.print(", ");
        x[k] = x[k-1];
    } 
    x[0] = (float)invalue;
    Serial.print(x[0]);Serial.println("]");

    //IIR
    Serial.print("y=[");
    for(uint32_t k = N; k > 0; k--){
        Serial.print(y[k-1]);Serial.print(", ");
        y[k] = y[k-1];
    }
    Serial.print("y[0]");Serial.println("]");

    Serial.print("[");
    double FIR = 0; 
    for(uint32_t i = 0; i <= M; i++){//Loop for the sum.
        Serial.print("(");Serial.print(b[i]);Serial.print("*");Serial.print(x[M-i]);Serial.print(")");
        if(i < M){Serial.print("+");}
        FIR += (b[i]*x[M-i]);
    }
    Serial.print("]+");

    Serial.print("[");
    double IIR = 0;
    for(uint32_t j = 1; j <= N; j++){
        Serial.print("(");/*Serial.print((a[j]*y[N-j]));*/Serial.print(a[j]);Serial.print("*");Serial.print(y[N-j]);Serial.print(")");
        if(j < N){Serial.print("+");}
        IIR += (a[j]*y[N-j]);
    }
    Serial.println("]");
    Serial.print("FIR=");Serial.println(FIR);
    Serial.print("IIR=");Serial.println(IIR);
    IIR = FIR+IIR;
    y[0] = IIR;
    Serial.print("y[0]=");Serial.println(y[0]);
    return (int)y[0];
}

But I get a really odd output that eventually overflows.

val: 754
invalue=754
x=[0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 754.00]
y=[0.00, 0.00, 0.00, 0.00, 0.00, 0.00, y[0]]
[(0.00*0.00)+(0.00*0.00)+(-0.00*0.00)+(0.00*0.00)+(0.00*0.00)+(-0.00*0.00)+ (-0.00*754.00)]+[(-5.93*0.00)+(14.64*0.00)+(-19.28*0.00)+(14.29*0.00)+(-5.65*0.00)+(0.93*0.00)]
FIR=-0.00
IIR=0.00
y[0]=-0.00
Filtered val: 0
val: 972
invalue=972
x=[0.00, 0.00, 0.00, 0.00, 0.00, 754.00, 972.00]
y=[0.00, 0.00, 0.00, 0.00, 0.00, -0.00, y[0]]
[(0.00*0.00)+(0.00*0.00)+(-0.00*0.00)+(0.00*0.00)+(0.00*0.00)+(-0.00*754.00)+(-0.00*972.00)]+[(-5.93*0.00)+(14.64*0.00)+(-19.28*0.00)+(14.29*0.00)+(-5.65*-0.00)+(0.93*-0.00)]
FIR=-0.00
IIR=0.01
y[0]=0.01
Filtered val: 0
val: 974
invalue=974
x=[0.00, 0.00, 0.00, 0.00, 754.00, 972.00, 974.00]
y=[0.00, 0.00, 0.00, 0.00, -0.00, 0.01, y[0]]
[(0.00*0.00)+(0.00*0.00)+(-0.00*0.00)+(0.00*0.00)+(0.00*754.00)+(-0.00*972.00)+(-0.00*974.00)]+[(-5.93*0.00)+(14.64*0.00)+(-19.28*0.00)+(14.29*-0.00)+(-5.65*0.01)+(0.93*0.01)]
FIR=0.00
IIR=-0.09
y[0]=-0.08
Filtered val: 0
val: 971
invalue=971
x=[0.00, 0.00, 0.00, 754.00, 972.00, 974.00, 971.00]
y=[0.00, 0.00, 0.00, -0.00, 0.01, -0.08, y[0]]
[(0.00*0.00)+(0.00*0.00)+(-0.00*0.00)+(0.00*754.00)+(0.00*972.00)+(-0.00*974.00)+(-0.00*971.00)]+[(-5.93*0.00)+(14.64*0.00)+(-19.28*-0.00)+(14.29*0.01)+(-5.65*-0.08)+(0.93*-0.08)]
FIR=0.01
IIR=0.57
y[0]=0.58
Filtered val: 0
val: 634
invalue=634
x=[0.00, 0.00, 754.00, 972.00, 974.00, 971.00, 634.00]
y=[0.00, 0.00, -0.00, 0.01, -0.08, 0.58, y[0]]
[(0.00*0.00)+(0.00*0.00)+(-0.00*754.00)+(0.00*972.00)+(0.00*974.00)+ (-0.00*971.00)+(-0.00*634.00)]+[(-5.93*0.00)+(14.64*-0.00)+(-19.28*0.01)+(14.29*-0.08)+(-5.65*0.58)+(0.93*0.58)]
FIR=0.00
IIR=-4.13
y[0]=-4.13
Filtered val: -4
val: 531
invalue=531
x=[0.00, 754.00, 972.00, 974.00, 971.00, 634.00, 531.00]
y=[0.00, -0.00, 0.01, -0.08, 0.58, -4.13, y[0]]
[(0.00*0.00)+(0.00*754.00)+(-0.00*972.00)+(0.00*974.00)+(0.00*971.00)+(-0.00*634.00)+(-0.00*531.00)]+[(-5.93*-0.00)+(14.64*0.01)+(-19.28*-0.08)+(14.29*0.58)+(-5.65*-4.13)+(0.93*-4.13)]
FIR=-0.00
IIR=29.50
y[0]=29.50
Filtered val: 29
val: 215
invalue=215
x=[754.00, 972.00, 974.00, 971.00, 634.00, 531.00, 215.00]
y=[-0.00, 0.01, -0.08, 0.58, -4.13, 29.50, y[0]]
[(0.00*754.00)+(0.00*972.00)+(-0.00*974.00)+(0.00*971.00)+(0.00*634.00)+(-0.00*531.00)+(-0.00*215.00)]+[(-5.93*0.01)+(14.64*-0.08)+(-19.28*0.58)+(14.29*-4.13)+(-5.65*29.50)+(0.93*29.50)]
FIR=-0.00
IIR=-210.56
y[0]=-210.56
Filtered val: -210
val: 599
invalue=599
x=[972.00, 974.00, 971.00, 634.00, 531.00, 215.00, 599.00]
y=[0.01, -0.08, 0.58, -4.13, 29.50, -210.56, y[0]]
[(0.00*972.00)+(0.00*974.00)+(-0.00*971.00)+(0.00*634.00)+(0.00*531.00)+(-0.00*215.00)+(-0.00*599.00)]+[(-5.93*-0.08)+(14.64*0.58)+(-19.28*-4.13)+(14.29*29.50)+(-5.65*-210.56)+(0.93*-210.56)]
FIR=-0.00
IIR=1503.07
y[0]=1503.07
Filtered val: 1503
val: 469
invalue=469
x=[974.00, 971.00, 634.00, 531.00, 215.00, 599.00, 469.00]
y=[-0.08, 0.58, -4.13, 29.50, -210.56, 1503.07, y[0]]
[(0.00*974.00)+(0.00*971.00)+(-0.00*634.00)+(0.00*531.00)+(0.00*215.00)+(-0.00*599.00)+(-0.00*469.00)]+[(-5.93*0.58)+(14.64*-4.13)+(-19.28*29.50)+(14.29*-210.56)+(-5.65*1503.07)+(0.93*1503.07)]
FIR=-0.00
IIR=-10729.75
y[0]=-10729.75
Filtered val: -10729
val: 541
invalue=541
x=[971.00, 634.00, 531.00, 215.00, 599.00, 469.00, 541.00]
y=[0.58, -4.13, 29.50, -210.56, 1503.07, -10729.75, y[0]]
[(0.00*971.00)+(0.00*634.00)+(-0.00*531.00)+(0.00*215.00)+(0.00*599.00)+(-0.00*469.00)+(-0.00*541.00)]+[(-5.93*-4.13)+(14.64*29.50)+(-19.28*-210.56)+(14.29*1503.07)+(-5.65*-10729.75)+(0.93*-10729.75)]
FIR=0.00
IIR=76594.85
y[0]=76594.85
Filtered val: 11058

Is there anyone that could tell me whats wrong?

\$\endgroup\$
  • 2
    \$\begingroup\$ Verify your filter with a known and verifiable sequence first before you connect your adc. \$\endgroup\$ – Humpawumpa Mar 20 at 10:18
  • \$\begingroup\$ Possible migrate to signal processing stack exchange? \$\endgroup\$ – loudnoises Mar 20 at 10:20
  • \$\begingroup\$ Try debugging it in excel. It's a useful tool to use because you can implement all sorts of algorithm modifications quite simply. I'm sure there are other tools that you can use to debug this problem. \$\endgroup\$ – Andy aka Mar 20 at 11:10
  • 4
    \$\begingroup\$ You need to fix the fundamental brokenness of doing Serial.print inside an ISR before you even start on trying to debug anything else. \$\endgroup\$ – brhans Mar 20 at 11:13
  • \$\begingroup\$ See this. \$\endgroup\$ – a concerned citizen Mar 20 at 18:15
2
\$\begingroup\$

I haven’t looked too closely at your code, but there is an obvious issue with it.

You are implementing your 6th order filter in a direct form, not as biquads (I.e., second order sections).

Although sixth order is really pushing the limit of numeric representation in a computer (it might still work for some pole combinations but not for others), this problem is further compounded by the lower resolution floating point number representations in an Arduino.

It’s never a good idea to represent any digital filter of this order as a single unit, as the magnitude of the poles lead to coefficients with very large magnitude disparities. When placed together in a single addition such disparities lead to large roundoff errors that in practice move the coefficients around and in a non-linear manner. Besides a seriously distorted response, this leads to placing some coefficients outside of the unit circle and thus to instability.

The solution is to break the filter into a sequence of smaller sections (normally biquads) by grouping similar magnitude poles and zeros together. This produces less roundoff error and thus fewer problems.

\$\endgroup\$
2
\$\begingroup\$

My bet would be numerical precision in the HPF or the probable brokenness that is serial_printf in an interrupt handler, I don't know the platform, but code that might block in an ISR?

On a UNO the double type is really a 4 byte float, not the 8 byte one you might expect if coming from a PC or grown up processor, so you only have 24 bits of precision and a range of about about $$10^{\pm38}$$

IIR filters generally become more demanding of numerical precision as the ratio of corner sample rate to frequency increases, and you have some very small coefficients in play $$10^{-22}$$ is not a number you ever really want to see in a filter kernel using what is effectively a float type.

I might be missing something among all the printfs but I don't see the IIR filter code itself there, I would expect a couple of biquads of some kind? The code also has comments that mention FIRs, which are a whole different animal.

Does your filter do the right thing if you run it as a matlab or octave script? How about if you run it on test data on a PC (Using a 4 byte float)?

Do you have any poles really close to the unit circle (Use the actual 32 bit floats when working this out, not the theoretical numbers)?

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.