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

- Scipy.signal documentation
- Feedback Systems online text by Astrom and Murray
- Linear Time Invariant (LTI) system theory
- Transfer function
- Convolution
- Time constant
- Laplace transform

what about h_times ?

LikeLike

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

LikeLike

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

LikeLike

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]

LikeLike

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

LikeLike

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

LikeLike

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?

LikeLike

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

LikeLike