# Implementing a continuous time transfer function in PC code

I have a transfer function like this:

$H(s) = \frac{1}{Ts + 1}$

I converted it into a difference equation to solve it iteratively:

$H(s) = \frac{Y(s)}{X(s)} = \frac{1}{Ts + 1} \\ TsY(s) + Y(s) = X(s) \\ \mbox{Assuming the initial conditions are zero,} \\ T\frac{dy(t)}{dt} + y(t) = x(t) \\ y(t) = x(t) - Ty'(t)$

I wrote a code like this;

const double dt = 0.001;
double x;       // Input
double y;       // Output
double yp=0;    // Previous output value
double yd;      // Derivative of output
for(double t=0; t<TIMEMAX; t+=dt)
{
// ...
yd = (y - yp) / dt;
yp = y;
y = x - T*yd;       // y(t) = x(t) - Ty'(t)
StoreNewOutput(y);
//...
}


T symbolizes sensor delay (the sensor which reads x(t)) in my code. If T=0 (no sensor delay) my code runs very good. But if I set even a very small sensor delay (e.g.; T=0.1s) output becomes unstable (y approaches to infinity after a few iterations).

Am I doing something wrong? This is my first time implementing a transfer function in a computer algorithm like this. Are realization of continuous time transfer functions done like this, or do people use a different method? Please, can you confirm validity of my method, or make corrections on it?

• As soon as you talk about sampling, you have a discrete time system rather than a continuous time one. As you are discovering, the rules are a bit different. – Chris Stratton Dec 24 '12 at 4:00

You are asking about the world of numerical solutions for ordinary differential equations. And I think what you are trying to get to in your code is essentially Euler's Method for solving this kind of equation. However, I think you've somewhat turned yourself around in analyzing the problem.

Instead of your equation, $y(t) = x(t) - Ty'(t)$, it's more common to analyze this problem by rearranging terms:

$y'(t) = \frac{1}{T}(x(t)-y(t))$

To apply Euler's method, we estimate the evolution of $y(t)$ by

$y(t+h) = y(t) + hy'(t)$,

where h is a timestep for the numerical solution (the dt in your code). We can put this in discrete time notation as:

$y_{n+1} = y_n + hy'_n$

This is a general formulation for any ordinary differential equation. For your problem, you'd have

$y_{n+1} = y_n + h\frac{1}{T}(x_n - y_n)$.

Which you can easily translate into C code.

That said, this method tends to not be particularly accurate. The overall error in the solution grows in proportion to the step size of the simulation, and the numerical method itself can be unstable (errors can go to infinity even if the real solution is finite) under some circumstances.

The usual method (a method that works for many problems; something to try first before looking for more advanced methods in case your problem turns out to be especially difficult) is fourth-order Runge-Kutta integration. The equations for this are somewhat more involved, but 1) many pre-written libraries are available to do it, and 2) if you work through the equations carefully you'll find they all fall into place very nicely. There are also numerous elaborations available in libraries, for example to automate finding the largest step size (and thus least processing time) to obtain the required accuracy in the solution.,

• Beat me to it by about 30 seconds! – Alfred Centauri Dec 24 '12 at 3:00
• @AlfredCentauri, It wouldn't have taken much to beat me. It took me quite a while to dredge my brain (and Wikipedia) to remember how to do this stuff. I was expecting "1 new answer to this question" to appear the whole time I was writing it. – The Photon Dec 24 '12 at 3:13
• @ThePhoton Thank you for your answer. My code runs as expected now. In my code, T is defined as a constant, I am changing it to different values and running the code to see the result. To use your method, the step size must be much smaller than T (i.e.; h<<T), isn't it? Also, I can't set T=0 (which would have yielded $y_{n+1} = y_n$ automatically), because h has a non-zero practical value. – hkBattousai Dec 24 '12 at 11:53
• @hkBattousai, setting T = 0 gives H(s) = 1, i.e., y = x. Also, the Euler method relies on the higher order terms in the Taylor expansion being insignificant. In your case, the 2nd order term is $\frac{\dot x - \dot y}{2T} (\Delta t)^2$. So, the time increment must be small enough that this term is much smaller than the linear term. Note that there is a dependence on the rate of change of x(t) so, if x changes slowly enough, the time step can be relatively large. – Alfred Centauri Dec 24 '12 at 13:42

Your filter would be more stable calculated as an IIR filter using two coefficients.

The transfer function

H(s)=1/(T * s + 1)

Is a first order low pass filter having cutoff frequency 1/T.

To calculate the result y[n] using discrete time step dt, at each time n * dt, as a discrete time IIR filter one general form is...

y[n] = Beta * y[n - 1] + Alpha * x[n]

y[n] is the filter output at time n*dt
y[n-1] is the value of the filter at time (n - 1) * dt, the previous time step.
x[n] is the filter input at time n*dt
Beta = exp(-dt/T)
Alpha = dt

Note that this form is unconditionally stable as long as Beta is less than 1, which it will be for all positive time steps and cutoff frequencies. When computed on a digital device with finite precision there is the additional constraint that your floating point numbers may need to round towards 0 for stability.

The derivation is as follows...

H(s)=1/(T * s + 1) = 1/T * 1/(s + 1 / T)

Taking the inverse laplace transform to find the time domain impulse response we get...
h(t) = 1/T * e^(-t/T)

In general...
To calculate the response of the transfer function h(t) to an input signal x(t) in continuous time take the convolution...

y(t) = integral(h(t - tau) * x(tau) * dtau), from tau = -infinity to t

in the discrete time case, using time step dt, the integral is approximated as a sum, and the convolution is...

y[n] = sum(h[n-i] * x[i] * dt), summed from i = -infinity to n, the current simulation time is n * dt

y[n] is the filter output at time n*dt.
h[n-i] is the impulse response of the filter at time (n-i)dt
x[i] is the value of the input signal at time i
dt

Given your specific transfer function, the discrete time convolution is...

y[n] = 1/T * sum(e^-((n-i)*dt/T) * x(i*dt) * dt), summed from i = -infinity to n
Note that the following recursion relationship exists...
y[n] = y[n-1] * e^(-dt/T) + 1/T * x[n] * dt

To generalize ever further, nearly any polynomial transfer function of the form...
Nu[n] * s^n + Nu[n-1] * s^(n-1) + ... + Nu[1] * s + Nu[0] / ( De[n] * s^n + De[n-1] * s^(n-1) + ... + De[1] * s + De[0] )

Where Nu[] are constant Numerator coefficients.
De[] are constant Denominator coefficients.

Can be factored into the form...

(Nu[n] * s^n + ... + Nu[0])) * ( C[n]/ (s-R[n]) + ... + C[0] / (s-R[0]))

Where C[n] are constants, and R[n] are the roots of the denominator in equation.

From the form of the equation its clear that in the time domain this is just a series of exponentials and their derivatives. For an equation having X exponential terms you get...

y[0][n] = Beta[0] * y[0][n-1] + Alpha[0] * x[n]
y[1][n] = Beta[1] * y[1][n-1] + Alpha[1] * x[n]
...
y[X][n] = Beta[X] * y[1][n-1] + Alpha[X] * x[n]

Where y[i][n] is the partial filter response y[i] at time n * dt due to one of the exponential terms.
Beta[i] is the constant decay factor for partial filter response y[i].
Alpha[i] is the constant input weighting factor for partial filter response y[i]
x[n] is the filter input at time n * dt.

The total result is the sum of the result from each partial filter response.

y[n] = y[0] + y[1] + ... + y[X]