Go Down

Topic: A FIR Finite Impuls Response filter routine (Read 6 times) previous topic - next topic

René Knuvers

Jul 14, 2010, 06:02 pm Last Edit: Jul 14, 2010, 06:09 pm by reneknuvers Reason: 1
Hello there, fellow Arduino Programmers.

I spent my day (well, part of it at least) programming a FIR (or Finite Impuls Response) filter on the Arduino.

It can be used as a building block for any program that requires some form of filtering. Because of its simple nature, you may use this for various purposes:
- filtering audio (duh...)
- smoothing of analog inputs
- smoothing of PWM outputs (ramp up, and down of motors, gradually turning lights/LED's on and off)
- taking a 'running average' on streaming data
- removing 'DC' or static components of a signal (requires a 'high pass'  or 'band pass' filter

I've written this to sharpen my limited C-programming skills, so there is probably lots of room for improvement. I hope people will criticize my code, so I can learn to write proper code, and optimize it. Once optimized, I hope this code could make it into the 'playground'. I guess if it handled INT's it would be more useful there.

I've commented the code extensively, so people could understand what it going on, and change it to their own needs.

Code: [Select]

/*

Simple FIR filter implementation

This routine is just written as an excersize in Arduino programming.
The routine implements a Finite Impuls Response filter.

Finite Impuls Response filters can be used to create all kinds of filters.
Filters can be used to partially remove noise, remove 'DC' offsets, or 'shape' analogue values. See http://en.wikipedia.org/wiki/Finite_impulse_response to read more aboute FIR filters. Several online sources are available to generate the coefficients. Google for the various sources.

About the routine: it is not really optimized for anything. It accepts exactly one float, spits out one float, and should be called for every sample. For the filter to work as expected, the samples should be 'equidistant' in time, so taken at identical intervals.

Making this an 'integer' function (to make it more usable with the onboard analog inputs) can be done. This would also limit execution time and memory usage.

Talking performance: using floats a 5 tap cycle takes about 180µs, so that is pretty fast and allows about 3~4kHz sample rate, taking some overhead in the main loop into consideration. Using ints or even better bytes should fraction the looptime by at least 4 to 8 times!

Rene Knuvers, 14JUL2010, arduino@reneknuvers.nl

*/

#define FILTERTAPS 5            // The number of taps in your filter, or better the number of coefficients

// declare array for values
float values[FILTERTAPS] = {0, 0, 0, 0, 0}; // to have a nice start up, fill the array with 0's
// NOTE: this could be 1 shorter, as we already have the input in a variable. This would save stack size/memory

// declare input and output variables
float input = 1; // without a real input, looking at the step respons (input at unity, 1) would be nice to see
float output = 0; // output as a 0, but that doesn't really matter

// our counter, real life applications won't need this
byte n = 0;

void setup() {
     Serial.begin(115200);      // open the serial port, I like them fast ;-)
}

float fir(float in){
     static byte k;                  // k stores a pointer to create a circular memory through the array. This variable remains its value the next time we call the fir routine (hence 'static')
     byte i = 0;                        // i is a counter to step through the filter coefficients and stored values
     float out = 0;                        // out is the return variable. It is set to 0 every time we call the filter!

     values[k] = input;                  // store the input in the array at the current pointer position

     
     // declare variables for coefficients
       // these should be calculated by hand, or using a tool
     // in case a phase linear filter is required, the coefficients are symmetric
     // for time optimization it seems best to enter symmetric values like below
     float coef[FILTERTAPS] = { 0.021, 0.096, 0.146, 0.096, 0.021};

     //declare gain coefficient to scale the output back to normal
     float gain = 0.38; // set to 1 and input unity to see what this needs to be

     
     for (i=0; i<FILTERTAPS; i++) {            // we step through the array
           out += coef[i] * values[(i + k) % FILTERTAPS];      // ... and add and multiply each value to accumulate the output
                                   //  (i + k) % FILTERTAPS creates a cyclic way of getting through the array
       }
     out /= gain;                        // We need to scale the output (unless the coefficients provide unity gain in the passband)

     k = (k+1) % FILTERTAPS;                  // k is increased and wraps around the FILTERTAPS, so next time we will overwrite the oldest saved sample in the array
     
       return out;                        // we send the output value back to whoever called the routine

}


void loop() {

     // This is the loop that takes care of calling the FIR filter for some samples

     for (n = 0; n < FILTERTAPS + 2; n++) {             // If you like to see the step response, take at least as many cycles as the length of your FIR filter (FILTERTAPS + 1 or 2)
           Serial.print("n= ");            // print the sample number
           Serial.println(n, DEC);
           Serial.println("Now calling fir...");
           output = fir(input);            // here we call the fir routine with the input. The value 'fir' spits out is stored in the output variable.
           Serial.print("fir presented the following value= ");
           Serial.println(output);            // just for debugging or to understand what it does, print the output value
     }
     while (true) {};                  // endless loop
}


Please, feel free to comment and post alternatives!

Rene Knuvers
Goirle, The Netherlands
arduino@reneknuvers.nl

davekw7x

#1
Jul 14, 2010, 06:29 pm Last Edit: Jul 14, 2010, 06:39 pm by davekw7x Reason: 1
I have done any testing, but I like it (so far).

As a matter of style, I'm sure this is a (minor) oversight:

Your fir function has a parameter in, but it actually works on the file scope variable input.  Since your main loop is expecting it to work on input the results are OK, but when you get ready to transfer the function to a library, you will have to make sure it works in its parameter.  (Of course a library function would (probably) use a pointer as a parameter for the array of taps and have a parameter for gain.  Etc.

I expect to have some fun with this later.  Maybe with more taps. Thanks for sharing.


Regards,

Dave

René Knuvers

Thanks for the reply. I guess I get what you are saying. Pointers, classes, constructors and deconstructors is a bit out of my league right now, but I will try to get that to work at some point.

Your remark about 'in' and 'input': you are pointing at the use of 'input' within the fir-routine? That's a typo, a leftover of the time when I changed the name to in for the routine. Thanks for finding that!

René Knuvers

The corrected code:

Code: [Select]

/*

Simple FIR filter implementation

This routine is just written as an excersize in Arduino programming.
The routine implements a Finite Impuls Response filter.

Finite Impuls Response filters can be used to create all kinds of filters.
Filters can be used to partially remove noise, remove 'DC' offsets, or 'shape' analogue values. See http://en.wikipedia.org/wiki/Finite_impulse_response to read more aboute FIR filters. Several online sources are available to generate the coefficients. Google for the various sources.

About the routine: it is not really optimized for anything. It accepts exactly one float, spits out one float, and should be called for every sample. For the filter to work as expected, the samples should be 'equidistant' in time, so taken at identical intervals.

Making this an 'integer' function (to make it more usable with the onboard analog inputs) can be done. This would also limit execution time and memory usage.

Talking performance: using floats a 5 tap cycle takes about 180µs, so that is pretty fast and allows about 3~4kHz sample rate, taking some overhead in the main loop into consideration. Using ints or even better bytes should fraction the looptime by at least 4 to 8 times!

Rene Knuvers, 14JUL2010, arduino@reneknuvers.nl

*/

#define FILTERTAPS 5            // The number of taps in your filter, or better the number of coefficients

// declare array for values
float values[FILTERTAPS] = {0, 0, 0, 0, 0}; // to have a nice start up, fill the array with 0's
// NOTE: this could be 1 shorter, as we already have the input in a variable. This would save stack size/memory

// declare input and output variables
float input = 1; // without a real input, looking at the step respons (input at unity, 1) would be nice to see
float output = 0; // output as a 0, but that doesn't really matter

// our counter, real life applications won't need this
byte n = 0;

void setup() {
     Serial.begin(115200);      // open the serial port, I like them fast ;-)
}

float fir(float in){
     static byte k;                  // k stores a pointer to create a circular memory through the array. This variable remains its value the next time we call the fir routine (hence 'static')
     byte i = 0;                        // i is a counter to step through the filter coefficients and stored values
     float out = 0;                        // out is the return variable. It is set to 0 every time we call the filter!

     values[k] = in;                        // store the input of the routine (contents of the 'in' variable) in the array at the current pointer position

     
     // declare variables for coefficients
       // these should be calculated by hand, or using a tool
     // in case a phase linear filter is required, the coefficients are symmetric
     // for time optimization it seems best to enter symmetric values like below
     float coef[FILTERTAPS] = { 0.021, 0.096, 0.146, 0.096, 0.021};

     //declare gain coefficient to scale the output back to normal
     float gain = 0.38; // set to 1 and input unity to see what this needs to be

     
     for (i=0; i<FILTERTAPS; i++) {            // we step through the array
           out += coef[i] * values[(i + k) % FILTERTAPS];      // ... and add and multiply each value to accumulate the output
                                   //  (i + k) % FILTERTAPS creates a cyclic way of getting through the array
       }
     out /= gain;                        // We need to scale the output (unless the coefficients provide unity gain in the passband)

     k = (k+1) % FILTERTAPS;                  // k is increased and wraps around the FILTERTAPS, so next time we will overwrite the oldest saved sample in the array
     
       return out;                        // we send the output value back to whoever called the routine

}


void loop() {

     // This is the loop that takes care of calling the FIR filter for some samples

     for (n = 0; n < FILTERTAPS + 2; n++) {             // If you like to see the step response, take at least as many cycles as the length of your FIR filter (FILTERTAPS + 1 or 2)
           Serial.print("n= ");            // print the sample number
           Serial.println(n, DEC);
           Serial.println("Now calling fir...");
           output = fir(input);            // here we call the fir routine with the input. The value 'fir' spits out is stored in the output variable.
           Serial.print("fir presented the following value= ");
           Serial.println(output);            // just for debugging or to understand what it does, print the output value
     }
     while (true) {};                  // endless loop
}

AlphaBeta

You do not set k to a value, but you use it as an indexer, this may cause some problems.

René Knuvers

Ah, of course, when it has a value that is greater than the array length... right?

René Knuvers

And this code initiates the indexer [k] at 0:

Code: [Select]

/*

Simple FIR filter implementation

This routine is just written as an excersize in Arduino programming.
The routine implements a Finite Impuls Response filter.

Finite Impuls Response filters can be used to create all kinds of filters.
Filters can be used to partially remove noise, remove 'DC' offsets, or 'shape' analogue values. See http://en.wikipedia.org/wiki/Finite_impulse_response to read more aboute FIR filters. Several online sources are available to generate the coefficients. Google for the various sources.

About the routine: it is not really optimized for anything. It accepts exactly one float, spits out one float, and should be called for every sample. For the filter to work as expected, the samples should be 'equidistant' in time, so taken at identical intervals.

Making this an 'integer' function (to make it more usable with the onboard analog inputs) can be done. This would also limit execution time and memory usage.

Talking performance: using floats a 5 tap cycle takes about 180µs, so that is pretty fast and allows about 3~4kHz sample rate, taking some overhead in the main loop into consideration. Using ints or even better bytes should fraction the looptime by at least 4 to 8 times!

Rene Knuvers, 14JUL2010, arduino@reneknuvers.nl

*/

#define FILTERTAPS 5            // The number of taps in your filter, or better the number of coefficients

// declare array for values
float values[FILTERTAPS] = {0, 0, 0, 0, 0}; // to have a nice start up, fill the array with 0's
// NOTE: this could be 1 shorter, as we already have the input in a variable. This would save stack size/memory

// declare input and output variables
float input = 1; // without a real input, looking at the step respons (input at unity, 1) would be nice to see
float output = 0; // output as a 0, but that doesn't really matter

// our counter, real life applications won't need this
byte n = 0;

void setup() {
     Serial.begin(115200);      // open the serial port, I like them fast ;-)
}

float fir(float in){
     static byte k = 0;                  // k stores a pointer to create a circular memory through the array. This variable remains its value the next time we call the fir routine (hence 'static')
     byte i = 0;                        // i is a counter to step through the filter coefficients and stored values
     float out = 0;                        // out is the return variable. It is set to 0 every time we call the filter!

     values[k] = in;                        // store the input of the routine (contents of the 'in' variable) in the array at the current pointer position

     
     // declare variables for coefficients
       // these should be calculated by hand, or using a tool
     // in case a phase linear filter is required, the coefficients are symmetric
     // for time optimization it seems best to enter symmetric values like below
     float coef[FILTERTAPS] = { 0.021, 0.096, 0.146, 0.096, 0.021};

     //declare gain coefficient to scale the output back to normal
     float gain = 0.38; // set to 1 and input unity to see what this needs to be

     
     for (i=0; i<FILTERTAPS; i++) {            // we step through the array
           out += coef[i] * values[(i + k) % FILTERTAPS];      // ... and add and multiply each value to accumulate the output
                                   //  (i + k) % FILTERTAPS creates a cyclic way of getting through the array
       }
     out /= gain;                        // We need to scale the output (unless the coefficients provide unity gain in the passband)

     k = (k+1) % FILTERTAPS;                  // k is increased and wraps around the FILTERTAPS, so next time we will overwrite the oldest saved sample in the array
     
       return out;                        // we send the output value back to whoever called the routine

}


void loop() {

     // This is the loop that takes care of calling the FIR filter for some samples

     for (n = 0; n < FILTERTAPS + 2; n++) {             // If you like to see the step response, take at least as many cycles as the length of your FIR filter (FILTERTAPS + 1 or 2)
           Serial.print("n= ");            // print the sample number
           Serial.println(n, DEC);
           Serial.println("Now calling fir...");
           output = fir(input);            // here we call the fir routine with the input. The value 'fir' spits out is stored in the output variable.
           Serial.print("fir presented the following value= ");
           Serial.println(output);            // just for debugging or to understand what it does, print the output value
     }
     while (true) {};                  // endless loop
}

davekw7x

Quote
You do not set k to a value

Variables with static storage duration are always initialized to zero (at the time the program is loaded), so it's not really necessary for the program to initialize it to zero.

However, as a matter of style (for the benefit of Humans reading the code), many people do show the initialization explicitly, and it's certainly not wrong to do so.

Regards,

Dave

AlphaBeta

Made a libraryesque class of this, use it if you want in any way you want.
Code: [Select]
#define FILTERTAPS 5            // The number of taps in your filter, or better the number of coefficients

/*
 LIBRARY CODE STARTS HERE
 
 http://www.arduino.cc/cgi-bin/yabb2/YaBB.pl?num=1279123369
 
 Made into a library by AlphaBeta (www.AlexanderBrevig.com)
 
 Do whatever you want with it.
*/
template<int filterTaps>
class FIR {
     public:
           //construct without coefs
           FIR() {
                 k = 0; //initialize so that we start to read at index 0
                 for (int i=0; i<filterTaps; i++) {      
                       values[i] = 0; // to have a nice start up, fill the array with 0's
                 }
                 //TODO calculate default gain?
                 //TODO calculate default coefs?
           }
           //construct with coefs
           FIR(float newGain, float *newCoefs) {
                 k = 0; //initialize so that we start to read at index 0
                 setGain(newGain);
                 for (int i=0; i<filterTaps; i++) {      
                       values[i] = 0; // to have a nice start up, fill the array with 0's
                 }
                 setCoefficients(newCoefs);
           }
           
           void setGain(float newGain) {gain = newGain;}

           void setCoefficients(float *newCoefs) {
                 for (int i=0; i<filterTaps; i++) {      
                       coef[i] = newCoefs[i];
                 }
           }
           //set coefficient at specified index
           void setCoefficient(int idx, float newCoef) { coef[idx] = newCoef; }
           
           float process(float in) {
                 float out = 0;                        // out is the return variable. It is set to 0 every time we call the filter!

                 values[k] = in;                        // store the input of the routine (contents of the 'in' variable) in the array at the current pointer position

                 for (int i=0; i<filterTaps; i++) {                              // we step through the array
                       out += coef[i] * values[(i + k) % filterTaps];      // ... and add and multiply each value to accumulate the output
                                                                                               //  (i + k) % filterTaps creates a cyclic way of getting through the array
                 }
                 
                 out /= gain;                        // We need to scale the output (unless the coefficients provide unity gain in the passband)

                 k = (k+1) % filterTaps;            // k is increased and wraps around the filterTaps, so next time we will overwrite the oldest saved sample in the array

                 return out;                              // we send the output value back to whoever called the routine
           }
           
     private:
           float values[filterTaps];
           
           float coef[filterTaps];
           
           //declare gain coefficient to scale the output back to normal
           float gain; // set to 1 and input unity to see what this needs to be
           
           int k; // k stores the index of the current array read to create a circular memory through the array
};
/*
 LIBRARY CODE ENDS HERE
*/


FIR<FILTERTAPS> fir;

void setup() {

     Serial.begin(115200);      // open the serial port, I like them fast ;-)
     
     // declare variables for coefficients
     // these should be calculated by hand, or using a tool
     // in case a phase linear filter is required, the coefficients are symmetric
     // for time optimization it seems best to enter symmetric values like below
     float coef[FILTERTAPS] = { 0.021, 0.096, 0.146, 0.096, 0.021};
     fir.setCoefficients(coef);

       //declare gain coefficient to scale the output back to normal
     float gain = 0.38; // set to 1 and input unity to see what this needs to be
     fir.setGain(gain);
}

void loop() {

     // declare input and output variables
     float input = 1; // without a real input, looking at the step respons (input at unity, 1) would be nice to see
     float output = 0; // output as a 0, but that doesn't really matter

     // This is the loop that takes care of calling the FIR filter for some samples

     for (byte n = 0; n < FILTERTAPS + 2; n++) {             // If you like to see the step response, take at least as many cycles as the length of your FIR filter (FILTERTAPS + 1 or 2)
           Serial.print("n= ");            // print the sample number
           Serial.println(n, DEC);
           Serial.println("Now calling fir...");
           output = fir.process(input);            // here we call the fir routine with the input. The value 'fir' spits out is stored in the output variable.
           Serial.print("fir presented the following value= ");
           Serial.println(output);            // just for debugging or to understand what it does, print the output value
     }
     
     while (true) {};                  // endless loop
}

René Knuvers

Ah, that looks a lot more professional. I sort of get the idea here, but I will print it out some time this weekend to see what exactly you did. Do you mind me asking for some explanation at a later time?

I see you use an integer for the indexes. Is this generally a good idea? I am wondering if this doesn't take too many processor cycles and a lot of register swapping on an AVR (I don't know C, but do know something about assembler language)

I guess a byte would be enough, since even for a demanding filter with impressive spec's, 255 steps would be enough. If you need more, you are either looking into the wrong filter design (look for an IIR filter) or doing something that ought not to be done on an Arduino (like detecting stars in a multi-megapixel picture of the sky...)

AlphaBeta

Feel free to contact me by PM or by mail :)

Talk to you later!

Honk

Hello!

This is just what I need for a quadcopter project measuring air pressure and IR distance sensors for height holding/regulating. I was gonna write it myself this evening, but this seems very much better than I suspect I could write now.

I wonder how to add/use it to/from a library? I don't see an #include in what I believe you made the "exapmple file"?

Thanks,
Henrik

René Knuvers

Hello Henrik,

I hope this works for you. There is many info on the Arduino website about putting routines in a library. Basically you take the routines that are in this post and put them in the appropriate files for the library. Take a closer look at http://www.arduino.cc/en/Hacking/LibraryTutorial to see what I mean.

Height control from different sources... hmm, that sounds like a job better done with a Kalman filter, or a PID regulator! Take a look at http://www.arduino.cc/cgi-bin/yabb2/YaBB.pl?num=1225283209 and http://starlino.com/imu_kalman_arduino.html

This is why: a filter helps with getting the data evened out, or removing stuff you don't need (noise, high frequencies, low frequencies, hum, etc) but it doesn't add relevance to the data. The Kalman technique does actually that. It takes into consideration what a sensor is good at (a barometer is OK for measuring long-time values that relate to height), and what it's not (a barometer will get noisy from external wind sources, movement of the copter, etc) Look for similarities with Gyro's and accelero's!

PID is the easiest source for stable control of height, although figuring out how to parametrize your PID regulator could be tricky. But in the end, a simple P action, and a little D would help you a lot getting things stable. Drift over time, or going to 'exact x meters' isn't relevant to you (or is it?) so the I action can be left out without a problem. Construction of your E (error) signal is the problem here. I would start off with simply take the set value (from the Remote), substract (0.5 × IR height estimation, and 0.5 × barometer height estimation) You can simply scale both to see what it does. You COULD use a filter (FIR low-pass) to remove HF noise from the barometer, and probably from the IR too. I wouldn't filter too much, though. Probably the best filter there, would be a simple RC filter on the analog outputs!

Go Up