3
\$\begingroup\$

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?

\$\endgroup\$
3
  • \$\begingroup\$ What is "e" in your transfer function? \$\endgroup\$
    – Andy aka
    Mar 31 '20 at 18:34
  • \$\begingroup\$ Sorry its, exponent to. 5e5 = 5*10^5 \$\endgroup\$
    – Leoc
    Mar 31 '20 at 18:38
  • \$\begingroup\$ I hate e notation because the natural logarithm exists with the same notation. \$\endgroup\$
    – DKNguyen
    Mar 31 '20 at 19:22
1
\$\begingroup\$

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.
\$\endgroup\$
14
  • \$\begingroup\$ Let me know if I am wrong here but: Initialize all the variables to 0 Set Ui = ADC First Iteration: Yi = a Constant. Second Iteration: Did the Yi before becomes Yi-1 now? After I I iterate to Yi-2 what now? Is it done? can I not sample "Continuously" now? \$\endgroup\$
    – Leoc
    Mar 31 '20 at 19:41
  • \$\begingroup\$ "Did the Yi before becomes Yi-1 now?" Yes. I did not understand your other question \$\endgroup\$
    – FrancoVS
    Mar 31 '20 at 21:04
  • \$\begingroup\$ Sorry, the other question was stating what do you have you solve the second iteration. To me it looks like you need to solve the two iterations and thats it, what happens after that? \$\endgroup\$
    – Leoc
    Mar 31 '20 at 21:45
  • \$\begingroup\$ When solving for the two iterations do I need to put this in a loop or can I first solve them then once I do put the equation in a loop? \$\endgroup\$
    – Leoc
    Mar 31 '20 at 21:48
  • \$\begingroup\$ every iteration is the same: you acquire Ui, calculate a new Yi, and advance time. \$\endgroup\$
    – FrancoVS
    Mar 31 '20 at 21:48
1
\$\begingroup\$

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);

}


\$\endgroup\$
6
  • \$\begingroup\$ Why did you make this a community wiki? \$\endgroup\$
    – pipe
    Mar 31 '20 at 19:01
  • \$\begingroup\$ So people can make it more complete with information regarding how to deal with the real life problems of implementing the algorithm, specially regarding the very limited FPU capability of the Arduino the OP is going to use and the errors due to quantization. \$\endgroup\$
    – jDAQ
    Mar 31 '20 at 19:47
  • \$\begingroup\$ For people downvoting, either point out the problem or feel free to add/correct things in the answer... \$\endgroup\$
    – jDAQ
    Mar 31 '20 at 20:51
  • \$\begingroup\$ I didn't see this edit and it was extremely helpful thank you. I have a few questions. When you say I am over budget for frequency, I understand why based on what you said, then how would one find out if a microcontroller is capable of needed frequency then? The way I did it was Arduino nano has a clock speed of 16MHz in the ADC register prescale it with 128, so 16MHz/128 = 125KHz, I would imagine the ADC is now sampling at Ts = 1/125KHz. If you are saying the loop() function is operating at a different speed then would using a timer ISR and shoving those equations in it would fix it? \$\endgroup\$
    – Leoc
    Apr 1 '20 at 17:59
  • \$\begingroup\$ "When you say I am over budget for frequency, I understand why based on what you said, then how would one find out if a microcontroller is capable of needed frequency then?' this can be a question itself, some similar ones are electronics.stackexchange.com/questions/151550/… electronics.stackexchange.com/questions/137742/… \$\endgroup\$
    – jDAQ
    Apr 1 '20 at 23:02

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.