Class: Psych 204B
Tutorial: Linear Systems
Author: Wandell
Date: 3.31.04
Duration: 30 minutes
Copyright: Stanford University, Brian Wandell
Checked: 10/13/09 JW
Checked: 10/13/10 BW
Translated to Python by Michael Waskom, October 2012
In [ ]:
# First set up some plotting things as usual
%pylab inline
import matplotlib as mpl
mpl.rcParams["figure.figsize"] = (9, 6)
mpl.rcParams["axes.grid"] = True
mpl.rcParams["lines.linewidth"] = 1.5
This tutorial introduces the student to the basic methods of linear systems theory, focusing on the tools of shift-invariant linear systems. The tutorial introduces the idea of superposition, shift-invariance, and impulse response functions. These ideas are fundamental to many methods and analyses used in functional neuroimaging and through engineering and science.
Linear systems theory is a method of characterizing certain types of common systems. A system is something that has an input and an output, and thus we can think of it as a function output = L(input).
For example, the system might be an optical element, like a lens, that takes an image as an input and produces another image as an output. Or, this system may be a biological system that we are measuring, for instance the BOLD response that is the consequence of neural activity.
Characterizing the complete input-output properties of a system by exhaustive measurement is usually impossible. When a system qualifies as a linear system, it is possible to use the responses to a small set of inputs to predict the response to any possible input. This can save the scientist enormous amounts of work, and makes it possible to characterize the system completely.
Not all systems can be described using linear systems theory. To use linear methods to describe a system, it must satisfy the basic principle of superposition.
The principle of superposition in words is simply this:
If the response to input A is Ra, and the response to input B is Rb, then the response to the sum of A and B (A+B) is the sum of the responses: Ra + Rb.
In equation form this is
There is one simple consequence of superposition that is often taken as
a special case and named the principle of homogeneity.
Suppose that the inputs are equal, namely, A = B.
It is easy to see that L(A + A) = Ra + Ra and thus L(2A) = 2Ra.
In general, if s is any scalar, a linear system will satisfy the rule:
An enormous number of simple systems satisfy superposition over a fairly large range of inputs. Linear systems theory is a basic tool used throughout engineering, physics, and statistical modeling.
In this tutorial, we will build up some of the tools of linear systems theory. These include the use of matrices to describe a linear system, the very important special case of a shift-invariant linear system, and the associated ideas of the impulse response function
All of these ideas are used in functional neuroimaging. In fact, nearly every field that uses linear systems theory develops its own terminology (jargon) to describe the general ideas of linear systems theory in the specific context of that field.
In neuroimaging, the hemodynamic response function (HRF) is the impulse response function. The analysis of the modulations in block-design experiments are closely connected to the system transfer function.
time varying stimuli can be thought of as a series of impulses at different levels. We can approximate any complex stimulus as if it were simply the sum of a number of shifted and scaled impulses.
A shift-invariant linear system is a special case when the system responds in the same way to every impulse, no matter when the impulse occurs in time. That is, it doesn't matter what happened just before the impulse or what happens just after the impulse. The response to the impulse will be the same.
The only difference in the response to different impulses is that the responses are shifted in time to correspond to the impulse and multiplied by a scalar that equals the size of the impulse.
The significance of the principle of superposition is this:
If we know the impulse response function, we can predict the response to any stimulus (combinations of impulses). We compute the response to each impulse and then add up all the results.
Practically, to characterize shift-invariant linear systems, then we need to measure only one thing: the way the system responds to an impulse of a particular intensity. This response is called the impulse response function of the system.
The analysis of BOLD data presupposes that the responses are shift-invariant linear systems. In this field, the impulse response function is called the hemodynamic response function, or HRF.
In neuroimaging, the principles of superposition and shift invariance are essential to the ideas used in rapid event-related designs: The underlying assumption is that if you have two neural (or cognitive) events A and B occurring one after the other, then the BOLD response of A and the BOLD response of B will be identical (but B will be shifted in time with respect to A) and thus total BOLD signal after these two events is a sum of two independent responses shifted in time:
BOLD(A+B) = BOLD(A) + BOLD(B).
Now, we do some simple calculations. But first, we define a function to plot a model hemodynamic response function.
In [ ]:
from scipy.special import gammaln
def spm_hrf(RT, fMRI_T=16):
"""Python implementation of spm_hrf"""
_spm_Gpdf = lambda x, h, l: exp(h * log(l) + (h - 1) * log(x) - (l * x) - gammaln(h))
dt = RT / float(fMRI_T)
u = np.arange(1, int(32 / dt) + 1)
hrf = np.concatenate(([0.0], _spm_Gpdf(u, 6 , dt) - _spm_Gpdf(u, 16, dt) / 6))
idx = arange(int((32 / RT)) + 1) * fMRI_T
hrf = hrf[idx]
hrf = hrf / sum(hrf)
return hrf
Let's take a look at this model (or "canonical" function)
In [ ]:
RT = 1 # Repetition time in seconds (sample every second)
hrf = spm_hrf(RT)
In [ ]:
nT = len(hrf) # Number of time steps
t = arange(nT) # Individual time samples of the HRF
plot(t, hrf)
xlabel('Time (sec)')
ylabel('HRF level');
Notice that the hrf values sum to 1.
This means that convolution with the HRF will preserve the mean of the input stimulus level.
In [ ]:
print(sum(hrf))
Also, notice that the default SPM HRF has a post-stimulus undershoot.
Next we will plot several instances of impulse functions with the corresponding estimated hemodynamic response using this model. To create the estimated response trace, we will use the convolve function from the numpy package (loaded automatically as we are using IPython in pylab mode; see above).
According to the principle of superposition the output should be the sum of the individual events. Thus, the output can be written as the sum of the outputs of stimuli 1 and 2.
Convolution is a mathematical function that applies the superposition principle: basically it sums over all impulse responses shifted in time
output = sum(hrf * input[t])
Since we'll be doing this a lot, let's write a function that takes a list of impulse times and makes the two plots.
In [ ]:
def plot_imp_hrf(times):
"""Make a plot of impulse functions and estimated hemodynamic response."""
# Get the predicted model
hrf = spm_hrf(1)
nT = len(hrf)
t = arange(nT)
# Make impulse vectors for each input stim
stims = []
for time in times:
stim = zeros(nT)
stim[time] = 1
stims.append(stim)
# Plot the stimulus
figure(figsize=(9, 4))
subplot(1, 2, 1)
plot(zeros(nT))
for time in times:
plot([time, time], [0, .5])
axis([0, 30, -0.5, 1])
xlabel("Time (s)")
title("Stimulus")
# Convolve the stimulus with the hrf model
sum_stim = sum(stims, axis=0)
bold = convolve(sum_stim, hrf)
# Plot the expected BOLD response
subplot(1, 2, 2)
plot(bold)
axis([0, 30, -0.5, 1])
xlabel("Time (s)")
title("BOLD Response")
Now let's see what this plot looks like for a simple single stimulus.
In [ ]:
plot_imp_hrf([1])
We see that if we move around the time of the stimulus impulse, it also changes the timing of the expected BOLD response
In [ ]:
plot_imp_hrf([1])
plot_imp_hrf([8])
Now, plot the stimulus and output of a linear system that responds to both stimulus impulses.
In [ ]:
plot_imp_hrf([1, 8])
We will start with blocks of 2 events and then continue to longer blocks. We will examine how the number of events and the spacing between events affects the predicted bold signal.
Example 1: Two events that are spaced 4 seconds apart
In [ ]:
plot_imp_hrf([1, 5])
Example 2: Two events that are spaced 2 seconds apart
Note that here you get one peak and not 2 peaks. Why?
In [ ]:
plot_imp_hrf([1, 3])
Example 3: 3 events that are spaced 2s apart
In [ ]:
plot_imp_hrf([1, 3, 5])
Example 4: 5 stimuli given 2 seconds -> block which is 10 seconds long
What changed as you increased the number of stimuli? Why?
In [ ]:
plot_imp_hrf([1, 3, 5, 7, 9])
Researchers presented a subject with an image of a face and measured the hemodynamic response to this image. Then they showed the same stimulus again and found that the response peaked at the same time, but the amplitude of response was half the amplitude of the first presentation. Is the response for the second presentation expected from a linear system? Explain.
Researchers presented stimuli in 2 blocks. In the first block they presented flashing checkerboards at a rate of 1Hz for 10 seconds. In the second block they presented flashing checkerboards at a rate of 1hz for 20 seconds. How will the amplitude and/or duration of the block-response change between conditions? Explain.
Researchers presented stimuli in 2 conditions. In the first condition they presented a tone for 1 second, waited for 8 seconds, and presented a second tone for 1 second. In the second condition, they presented a tone for 1 second, waited 1 second, and presented a second tone for 1 second. How will the response differ in amplitude and/or duration between the two conditions? A plot/graph will be useful in explaining the answer.