FIR Implementation

Projects Programming
1

Hello,

I am trying to implement an FIR, using the Arduino FIR Library. I first designed the FIR filter in MATLAB's Filter Design and Analysis Toolbox (fdatool), then validated it in Simulink. I am now trying to implement it on an Arduino Due to filter real-time data.

The filter is a 48-tap FIR low-pass filter sampled at 20 Hz (I know its slow, but its apart of my application). I have a 20 Hz timer calling the function to calculate the FIR. Majority of the following code is just to generate a test waveform that contains 2 frequencies (aptly named 'FASTER' and 'SLOWER'). The FIR is designed to attenuate the FASTER frequency.

#define FASTER_FREQUENCY         75 // BPM
#define FASTER_AMPLITUDE         12 
#define SLOWER_FREQUENCY     12 // BPM
#define SLOWER_AMPLITUDE     20 
#define OFFSET                    0 

#include <stdio.h>

#include <FIR.h>
#define FILTERTAPS 48
FIR fir;

double FASTER_motion [] = {
  0.0001234098, 0.0002226299, 0.0003936690, 0.0006823281, 0.0011592292, 0.0019304541, 0.0031511116, 0.0050417603, 0.0079070541, 0.0121551783, 0.0183156389, 0.0270518469, 0.0391638951, 0.0555762126, 0.0773047404, 0.1053992246, 0.1408584209, 0.1845195240, 0.2369277587, 0.2981972794, 0.3678794412, 0.4448580662, 0.5272924240, 0.6126263942, 0.6976763261, 0.7788007831, 0.8521437890, 0.9139311853, 0.9607894392, 0.9900498337, 1.0000000000, 0.9900498337, 0.9607894392, 0.9139311853, 0.8521437890, 0.7788007831, 0.6976763261, 0.6126263942, 0.5272924240, 0.4448580662, 0.3678794412, 0.2981972794, 0.2369277587, 0.1845195240, 0.1408584209, 0.1053992246, 0.0773047404, 0.0555762126, 0.0391638951, 0.0270518469, 0.0183156389, 0.0121551783, 0.0079070541, 0.0050417603, 0.0031511116, 0.0019304541, 0.0011592292, 0.0006823281, 0.0003936690, 0.0002226299, 0.0001234098
};
int FASTER_frequency(FASTER_FREQUENCY), FASTER_period(0), FASTER_amplitude(FASTER_AMPLITUDE), FASTER_index(0);
int FASTER_size = sizeof(FASTER_motion) / sizeof(FASTER_motion[0]);
boolean FASTER_state = false;

double SLOWER_motion [] = {
  0.0019305, 0.0020976, 0.002278, 0.0024726, 0.0026822, 0.002908, 0.0031511, 0.0034126, 0.0036938, 0.0039958, 0.0043202, 0.0046684, 0.0050418, 0.005442, 0.0058707, 0.0063297, 0.0068208, 0.0073459, 0.0079071, 0.0085063, 0.0091459, 0.0098282, 0.010555, 0.01133, 0.012155, 0.013033, 0.013966, 0.014958, 0.016012, 0.01713, 0.018316, 0.019573, 0.020905, 0.022315, 0.023807, 0.025385, 0.027052, 0.028813, 0.030671, 0.032631, 0.034697, 0.036873, 0.039164, 0.041574, 0.044108, 0.046771, 0.049566, 0.0525, 0.055576, 0.0588, 0.062177, 0.06571, 0.069406, 0.073269, 0.077305, 0.081517, 0.085911, 0.090491, 0.095263, 0.10023, 0.1054, 0.11077, 0.11635, 0.12215, 0.12816, 0.1344, 0.14086, 0.14755, 0.15447, 0.16162, 0.16901, 0.17665, 0.18452, 0.19264, 0.201, 0.20961, 0.21847, 0.22757, 0.23693, 0.24653, 0.25638, 0.26647, 0.2768, 0.28738, 0.2982, 0.30925, 0.32053, 0.33204, 0.34377, 0.35572, 0.36788, 0.38024, 0.3928, 0.40555, 0.41849, 0.43159, 0.44486, 0.45828, 0.47184, 0.48554, 0.49935, 0.51327, 0.52729, 0.54139, 0.55556, 0.56978, 0.58404, 0.59833, 0.61263, 0.62691, 0.64118, 0.65541, 0.66958, 0.68367, 0.69768, 0.71157, 0.72534, 0.73897, 0.75243, 0.76572, 0.7788, 0.79167, 0.8043, 0.81669, 0.8288, 0.84062, 0.85214, 0.86334, 0.8742, 0.88471, 0.89484, 0.90459, 0.91393, 0.92286, 0.93136, 0.93941, 0.94701, 0.95414, 0.96079, 0.96695, 0.9726, 0.97775, 0.98238, 0.98648, 0.99005, 0.99308, 0.99557, 0.9975, 0.99889, 0.99972, 1, 0.99972, 0.99889, 0.9975, 0.99557, 0.99308, 0.99005, 0.98648, 0.98238, 0.97775, 0.9726, 0.96695, 0.96079, 0.95414, 0.94701, 0.93941, 0.93136, 0.92286, 0.91393, 0.90459, 0.89484, 0.88471, 0.8742, 0.86334, 0.85214, 0.84062, 0.8288, 0.81669, 0.8043, 0.79167, 0.7788, 0.76572, 0.75243, 0.73897, 0.72534, 0.71157, 0.69768, 0.68367, 0.66958, 0.65541, 0.64118, 0.62691, 0.61263, 0.59833, 0.58404, 0.56978, 0.55556, 0.54139, 0.52729, 0.51327, 0.49935, 0.48554, 0.47184, 0.45828, 0.44486, 0.43159, 0.41849, 0.40555, 0.3928, 0.38024, 0.36788, 0.35572, 0.34377, 0.33204, 0.32053, 0.30925, 0.2982, 0.28738, 0.2768, 0.26647, 0.25638, 0.24653, 0.23693, 0.22757, 0.21847, 0.20961, 0.201, 0.19264, 0.18452, 0.17665, 0.16901, 0.16162, 0.15447, 0.14755, 0.14086, 0.1344, 0.12816, 0.12215, 0.11635, 0.11077, 0.1054, 0.10023, 0.095263, 0.090491, 0.085911, 0.081517, 0.077305, 0.073269, 0.069406, 0.06571, 0.062177, 0.0588, 0.055576, 0.0525, 0.049566, 0.046771, 0.044108, 0.041574, 0.039164, 0.036873, 0.034697, 0.032631, 0.030671, 0.028813, 0.027052, 0.025385, 0.023807, 0.022315, 0.020905, 0.019573, 0.018316, 0.01713, 0.016012, 0.014958, 0.013966, 0.013033, 0.012155, 0.01133, 0.010555, 0.0098282, 0.0091459, 0.0085063, 0.0079071, 0.0073459, 0.0068208, 0.0063297, 0.0058707, 0.005442, 0.0050418, 0.0046684, 0.0043202, 0.0039958, 0.0036938, 0.0034126, 0.0031511, 0.002908, 0.0026822, 0.0024726, 0.002278, 0.0020976, 0.0019305
};
int SLOWER_frequency(SLOWER_FREQUENCY), SLOWER_period(0), SLOWER_amplitude(SLOWER_AMPLITUDE), SLOWER_index(0);
int SLOWER_size = sizeof(SLOWER_motion) / sizeof(SLOWER_motion[0]);
boolean SLOWER_state = false;

volatile boolean hz100_flag(false), hz20_flag(false);
long elapsed_time = 0;
int offset(OFFSET), scaling_factor(1);

/* FIR Variables */
double FIR_output(0), FIR_input(0);
float FIR_gain(1);
float FIR_coef[FILTERTAPS] = {
  0.065712, -0.017099, -0.016271, -0.0161, -0.016267, -0.016425, -0.016262, -0.015493, -0.013866, -0.011182, -0.0073126, -0.0022071, 0.0041003, 0.011493, 0.019777, 0.028693, 0.037929, 0.04714, 0.055962, 0.064028, 0.070981, 0.076483, 0.080236, 0.082223, 0.082223, 0.080236, 0.076483, 0.070981, 0.064028, 0.055962, 0.04714, 0.037929, 0.028693, 0.019777, 0.011493, 0.0041003, -0.0022071, -0.0073126, -0.011182, -0.013866, -0.015493, -0.016262, -0.016425, -0.016267, -0.0161, -0.016271, -0.017099, 0.065712
};

void setup() {
  Serial.begin(115200);

  /* Initialize Timers */
  startTimer(TC1, 0, TC3_IRQn, 100);
  startTimer(TC1, 2, TC5_IRQn, 20);

  /* Determining the period from BPM */
  FASTER_period = roundUP(60 / (double)FASTER_frequency * 1000, 10);
  SLOWER_period = roundUP(60 / (double)SLOWER_frequency * 1000, 10);

  analogWriteResolution(12);
  analogReadResolution(12);

  /* Initialize FIR */
  fir.setCoefficients(FIR_coef);
  fir.setGain(FIR_gain);
}

void loop() {
  /* Controlling indexing of waveforms */
  if (elapsed_time % (int)FASTER_period == 0) {
    FASTER_index = 0;
    FASTER_state = true;
  }
  if (elapsed_time % (int)SLOWER_period == 0) {
    SLOWER_index = 0;
    SLOWER_state = true;
  }
  if (FASTER_index >= FASTER_size - 1)
    FASTER_state = false;
  if (SLOWER_index >= SLOWER_size - 1)
    SLOWER_state = false;

  /* Timer @ 100 Hz - indexing through waveform arrays */
  if (hz100_flag) {
    hz100_flag = false;
    if (FASTER_state)
      FASTER_index++;
    if (SLOWER_state)
      SLOWER_index++;
    elapsed_time += 10;
  }

  FIR_input = scaling_factor * (FASTER_amplitude * FASTER_motion[FASTER_index] + SLOWER_amplitude * SLOWER_motion[SLOWER_index] + offset);

  /* Timer @ 20 Hz - calculating FIR */
  if (hz20_flag) {
    FIR_output = fir.process(FIR_input);
    serial_output();
    hz20_flag = false;
  }
}

Additional functions:

void serial_output() {
  Serial.print(FIR_input); Serial.print('\t');
  Serial.println(FIR_output);
}

int roundUP(int numToRound, int multiple) {
  if (multiple == 0)
    return numToRound;

  int remainder = numToRound % multiple;
  if (remainder == 0)
    return numToRound;

  return numToRound + multiple - remainder;
}

void TC3_Handler() {
  TC_GetStatus(TC1, 0);
  hz100_flag = true;
}

void TC5_Handler() {
  TC_GetStatus(TC1, 2);
  hz20_flag = true;
}

void startTimer(Tc *tc, uint32_t channel, IRQn_Type irq, uint32_t frequency) {
  pmc_set_writeprotect(false);
  pmc_enable_periph_clk((uint32_t)irq);
  TC_Configure(tc, channel, TC_CMR_WAVE | TC_CMR_WAVSEL_UP_RC | TC_CMR_TCCLKS_TIMER_CLOCK4);
  uint32_t rc = VARIANT_MCK / 128 / frequency; //128 because we selected TIMER_CLOCK4 above
  TC_SetRA(tc, channel, rc / 2); //50% high, 50% low
  TC_SetRC(tc, channel, rc);
  TC_Start(tc, channel);
  tc->TC_CHANNEL[channel].TC_IER = TC_IER_CPCS;
  tc->TC_CHANNEL[channel].TC_IDR = ~TC_IER_CPCS;
  NVIC_EnableIRQ(irq);
}

Here is the expected filtered output (Simulink) from the same input and this is what I am getting.

I am convinced am using this library incorrectly, I have some experience designing and testing discrete Fourier transforms in MATLAB then implementing them in Arduino - I have tested the FIR filter however this is the first time I am implementing it.

Any help would be much appreciated, thank you.

EDIT: Added entire working code.

2

Where is getdata()?
Are you sure your timer is not conflicting with motors PWM or inversely?
Have you done a small test code with just that part?

3

J-M-L:
Where is getdata()?
Are you sure your timer is not conflicting with motors PWM or inversely?
Have you done a small test code with just that part?

JML,

I am positive that nothing is conflicting with the other parts of the code. I understand how naive it can be if I assume that getData() can be taken for granted. I edited my original post with code that can compile - and I am getting the same output. Majority of this code is just to generate the data I want to filter.

4

I don't understand why everything in your loop() method is not conditional on

  /* Timer @ 20 Hz - calculating FIR */
  if (hz20_flag) {

I would expect to see

void loop() {
 /* Timer @ 20 Hz - calculating FIR */
  if (hz20_flag) {

  /* Controlling indexing of waveforms */
  if (elapsed_time % (int)FASTER_period == 0) {
    FASTER_index = 0;
    FASTER_state = true;
  }
  if (elapsed_time % (int)SLOWER_period == 0) {
    SLOWER_index = 0;
    SLOWER_state = true;
  }
  if (FASTER_index >= FASTER_size - 1)
    FASTER_state = false;
  if (SLOWER_index >= SLOWER_size - 1)
    SLOWER_state = false;

  /* Timer @ 100 Hz - indexing through waveform arrays */
  if (hz100_flag) {
    hz100_flag = false;
    if (FASTER_state)
      FASTER_index++;
    if (SLOWER_state)
      SLOWER_index++;
    elapsed_time += 10;
  }

  FIR_input = scaling_factor * (FASTER_amplitude * FASTER_motion[FASTER_index] + SLOWER_amplitude * SLOWER_motion[SLOWER_index] + offset);

 
    FIR_output = fir.process(FIR_input);
    serial_output();
    hz20_flag = false;
  }
}

of course removing the redundant conditionals!

5

stowite:
I don't understand why everything in your loop() method is not conditional on

  /* Timer @ 20 Hz - calculating FIR */

if (hz20_flag) {

You are right, for the most part. The if-statements relating to generating the input data (from the array) can be conditioned, however they would fall under the 100 Hz timer. My test waveforms (that input to the FIR) are sampled at 100 Hz. This setup it just simulating my real-world system.

In the real-world system, this data will be sampled at 20 Hz, which is why I designed the FIR at 20 Hz.

This may look ridiculous, but this is how I validate my discrete transfer functions (IIRs) in the past.

None-the-less, I dont think this will solve my problem.

6

The other potential problem I see is the loop taking more than 50 ms. You are doing a lot of serial IO and floating point operations (multiplications and additions) in that 50 ms. In your position I would would set a spare output high at the start of the calculation and low and the end and use an oscilloscope or logic analyser to time the calculation.

Out of interest, why is it necessary to use an FIR rather than an IIR ? IIR usually takes less computational resources than FIR for a similar filter characteristic.

7

stowite:
The other potential problem I see is the loop taking more than 50 ms. You are doing a lot of serial IO and floating point operations (multiplications and additions) in that 50 ms. In your position I would would set a spare output high at the start of the calculation and low and the end and use an oscilloscope or logic analyser to time the calculation.

Out of interest, why is it necessary to use an FIR rather than an IIR ? IIR usually takes less computational resources than FIR for a similar filter characteristic.

I will definitely check this out. I was worried about having so many taps - but I thought the Due could handle it without any worry.

I have been using IIRs for a while, they are fantastic. However I need the zero phase distortion that comes with FIRs. If the timing is indeed a problem, I think I can design my FIR with few taps.

Thanks

8

dgelman:
IHowever I need the zero phase distortion that comes with FIRs. If the timing is indeed a problem, I think I can design my FIR with few taps.

Presumably you only need zero phase distortion in the pass band. If so, why not use a Bessel IIR ? Of the 4 major filter types (Butterworth, Chebyshev, Bessel and elliptic) Bessel has the most linear phase response in the pass band.

I would be interested to know what post filter processing you are doing that requires zero phase distortion.

9

The filter is a 48-tap FIR low-pass filter sampled at 20 Hz

A designed sampling rate for the FIR of 20Hz means that the signal it is filtering must not have any frequency components above 10Hz. Even your slow frequency is above that.

My test waveforms (that input to the FIR) are sampled at 100 Hz.

Wouldn't they have to be sampled at 20Hz to match the filter's design rate?

I only know enough about FIR etc. to be dangerous, but for what you are doing I would have expected that the sampling rate would have to be at least 150Hz to be able to avoid aliasing with the "faster" signal. A sampling rate of about 200Hz would be appropriate along with a low-pass filter designed for 200Hz and which passes everything below, say, 20Hz.

I can't get your code to compile so I can't test it.

Pete

10

stowite:
Presumably you only need zero phase distortion in the pass band. If so, why not use a Bessel IIR ? Of the 4 major filter types (Butterworth, Chebyshev, Bessel and elliptic) Bessel has the most linear phase response in the pass band.

I would be interested to know what post filter processing you are doing that requires zero phase distortion.

The FIR filter will be used in a control loop. I am not using PID rather an error-based control system designed for processes with significant dead-time. The FIR will be used to isolate the SLOWER frequency in order to reject its disturbance. Unfortunately, the SLOWER frequency is not fixed and can slightly vary. The non-linear phase-delay will negatively impact the control system - it requires accurate measurements of delay. Using FIR, I see (in my simulated system) that I can reject the HIGHER frequency and measure the delay accurately.

I will take a look into Bessel IIR. The IIRs I have been using up to now are 2nd-order configurations (generally low-pass) - I discretize it, and implement the h[n] function into Arduino. I am sure it will be the same process for Bessel IIR.

11

el_supremo:
A designed sampling rate for the FIR of 20Hz means that the signal it is filtering must not have any frequency components above 10Hz. Even your slow frequency is above that.

Yes, the frequency components in my input are < 10 Hz. In this example, my SLOWER frequency is 0.3 Hz and my HIGHER frequency is 1.25 Hz.

el_supremo:
Wouldn't they have to be sampled at 20Hz to match the filter's design rate?

Yes, my real-world system, it is sampled at 20 Hz. The generated waveforms in the code are just used to test the FIR. In the real system, I am receiving the data from an external source at 20 Hz. The test waveform were made at a sampling frequency of 100 Hz, but they are being filtered at 20 Hz.

I can re-sample the test waveform for 20 Hz but it will not make a difference.

12

I think you ran off the end of your arrays. You check if the elapsed time modulo is zero, to reset the array indices back to zero if necessary, but then you go ahead and increment them before using them. I would always increment them, then check if they wrapped back to zero before using them.

Since this is still just a simulation, why not ignore real time altogether? Just increment elapsedTime once per loop and then index into the arrays as needed.

13

I am suspicious of

#include <FIR.h>
#define FILTERTAPS 48
FIR fir;

FILTERTAPS is also defined and used within FIR.h . You are trying to overwrite it. How do you know which FILTERTAPS is used when FIR source is compiled?

14

dgelman:
The FIR filter will be used in a control loop. I am not using PID rather an error-based control system designed for processes with significant dead-time. The FIR will be used to isolate the SLOWER frequency in order to reject its disturbance. Unfortunately, the SLOWER frequency is not fixed and can slightly vary. The non-linear phase-delay will negatively impact the control system - it requires accurate measurements of delay. Using FIR, I see (in my simulated system) that I can reject the HIGHER frequency and measure the delay accurately.

So all you need is a linear phase characteristic over the bandwidth represented by the range over which the SLOWER frequency can vary. How big is this range?

15

stowite:
So all you need is a linear phase characteristic over the bandwidth represented by the range over which the SLOWER frequency can vary. How big is this range?

Yes, I need a linear phase lag from 0.166 - 0.33 Hz, while filtering out the >1.25 Hz, the 48-tap FIR does this.

16

stowite:
I am suspicious of

#include <FIR.h>

#define FILTERTAPS 48
FIR fir;




FILTERTAPS is also defined and used within FIR.h . You are trying to overwrite it. How do you know which FILTERTAPS is used when FIR source is compiled?

Oh wow, I didn't see it begin define again in FIR.h.

I changed it to 48 and now I am getting something better - not perfect though.

17

That library cries out to be 'templatized' making FILTERTAPS a template parameter. I would not implement a default value; too dangerous.

18

dgelman:
Oh wow, I didn't see it begin define again in FIR.h.

I changed it to 48 and now I am getting something better - not perfect though.

The FIR library uses float. If on the Due float is 32 bits then you will be calculating FIR with 32 bits and I'm betting that all your PC based calculations use 64 bits floating point.

I use MATLAB and I'm pretty sure it uses 64 bit floating point.