How to implement a discrete transfer function into an arduino?

Question

I am here asking how does one transfer a difference equation into a MCU? I have never done it personally and looking into this topic I was never able to find a good answer. Its usually used with Matlab but never a MCU.

I have a continuous transfer function:

$$LPF_A = \frac{1}{6.416\cdot 10^{-11}s^2+1.133\cdot 10^{-5}s+1}$$

Its a 2nd Order Butterworth filter with a \$F_c = 20\text{kHz}\$.

Using the Tustin method I have the discrete equation:

$$LPF_D=\frac{0.1441z^2+0.2881z+0.1441}{z^2-0.6777z+0.254}$$ With a sample of \$T_s = \$ 8e-6 seconds, in scientific E notation.

The difference equation:

\$Y\$ = Output

\$U\$ = Input

$$Y_i =0.1441U_{i}+0.2281U_{i-1}+0.1441U_{i-2}+0.6777Y_{i-1}-0.254Y_{i-2}$$

I know I have to make my ADC output = \$U_i\$ and initialize \$U_{i-1}\$, \$U_{i-2}\$, \$Y_{i-1}\$, \$Y_{i-2} = 0\$, But after that no idea.

I have no idea what to do with this after. I have always used Matlab's embedded ways of doing it but never have transferred a difference equation into a MCU before. Any help would be appreciated.

Edit: Here's the code based on what I think

  unsigned short x;
double Y_i[10] = {0,0,0,0,0,0,0,0,0,0};
double U_i[10] = {0,0,0,0,0,0,0,0,0,0};
double Y_i_1 = 0;
double Y_i_2 = 0;
double U_i_1 = 0;
double U_i_2 = 0;

void setup() {

  ADMUX |= 0xC0;
  ADCSRA |= 0x8F;
  ADCSRA |= 0x40;
  Serial.begin(9600);
}

ISR(ADC_vect) {
 x = ADCL << 2;
 x |= ADCH;
 ADCSRA |= 0x40;
}

void loop() {
  //Yi = 0.1441Ui + 0.2281Ui-1 + 0.1441Ui-2 + 0.6777Yi-1 - 0.254Yi-2 
  for (int i = 1; i < 4; i++){

    if (i == 1) {
       Y_i[i] = 0.144*x + 0.2281*U_i_1 + 0.1441*U_i_2 + 0.6777*Y_i_1 - 0.254*Y_i_2;
       U_i[i] = x;
    } else if (i == 2) {
      Y_i[i] = 0.144*x + 0.2281*U_i[i-1] + 0.1441*U_i[i-2] + 0.6777*Y_i[i-1] - 0.254*Y_i[i-2];
      U_i[i] = x;
    } else {
      Y_i[i] = 0.144*x + 0.2281*U_i[i-1] + 0.1441*U_i[i-2] + 0.6777*Y_i[i-1] - 0.254*Y_i[i-2];
      U_i[i] = x;
    }

  }

}

What happens after the iteration is done? It resets? thats it?

Answer 1

Maybe you are having trouble with the "advance time" step. It helps to separate the logic as follows:

1. Declare your variables: \$Y_{i}\$, \$Y_{i-1}\$, \$Y_{i-2}\$, \$U_{i}\$, \$U_{i-1}\$ and \$U_{i-2}\$. Let's also initialize a counter \$i = 0\$, which represents the current value of the subscript in your variables (this counter is completely unnecessary, but helps getting across what we are doing). Here is how \$i\$ works: if \$i = 10\$, this means the variable for \$Y_{i-2}\$ currently holds \$Y_{8}\$. Initialize everyone to zero, but keep in mind \$Y_{i}\$ and \$U_{i}\$ are not valid yet - we'll deal with them in the next two steps.
2. Acquire new data. That is, acquire \$U_{i}\$ from the ADC.
3. Recalculate your output. You know how to do it. Store it in \$Y_i\$.
4. Profit. Now all your variables are correct and reflect the system at the timestep \$i\$ (if this is your first pass, now \$i = 0\$). Do whatever you want with them.
5. Now you have to advance time. This means updating your variables as if time passed. This is necessary for calculating the next timestep. To do this, have \$Y_{i}\$ become \$Y_{i-1}\$, and so on (I won't spoil everything you have to do). You'll notice old values of \$U\$ and \$Y\$ will be lost as time advances. Also increment \$i\$ to keep track of which sample we are talking about.
6. Now, all your variables except \$Y_i\$ and \$U_i\$ reflect the state of the sistem at the current \$i\$ (if this is your first pass, now \$i = 1\$). \$Y_i\$ and \$U_i\$ are still unknown because you need another sample to calculate them.
7. Wait for the time to sample again and go to step 2.

Answer 2

Well, every iteration update your variables like this:

$$ U_{i-2} \leftarrow U_{i-1}$$ $$ U_{i-1} \leftarrow U_{i}$$ $$ U_{i} \leftarrow \text{curent ADC reading}$$ $$Y_{i-2} \leftarrow Y_{i-1}$$ $$Y_{i-1} \leftarrow Y_i$$ $$Y_i \leftarrow \text{calculate current output}$$

In you code you are using double variables for precision, but, for the Arduino Uno and other that use AVR architecture microcontrollers double is the same as float. If you intend on keep using double, notice you will not have the desired precision.

Also, if you are using the Arduino Uno you will not be able to do many floating point operations per cycle, someone got a figure of 9 us per operation (except division), so just the part where you calculate your output \$Y_i\$ would use 5 multiplications and 5 additions, taking 90 us. Your filter is designed for a sampling time Ts= 8e-6 seconds, so you would not be able to finish the first floating point multiplication before getting your second reading from the ADC, and would have done 10 more readings before calculating the first output \$Y_i\$, long story short, the microcontroller will not be able to keep up and run that filter at that speed.

So you are way over your budget in frequency, if you where doing only that calculation, you would be barely able to achieve a 10 kHz frequency.

Finally, you could improve your code in two ways, first

//you do not need all those values, 
//from your expression for the filter, you could use just this
//double Y_i[10] = {0,0,0,0,0,0,0,0,0,0};
//double U_i[10] = {0,0,0,0,0,0,0,0,0,0};
double Y_i[3] = {0,0,0};
double U_i[3] = {0,0,0};



//not needed, just load the values you want as initialization in the array above 
//double Y_i_1 = 0;
//double Y_i_2 = 0;
//double U_i_1 = 0;
//double U_i_2 = 0;

notice that

• \$ U_{i-2}\$ = U_i[2]
• \$ U_{i-1}\$ = U_i[1]
• \$ U_{i}\$ = U_i[0]
• \$Y_{i-2}\$ = Y_i[2]
• \$Y_{i-1}\$ = Y_i[1]
• \$Y_{i}\$ = Y_i[0]

And that filter expression with a bunch of if could be just

void loop() {

    U_i[2] = U_i[1];
    U_i[1] = U_i[0];
    U_i[i] = x;
    Y_i[2] = Y_i[1];
    Y_i[1] = Y_i[0];
    Y_i[0] = 0.144*U_i[0] + 0.2281*U_i[1] + 0.1441*U_i[2] + 0.6777*Y_i[1] - 0.254*Y_i[2];
    //send or store any of the values that you need, 
    //Y_i[0] will always have the latest update of the filter, and so on...


    //either put this filter in some non-critical interruption to execute it with the desired periodicity (best way) 
    //or use some precise delay here right after calculating the filter inside the loop
    delayMicroseconds(50);

}