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.

Microsoft Word is not the best tool for doing technical writing.  However, sometimes we are required to use Word because we need to collaborate with others who want to use Word.  Right now, I am using Word 2007  on Windows XP to write several mathematical papers.  In general, it is a big improvement from previous versions.  The new integrated equation editor is outstanding–except for the major bug I’ll discuss below.  Here is a brief “FAQ” you will want to read if you are using Word for technical writing. **Q:**How do I enter multi-line equations in Word 2007?**A:Press Shift-Enter where you want a line break to appearQ:**How do I prevent a page break from splitting a table into two parts?**A:**Select the table.  On the Home tab, click the little box in the lower-right corner of the Paragraph box.  On the paragraph dialog, choose the “Lines and Page Breaks” tab.  Check the box for “Keep Lines Together.”**Q:**How do I automatically number headings in Word 2007?  For example: Section 1, Subsection 1.1, 1.2, 1.3, Section 2, etc.**A:**It’s not obvious.  See this page from dummies.com about number headings.**Q:**In Word 2007, why do equations sometimes appear as blank spaces or question marks when I print or save my document as a PDF file?**A.**This occurs when Word is installed on Windows XP Pro.  See the following Microsoft tech support item to find out how to install missing scripts: The characters in an equation are not printed… You may also have an outdated printer driver: Microsoft Support Item 920228**Q:**Why do equations created in Word 2007 disappear when I open the document in Word 2008 for Mac?A. Word for Mac does not support equations written in the new Word 2007 equation editor. Unfortunately, neither does PowerPoint 2007 on the PC.  You can work around this by inserting equations into Word the old-fashioned way: go to the Insert tab, click on Object (found in the “Text” box towards the right side of the tab), and choose “Microsoft Equation 3.0” from the list in the dialog box.

Scipy contains a good least-squares fitting routine, leastsq(), which implements a modified Levenberg-Marquardt algorithm.  I just learned that it also has a constrained least-squared routine called fmin_slsqp().   I am using simple upper and lower bound constraints, but it’s also possible to specify more complex functional constraints. What I did not realize, at first, is that fmin_slsqprequires a different type of objective function thanleastsq.   leastsqrequires you to write a function that returns a vector of residuals, andleastsqautomatically squares and sums the residuals.fmin_slsqp is actually more flexible, in that it can use any objective function that returns a single scalar value.  To implement least-squares curve fitting, your objective function will need to find the residual at each data point, square the values, and sum them up.  Hopefully this tip will save you some time. Check out the scipy optimization tutorial for more examples.  Here is the original paper by Dieter Kraft which introduces the algorithm used by fmin_slsqp.

In a previous post, I described how Python’s Pickle module is fast and convenient for storing all sorts of data on disk. More recently, I showed how to profile the memory usage of Python code.  In recent weeks, I’ve uncovered a serious limitation in the Pickle module when storing large amounts of data: Pickle requires a large amount of memory to save a data structure to disk. Fortunately, there is an open standard called HDF, which defines a binary file format that is designed to efficiently store large scientific data sets. I will demonstrate both approaches, and profile them to see how much memory is required. I am writing the HDF file using the PyTables interface. Here’s the little test program I’ve been using: [code language=“Python”]#!/usr/bin/env python from numpy import array

I’m using GNU Fortran (gfortran) to build several shared libraries, and then dynamically linking to them from a Fortran program.  The process is a little different than what I’m used to for C libraries, so I thought I’d explain it.  Unlike C, there is no need to #include header files when compiling code that relies on functions defined in an external library.  Likewise, there is no need to use -l or -L linking flags to tell the linker about s hared libraries (at least when they’re in the same directory).  In fact, the whole process requires a lot less command-line options than I had expected.

In a previous post, I explained how to use the Python profiler.  The profile is great for finding out which parts of the code run the slowest, or are called most often.  However, the profiler doesn’t give any information about how much RAM is being consumed, or where it’s being consumed.  If your program needs so much memory that it starts swapping to disk, its speed can be reduced by orders of magnitude.  On the positive side, your code may run much faster if it fits entirely in the processor cache.  In this post, I will introduce two tools that can help you understand the RAM usage of your Python code.

Lookup tables and spline fitting are widely used by scientific programmers.  A particular function may not have an analytic solution–in other words, it can’t be expressed as an equation of elementary functions.  This might happen if the function were empirically determined from experimental data, or if the equation can’t be algebraically solved for one variable.  If the analytical function is available, but takes a long time to evaluate, a lookup table or spline approximation can be considerably faster.  In a previous post, I showed how to use the function interp1d from scipy.interpolate as a lookup table.  In a later post, I showed that interp1d is actually rather slow, and scipy.interpolate.UnivariateSpline is much faster.  Now, I will show some benchmark results, and explain a potential pitfall when using UnivariateSpline. The two functions I will use for this demonstration are not very complicated.  Each function is a piecewise approximation of a more complex function.  Here is one of the functions.  For small z, the function uses one approximation, and for large z, it uses a far-field approximation.  For intermediate values, a polynomial is used to smoothly match the two solutions.

There are lots of good open-source tools that you can use to make high-resolution, publication-quality 2D plots.  Personally, I like to use Python, numpy, and matplotlib.  Unfortunately, it is much harder to find a good tools to make 3D plots.  Older versions of matplotlib had rudimentary 3D support, but this was removed in version 0.98.  In this post, I will review a Python 3D plotting library called PyX.

Preparing the data

Figuring out how to store the data to be plotted was actually the hardest part of learning to use PyX.  The data format for 3D plots is not well documented.  PyX requires a list of (x,y,z) lists like this:

I have released a new version of my Python class that generates VTK data files in the .vtu format, which is compatible with Paraview and other VTK applications.  If you have downloaded the old one, please get the latest version, which incorporates some bug fixes and has been more thoroughly tested.