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:

LTI System Step Response

LTI System Step Response

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

  1. Scipy.signal documentation
  2. Feedback Systems online text by Astrom and Murray
  3. Linear Time Invariant (LTI) system theory
  4. Transfer function
  5. Convolution
  6. Time constant
  7. Laplace transform

Comments

Comment #1 by Gabriel

Gabriel - Jul 5, 2011

what about h_times ?

Comment #2 by Andrew

Andrew - Sep 4, 2011

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

Comment #3 by craig

craig - Feb 3, 2012

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

Comment #4 by Sivaram

Sivaram - Sep 6, 2012

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]

Comment #5 by craig

craig - Sep 1, 2012

Sivaram, The method sys.step(T) returns the tuple: (T, yout) The [1] selects the second element of the tuple, yout, which is the step response of the LTI system that we have defined. Here is more about the step method: http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.step.html#scipy.signal.step …and here is more about tuples: http://docs.python.org/tutorial/datastructures.html#tuples-and-sequences

Comment #6 by MoMo

MoMo - Mar 6, 2013

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?

Comment #7 by Igneel Dragneel

Igneel Dragneel - Nov 0, 2014

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

Comment #8 by Robert Papanicola

Robert Papanicola - Aug 2, 2016

Hello is it possible to take into account non-linear elements (eg saturation) with scipy As in the example: http://stackoverflow.com/questions/38999317/how-to-simulate-saturations-and-thresholds-with-scipy