4 Ways to Implement a Transfer Function in Code | Control Systems in Practice

In some situations, it is easier to design a controller or a filter using continuous, s-domain transfer functions. We have a lot of mathematical tools that make analyzing and manipulating them easier than, say, their equivalent time-domain differential equations. However, while we like to work with transfer functions for design and analysis, it’s often the case that we need to implement them in software that runs on some processor or microcontroller, and it might not be obvious how to go about doing this. So, in this video we’ll cover four simple ways to implement a transfer function in code. I’m Brian, and welcome to a MATLAB Tech Talk.

Before we get into any of the methods, we need to talk briefly about continuous-time and discrete-time systems. Transfer functions are continuous-time functions, so they produce continuously changing signals in the time domain. But software runs in discrete time. At each time step, the code runs and variables are updated. Those values are then held until the next step when the code runs again and they’re updated once more. So, the variables in a discrete system are only updated periodically.

Why is this a big deal? Well, when we implement a transfer function in code, we’re basically implementing a differential equation, which means we need to perform differentiation and integration. And since we only have periodic information in software, we have to figure out how to integrate the signal and differentiate the signal using these discrete points. There are many different ways to do this. To illustrate a few, let’s say we want to find the area under this continuous curve, and we’re sampling the signal at discrete time steps. The question is, how can we sum up the area under this curve if all we have access to are these data points?

Well, a simple way is to assume the value of the signal is held constant between time steps and so the area is just the sum of all of the signal value times the delta time. This is a first order method called forward Euler integration and it typically produces the most error as you can see by the fact that it grossly over or under estimates the area.

We can decrease the error by lowering the sample time, which means that we’re holding the value constant for a shorter time before the code runs again and updates the variable. But the lower the sample time, the less time we have to run the code since it has to finish running in this time. So, there is a limit to how short we can make this. In addition to lowering the sample time, we can also decrease error by choosing a different integration method. Other methods make different assumptions as to how the signal is changing between the discrete points. For example, second order trapezoidal integration assumes the signal linearly transitions from the previous time step to the current time step and then calculates the area of these polygons. And there are a number of higher order methods that use more of the past data points to try to create a more accurate representation of what the signal is doing between time steps. Typically, the higher the order the better the result.

Discrete differentiation is similar. A simple way to find the slope between two discrete points is to just subtract the two values and divide by the sample time. And just like with integration, we can improve this with higher order methods or by lowering the time between discrete points.

So, hopefully you can see how the accuracy of our code will depend on both the sample time, you know how fast the discrete points are updated, and the integration and differentiation methods that we choose. So, these are two things you need to consider when writing a continuous transfer function with discrete-time software. Now, with that being said, in the code I’m going to show you, I’ve just picked an arbitrary sample time and an arbitrary differentiation and integration method. When you write your own code, you can choose a different sample time and you can swap out the discrete methods for something else. It all depends on how accurate you want the result and how much processing power and memory the microcontroller has that you’re working with. So, let’s get to it!

Method 1: Numerically solve the differential equations

A transfer function is a differential equation that is represented in the s-domain rather than the time domain. And since our code is going to execute in the time domain, we will want to get back to the differential equations with the inverse Laplace transform. For example, we can multiply out the numerator and denominator and take the inverse Laplace transform to get this differential equation. The second derivative of y plus the first derivative of y plus 2 times y equals the first derivative of u + 2 times u. Where, u is the input signal, and y is the output signal. Taking the inverse Laplace transform like this is really easy because it can be done by inspection. Every s multiplied by the variable is just another derivative of that variable. So, s^2 times y(s) is the second derivative of y(t), and s time y(s) is the first derivative. In a differential equations class you may be asked to solve for y(t) analytically, but with software we can find an approximate solution, numerically. And that’s what we’re going to do.

To do this, we can rewrite the differential equation so that the highest derivative of the output is on the left side and the rest is on the right side. The basic idea for how we’re going to approach this code is this. In the first time step, we’ll solve for y double dot with our differential equation. Then in the next time step, we will integrate that y double dot to get y dot and then integrate y dot to get y. The input u is what is passed to the transfer function so we will always have a current u(t), and we can take the derivative of the input to get u dot. OK, now we have all of the information we need to solve for y double dot and then we can move onto the next time step and the process just repeats itself.

This works perfectly fine for every time step but, as you can see, it doesn’t work for the first one since at the beginning we don’t have y double dot to integrate which means we don’t have y or y dot. We also can’t determine how the input is changing from a single discrete input point, so we’re missing three-quarters of our equation. But that’s OK, we can just set these three to an initial value of our choosing and we’re all set.

So, there it is, a simple way to use discrete differentiation and integration to numerically solve a differential equation. Let me show you one version of this code in MATLAB. And if you are looking for C code, or something else that you can run embedded on a microcontroller, this will give you a good starting point to understanding the logic.

To get the true response, I’ll use MATLAB’s ability to handle continuous functions and plot the step response for our s-domain transfer function s + 2 over s^2 + s + 2. We’ll compare each of our discrete methods to this one to see how they do. In the code for Method 1, we build the time vector with the sample time of .1 seconds, then use that time vector to build the step input. We set the initial conditions, and predefine the size of the input/output vectors so that they aren’t changing size each iteration in the for loop.

Inside the for loop is where we code the differential equation. For the first step, we have to set the values to the initial conditions, like we said, and then we can solve for the second derivative of y, but for every other time step when can solve for u dot, y dot, and y, using a derivative, and two integrals. I’m using first order methods for differentiation and integration but again, this is where you can replace this with a higher order method of your choosing.

We can see that the result is OK. It follows the general shape of the true response, but the magnitude is a bit off. Let me change the sample time to 10 milliseconds just to show you the impact sample time has on the accuracy of the result. I’ll rerun the script and it’ll update the graph. There we go; this looks much better. So, numerically solving the differential equation like this is one way to code a transfer function.

For the second method we’re going to convert the transfer function to a state space model first

Setting up multiple derivatives and integrals in your code is straightforward but each time they run they add a little more error into the solution, so we’d do well to try to remove them where we can. It turns out, if we set up the problem differently, we don’t need the derivative at all. We can solve the problem with just a single integral that operates on the two-element state vector. To do this, we need to convert our transfer function into state space representation. There are simple ways to go from a transfer function to state space with hand calculations, but for us we’re just going to do that transformation with the ss() command in MATLAB.

This representation takes the high order differential equation and breaks it up into a number of first order differential equations. If the differential equations are linear (which they will be since transfer functions are linear) then we can package them into matrix form. I’ll come back to the values of these matrices in a second. Just realize that for state space representation, we only need four matrices: A, B, C, and D.

You can find a MATLAB Tech Talk video on state space in the description if you want more information. But for the purpose of generating code I’ll give you a fast introduction. We have two equations. The first is the state equation and it shows how the state vector changes over time, so how x dot varies based on the state vector, x, and the input vector, u. So, knowing x and u, we can solve for x dot. The states, however, aren’t necessarily the variables that we want to output so the output equation gives us the y. And y is a linear combination of the state vector and the input vector as well.

With these two equations, it’s pretty straightforward to write the code. To begin, we initialize the state vector to some initial condition. Then, in the first step, we use the initial states and the input to solve for x dot and the output y. Then in the next step, we integrate x dot to give us the new state vector for this time step and we use it and the input to again solve for x dot and y. Then we rerun this whole process over and over.

Let’s check out the MATLAB code. I’ve defined the A, B, C, and D matrices using the values that we got earlier. Then we set the initial states, and once again predefine the vector sizes. Now we can step through solving the set of first order differential equations. In the first step, we set everything to the initial conditions, but for every other step, we integrate x dot to get x, and then we use that and the current input, u, to get a new x dot and the output y. So, in this way, we only need to code a single integrator, which makes the code easier to read, and we are only integrating the two states in our state vector so that’s one less differential operation than we needed in method one. And as you can see, the result is very nearly the same as the other method, so our math is working. But if we zoom in, you’ll see that this method produces a result that is actually slightly closer to the truth than the first method and that is all because we’re not adding the additional error that comes from the extra discrete derivative operation.

OK, for the third method we’re going to convert the transfer function to the z-domain.

Both of the previous two methods converted the transfer function to continuous differential equations that were then discretized when we converted it to code. However, in this third method, we’re going to do the discretization first by converting the transfer function to the z-domain. There are a number of different methods for going from an s-domain transfer function to a z-domain transfer function and I’ve linked to a video series that goes over several of the more popular methods if you want to check that out. To convert our transfer function, we’re going to use the c2d function, or continuous to discrete function in MATLAB.

With c2d, we have to pass it the function we want to convert, of course. But we also have to select the sample time and the discretization method, which is effectively the integration method we want to use. This means that the sample time and integration method are already built in when we have a z transfer function. So, we don’t have to worry about writing our own integrator in code but we do have to make sure we run the resulting difference equation with the sample time that we used to create it. If we want to change the sample time or discretization method, we have to create a new z transfer function. This particular one is the result of a 10 millisecond sample time and the Tustin discretization method. The Tustin method is equivalent to saying that we’re using second order trapezoidal integration. Again, the linked videos below show why this is the case.

Now, we want to get this transfer function back into the time domain to write it in code, but luckily this is just as easy as it was with the inverse Laplace transform in the first method. The inverse Z-transform of 1/z is a one sample delay. So, if you have Y(z)/(U(z) = 1 + 1/z. We could take the inverse Z-transform to get that the output y is the current input plus the input from one step earlier. Since it’s this easy, we can essentially just do the conversion in our head and write out what the time-domain equivalent is directly. For our transfer function, we don’t have any obvious 1/z’s anywhere, but we can divide the numerator and denominator each by z^2 to write this as rational expressions of inverse z. Alternatively, in MATLAB we could change the displayed transfer function to use z inverse as the variable with this command.

To convert to a time domain difference equation, we can easily take the inverse Z-transform like this. We just replace each z inverse with a unit delay and solve for y. Or in other words, the output y is equal to .005 times the current input plus 9.95 e-5 times the input from one sample ago since there is a 1 over z, minus 0.005 times the input from two sample times ago since there are two 1 over z’s, and so on. Now, in code we can solve this difference equation in almost exactly the same way that we did in the first two methods when we solved continuous differential equation. We just choose the initial conditions and solve for y.

Let’s look at the code for this method.

We set up the time vector and step input, and then we initialize the variables. Since, our particular difference equation is using two past values for both the input and the out we have to initialize all four of these values. The for loop is still where we are solving the difference equations each time step, but we have to handle the initial conditions for the first two steps and then for each step afterward. We’re just multiplying and adding values together. You can see that the output is just a function of the last three inputs and the last two outputs. Again, we don’t have to explicitly code any integrators with this method since it’s already taken care of in the Z-transform process.

OK, check out this response, It looks really good. If I zoom in, you can see that it’s almost perfectly on top of the true response. Now, I don’t want this response to fool you into thinking this is the most accurate way to implement a transfer function in code. It’s producing a better response than the other two simply because the Tustin method that we use is equivalent to a second order integrator, which we know will produce a better result than the first order integrators that we used in the first two methods. If we upgraded the integrators in the previous two methods, they would have much better results as well.

Alright, for the fourth method we’re going to about autocode from MATLAB or Simulink

I’m not going to go into too much depth on this method since there are other good resources for generating code automatically using MATLAB and Simulink. I’ve linked to some in the description. I just want to point out that now that you know several different ways to code a transfer function by hand, you’ll have a better understanding of what auto coding is doing for us. Just like with writing code by hand, we still have to pick a sample time and we still have to pick an integration method. But with auto code the act of generating the C code is taken care of for you. This removes the possibility of making a mistake in the code. But also it allows you to more easily choose a higher order integration method with no extra effort on your part. You can just select from a number of built-in integration methods and you’re all set.

So, that’s where I’m going to leave this video. If you have a controller or a filter that you designed as a continuous function and you want to build C code or other code that is going to run on some other target hardware, hopefully you find one of these four methods helpful.

If you don’t want to miss future Tech Talk videos, don’t forget to subscribe to this channel. And you want to check out my channel, control system lectures. I cover more control theory topics there as well. Thanks for watching, and I’ll see you next time.