# Linear system simulation with Python

Linear time-invariant (LTI) systems are widely used in the field of signal processing.  Scipy contains powerful tools for simulating LTI systems in the scipy.signal package, but they are not well documented.  I will provide a simple example that demonstrates how to use a few of the core classes and functions in scipy.signal for simulating LTI systems with Python.

## Define an LTI system

You will need to have Scipy installed, and you will need to have Matplotlib as well to make plots. We will start with an example of a first-order LTI system, which is characterized by a single parameter known as the time constant.  scipy.signal defines the class lti to represent a linear system in Python.

```#!/usr/bin/env python
from numpy import *
from matplotlib import pyplot as plt
from scipy import signal
tau = 5.0*60   # 5 minutes
h_times = arange(0.0, 10*tau, 0.1)
sys = signal.lti(1,[1,1.0/tau])
step_response = sys.step(T=h_times)[1]
plt.plot(h_times, step_response/step_response.max())    # normalized
plt.axhline(0.63, color='red')
plt.axvline(tau, color='red')
plt.xlabel('t')
plt.title('Step response')
plt.show()```

An lti object can be defined in three ways:

• If two arguments are provided, they are assumed to be the numerator and denominator of the transfer function of the system.
• If three arguments are provided, they represent poles, zeros, and gain.
• If four arguments are given, they must represent the state space (A,B,C,D) of the system.

In this case, we used the transfer function representation.  The transfer function is 1/(s+1/tau) where tau is the delay coefficient and s is the independent variable in the Laplace domain.  In Python, the denominator is represented as a list of coefficients, starting with the highest-order coefficient.  In our example, the coefficient of s to the first power is one, and the coefficient of s to the zero power is two.

## Step response

We can use the step_response method of an lti object to get the response of the system to a step function:
The blue curve shows the step response, while the red lines illustrate that the time constant is the amount of time required for the step response to reach a value of 1-1/e, or about 63.2% of its final value.
In a future post, I will cover some of the other things you can do with the tools in scipy.signal.

## References

### 8 thoughts on “Linear system simulation with Python”

1. I was wondering that too, it seems like it should be a time array, and the [1] is a reference to the documentation, noted in the “Refrence” heading below. Let’s see… if I say h_times = linspace(1,900,100) I get 100 points from zero to three time constants, which is accepted by sys.step(T=h_times), which returns both the time vector and the response.
The “T= ” portion of the call seems to be cuing the function as to which argument is being supplied. Coming from a c/c++ background I would have expected overloaded prototypes for each available combination of arguments, but this seems pretty slick.
http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.step.html#scipy.signal.step

2. Thank you for the comments,and I apologize for the mistakes in the original post. I added in a couple of lines of code that were missing and clarified some of the explanation.
@Andrew: you are correct that the construct “T=” in the call to sys.step(T=h_times) is a way of telling the function which argument is being supplied. This is a standard Python construct called a keyword argument.
http://docs.python.org/tutorial/controlflow.html#keyword-arguments

3. hi Craig, Thank you for the post . It was very helpful for me.
In your code step_response = sys.step(T=h_times)[1]
could you explain what is the significance of [1]

1. step is a function of the signals module in scipy.When you apply step to a system, this function will return 2 arrays, the time being the first and the corresponding signal being the second. Since we want the signal, we extract the 2nd array( represented by [1], 1st is represented by [0])

4. Thank you Craig. Very helpful post. Just a question, if i had to introduce a system with a time delay, how would i go about doing it?

This site uses Akismet to reduce spam. Learn how your comment data is processed.