2.1: Relevant Linear Systems Theory

Some common groundwork is necessary before investigating natural and biomimetic sensory systems and signal processing. A review of the salient aspects of linear systems (Section 2.1) is covered first. This is followed by fundamentals of neuronal systems (Section 2.2), neuronal processing (Section 2.3), an electric circuit model of a neuron (Section 2.4), and basic neuronal motion detection models (Section 2.5). The following free texts are recommended for more thorough treatment of linear systems theory and image processing:

Ulaby, F. and Yagle, A, Signals and Systems: Theory and Applications, Michigan Publishing, ISBN 978–1–60785–487–6, 2018. Available at ss2.eecs.umich.edu.

Yagle, A. and Ulaby, F., Image Processing for Engineers, Michigan Publishing, ISBN 978–1–60785–489–0, 2018. Available at ip.eecs.umich.edu.

2.1 Relevant Linear Systems Theory

This section summarizes the salient points of linear systems theory relevant to biomimetic sensory systems and signal processing. Many facets of natural vision processing can be modeled as spatial-temporal filters operating on input signals or image sequences. The two-dimensional spatial filters operate on each image frame, which is subsequently modified to account for the temporal history of previously filtered image frames. These operations are extensions of 1D discrete-time convolutions and other signal processing operations. The topics covered here include the motivation for LTI system modeling, continuous-time convolution and impulse response, δ(t), discrete-time convolution and unit pulse function, δ[n], 2D discrete-time convolution, and the Fourier Series and Fourier Transform.

2.1.1 Motivation for linear time-invariant (LTI) system modeling

Biological systems are naturally non-linear and time-varying. However, there is much practicality in approximating portions of the system as piece-wise linear over a nominal range of input values and time-invariant for relatively short periods of time. Such models many times can be "close enough" to be very useful. To simplify signal processing, we desire to use models of biological information transforms that are linear and time-invariant, or LTI transforms.

A system is linear if the law of superposition applies. In electrical engineering, superposition implies mathematical homogeneity and additivity. If the system output is y₁ = f(x₁) for an input x₁ and y₂ = f(x₂) for an input x₂, then the system is homogeneous if f(αx₁) = αy₁, additive if

f(x₁ + x₂) = f(x₁) + f(x₂) = y₁ + y₂

and therefore, linear if f(αx₁ + βx₂) = f(αx₁) + f(βx₂) = αy₁ + βy₂

A system is time-invariant if for the same response is given for the same set of inputs, regardless of when the inputs are presented. That is, for y(t) = f(x(t)), then f(x(t - t₀)) = y(t - t₀) for a constant time interval t₀.

Many natural (biological) signal processing functions can be modeled as a sequence of LTI subsystems we refer to as filters. The signal is sent through a filter and looks different at the output. For example, a filter simulating the layer of photoreceptors in mammalian vision systems may include a logarithmic conversion of light intensity followed by a blurring effect due to interactions with nearby photoreceptors. This model would include a logarithmic filter followed by a Gaussian blurring filter:

If h(t) is the impulse response to the system, then, for a continuous-time input signal f(t)

\(h(t) = f(t) * \delta (t) = \int_{0}^{t} f(\tau) \delta (t - \tau) \,d\tau\)

and for a discrete-time input sequence f[n]

\(h[n] = f[n] * \delta [n] = \sum_{k = -\infty}^{\infty} f[k] \delta [n - k]\)

The symbol * is used to denote convolution, and δ denotes the Direct delta function, both of which are discussed later. For the discrete-time case, time-invariance implies that if h[n] is the response to δ[n], then h[n-k] is the response to δ[n-k]. The impulse response (h[n] for discrete-time systems and h(t) for continuous-time systems) of a linear time-invariant (LTI) system completely characterizes that system. The output, y[n], of an LTI system is the sum of individual impulse responses weighted by the current input value, x[n]:

\(y[n] = \sum_{k = -\infty}^{\infty} x[k] h[n - k]\)

For a continuous-time signal, s(t), the Fourier Transform of s(t), represented as S(ω), shows the frequency content of s(t) as a function of radian frequency ω. The Inverse Fourier Transform of S(ω) gives back the original time-domain representation. The two functions, both representing the same signal from different domain perspectives form a Fourier Transform pair, denoted as: s(t) <==> S(ω).

The response of a system to an impulse is called the impulse response and is denoted as h(t). The frequency representation of the impulse response, H(ω), is called the frequency response. Given an input signal to a known LTI system, the output can be determined:

in the time domain as the convolution of input and LTI impulse response, or

in the frequency domain as the product of the input spectrum and the LTI frequency response.

To put it another way, given an input, x(t), whose spectrum is X(ω), to an LTI systems whose impulse response is h(t) (whose spectrum is H(ω), the frequency response) then the output, y(t), and its spectrum, Y(ω), is given as

Time domain: y(t) = x(t) * h(t)

Frequency domain: Y(ω) = X(ω) H(ω)

2.1.2 Continuous-time convolution

For continuous functions, f(t) and g(t), the convolution is defined by

\(f(t) * g(t) = \int_{0}^{t} f(\tau) g(t - \tau) \,d\tau\)

and has these properties:

f * g = g * f	commutative property
f * (g1 + g2) = f *g1 + f *g2	distributive property
(f * g) * v = f * (g * v)	associative property

2.1.3 Continuous-time unit impulse function

The unit impulse function or Dirac delta function is defined as

δ(t) = ∞, t = 0 δ(t-a) = ∞, t = a

= 0, t ≠ 0 = 0, t ≠ a

and

\(\int_{-\infty}^{\infty} \delta(t) \,dt =\int_{-\infty}^{\infty} \delta(t - a) \,dt = 1\)

since δ(t-a) has infinite strength at t = a, has zero duration at t = a, and has unity area.

A continuous-time signal x(t) can be replicated by the convolution with the unit impulse function:

\(x(t) = \int_{-\infty}^{\infty} x(\tau) \delta (t - \tau) \,d\tau \)

2.1.4 Discrete-time unit pulse function

A continuous-time signal can be discretized into a sequence, x[0], x[1], x[2] ..., by collecting the values of this integral at the specific times t₀, t₀+T, t₀+2T, t₀+3T, etc. so that

x[0] = x(t)*δ(t₀) = x(t₀)

x[1] = x(t)*δ(t₀+T) = x(t₀+T)

x[2] = x(t)*δ(t₀+2T) = x(t₀+2T)

_{: : :}

x[n] = x(t)*δ(t₀+nT) = x(t₀+nT)

Brackets, [ ], are used to denote discrete-time sequences while parentheses, ( ), are used to denote continuous-time functions. In this example, the sequence x[0], x[1], ...x[n]... is a discrete-time representation of the continuous-time function x(t).

The discrete-time version of the unit impulse sequence is called the unit sample sequence or the unit pulse function. The unit pulse function is defined as

δ[n] = 1, n = 0,

0, n ≠ 0

To properly interpret between continuous-time functions and discrete-time sequences, it should be noted that the unit pulse function must have unity area. To maintain this definition for a pulse of duration T, then the amplitude must be defined as 1/T. The definition above presumes T = 1.

A discrete-time signal x[n] can be replicated by the convolution with the unit pulse function:

\(x[n] = x[n] * \delta [n] = \sum_{k = -\infty}^{\infty} x[k] \delta [n - k]\)

2.1.5 One-dimensional discrete-time convolution

The discrete-time convolution operation is defined as

\(f[n] * g[n] = \sum_{k = -\infty}^{\infty} f[k] g[n - k]\)

Equation 2.1-1

Example 2.1-1

Use Equation 2.1-1 to calculate y[n] = f[n] * g[n] where f[n] = {1 2 3} and g[n] = {1 1 1 1 1 1}.

Example 2.1-1

Use Equation 2.1-1 to calculate y[n] = f[n] * g[n] where f[n] = {1 2 3} and g[n] = {1 1 1 1 1 1}.

Solution:

f[n] = {1 2 3} implies f[0] = 1, f[1] = 2, f[2] = 3 and

g[n] = {1 1 1 1 1 1} implies g[0] = g[1] = g[2] = g[3] = g[4] = g[5] = 1

For n = 0, Equation 2.1-1 gives the summation y[0] = …f[−1]g[1] + f[0]g[0] + f[1]g[−1] + …

but all values of f[n] and g[n] are zero for n < 0, so y[0] = f[0]g[0] = 1.

\(y[0] = \sum_{k = 0}^{5} f[k] g[0 - k] = f[0]g[0] = 1\)

\(y[1] = \sum_{k = 0}^{5} f[k] g[1 - k] = f[0]g[1] + f[1]g[0] = 1 + 2 = 3\)

\(y[2] = \sum_{k = 0}^{5} f[k] g[2 - k] = f[0]g[2] + f[1]g[1] + f[2]g[0] = 1 + 2 + 3 = 6\)

\(y[3] = \sum_{k = 0}^{5} f[k] g[3 - k] = f[0]g[3] + f[1]g[2] + f[2]g[1] = 1 + 2 + 3 = 6\)

\(y[4] = \sum_{k = 0}^{5} f[k] g[4 - k] = f[0]g[4] + f[1]g[3] + f[2]g[2] = 1 + 2 + 3 = 6\)

\(y[5] = \sum_{k = 0}^{5} f[k] g[5 - k] = f[0]g[5] + f[1]g[4] + f[2]g[3] = 1 + 2 + 3 = 6\)

\(y[6] = \sum_{k = 0}^{5} f[k] g[6 - k] = f[1]g[5] + f[2]g[4] = 2 + 3 = 5\)

\(y[7] = \sum_{k = 0}^{5} f[k] g[7 - k] = f[2]g[5] = 3\)

Where products are omitted when they are zero for specific k values. The result is y[n] = { 1 3 6 6 6 6 5 3 }.

■

We could also visualize the answer graphically by reversing the sequence g[n] and placing it below f[n] and offsetting by the value of n in Equation 2.1-1. The first three values determined using this graphical approach are shown here:

n =0:

k: -5 -4 -3 -2 -1 0 1 2 3 4 5

f[k]: 0 0 0 0 0 1 2 3 0 0 0

g[n-k]: 1 1 1 1 1 1 0 0 0 0 0

\(y[0]=\sum_{k=0}^{5} f[k] g[0-k]=f[0] g[0]=1\)

n =1:

k: -5 -4 -3 -2 -1 0 1 2 3 4 5

f[k]: 0 0 0 0 0 1 2 3 0 0 0

g[n-k]: 0 1 1 1 1 1 1 0 0 0 0

\(y[1]=\sum_{k=0}^{5} f[k] g[1-k]=f[0] g[1] + f[1] g[0]=1 + 2 = 3\)

n =2:

k: -5 -4 -3 -2 -1 0 1 2 3 4 5

f[k]: 0 0 0 0 0 1 2 3 0 0 0

g[n-k]: 0 0 1 1 1 1 1 1 0 0 0

\(y[2] = \sum_{k = 0}^{5} f[k] g[2 - k] = f[0]g[2] + f[1]g[1] + f[2]g[0] = 1 + 2 + 3 = 6\)

And so on.

Example 2.1-2

Solve the same convolution (Example 2.1-1) using Equation 2.1-1 but reverse the roles of f[n] and g[n].

Solution:

Due to the commutative property the result should be the same. Using Equation 2.1-1, the order could also be reversed, as y = g * f, giving:
\(y[0]=\sum_{k=0}^{5} g[k] f[0-k]=g[0] f[0]=1 \)
\(y[1]=\sum_{k=0}^{5} g[k] f[1-k]=g[0] f[1]+g[1] f[0]=2+1=3 \)
\(y[2]=\sum_{k=0}^{5} g[k] f[2-k]=g[0] f[2]+g[1] f[1]+g[2] f[0]=3+2+1=6 \)
\(y[3]=\sum_{k=0}^{5} g[k] f[3-k]=g[1] f[2]+g[2] f[1]+g[3] f[0]=3+2+1=6 \)
\(y[4]=\sum_{k=0}^{5} g[k] f[4-k]=g[2] f[2]+g[3] f[1]+g[4] f[0]=3+2+1=6 \)
\(y[5]=\sum_{k=0}^{5} g[k] f[5-k]=g[3] f[2]+g[4] f[1]+g[5] f[0]=3+2+1=6 \)
\(y[6]=\sum_{k=0}^{5} g[k] f[6-k]=g[4] f[2]+g[5] f[1]=3+2=5 \)
\(y[7]=\sum_{k=0}^{5} g[k] f[7-k]=g[5] f[2]+g[5] f[1]=3 \)

so that once again y[n] = { 1 3 6 6 6 6 5 3}.

■

As before we could also visualize the answer graphically, but this time by reversing the sequence f[n] and placing it below g[n] and offsetting by the value of n in Equation 2.1-1. The first three values determined using this graphical approach are shown here:

n =0:

k: –2 –1 0 1 2 3 4 5 6

g[k]: 0 0 1 1 1 1 1 1 0

f[n-k]: 3 2 1 0 0 0 0 0 0

\(y[0]=\sum_{k=0}^{5} g[k] f[0-k]=g[0] f[0]=1\)

n =1:

k: –2 –1 0 1 2 3 4 5 6

g[k]: 0 0 1 1 1 1 1 1 0

f[n-k]: 0 3 2 1 0 0 0 0 0

\(y[1]=\sum_{k=0}^{5} g[k] f[1-k]=g[0] f[1]+g[1] f[0]=2+1=3 \)

n =2:

k: –2 –1 0 1 2 3 4 5 6

g[k]: 0 0 1 1 1 1 1 1 0

f[n-k]: 0 0 3 2 1 0 0 0 0

\(y[2]=\sum_{k=0}^{5} g[k] f[2-k]=g[0] f[2]+g[1] f[1]+g[2] f[0]=3+2+1=6 \)

Convolutions can be computed quickly in MatLab, which stands for “matrix laboratory”, a software product of Mathworks, Inc. [MatLab]. The sequence y = f * g {implying y[n] = f[n]*g[n]} as defined in Example 2.1-1 is computed in Matlab as

>> f = [1 2 3];
>> g = [1 1 1 1 1 1];
>> y = conv(f,g)

y = 1 3 6 6 6 6 5 3

Note that the resulting sequence is longer that both input sequences; this is a natural consequence of the convolution operation. Also, MatLab does not use a zero^th element; the first element of y is referenced in MatLab as y[1], which is the same value as y[0] in Equation 2.1-1. Given these considerations the simulation confirms the hand-written results in Example 2.1-1.

Exercise 2.1-1

Use Equation 2.1-1 to calculate y[n]= x[n] * h[n] for x[n] = {1 2 2} and h[n] = {1 3}. Check your answer using graphical convolution.

Answer: y[n] = {1 5 8 6}

■

Exercise 2.1-2

Use Equation 2.1-1 to calculate y[n]= x[n] * h[n] for the given sequences. Check your answer using graphical convolution.

1. x[n] = {1 1 1 1 1} and h[n] = { 0.25 0.5 0.25}.

2. x[n] = {1 1 1 1 1} and h[n] = { 0.25 -0.5 0.25}.

Answers:

1. y[n] = {0.25 0.75 1.0 1.0 1.0 0.75 0.25}

2. y[n] = {-0.25 0.25 0 0 0 0.25 -0.25}

■

2.1.6 Two-dimensional discrete-time convolution

In signal processing applications, one sequence may represent a filter and another a given input signal. Convolving the two signals would give an output that is a filtered representation of the given input signal. Image processing is two-dimensional (2D) signal processing, where one 2D signal (image) represents a filter and the other the input image. For biomimetic applications, the filter could represent a model of a natural phenomenon as it affects the raw imagery. Assuming 2D variables f(x,y) and g(x,y), the 2D convolution operation is given as:

\(f[x, y] * g[x, y]=\sum_{m=0}^{M-1} \sum_{n=0}^{N-1} f[m, n] g[x-m, y-n]\)

Example 2.1-3

Using MatLab computer the 2D convolution y = f * g given the 2D variables f and g defined here:

f = 1 1 g = 2 2 2

1 1 2 2 2

2 2 2

Solution:

>> f = [1 1; 1 1];

>> g = [2 2 2; 2 2 2; 2 2 2];

>> y = conv2(f,g) y = 2 4 4 2

4 8 8 4

4 8 8 4

2 4 4 2

■

Notice that for both 1D and 2D convolutions the result tends to grow. Sometimes it is necessary to ‘crop’ out the internal pixels so that the filtered version is the same size as the original. For example, defining f below as a 5x5 2D unit pulse and performing the convolution with the 3x3 2D variable g results in a 7x7 result:

>> f = zeros(5);

>> f(3,3) = 1; f = 0 0 0 0 0

0 0 0 0 0

0 0 1 0 0

0 0 0 0 0

0 0 0 0 0

>> g = [1 2 3; 4 5 6; 7 8 9]; g = 1 2 3

4 5 6

7 8 9

>> y=conv2(f,g); y = 0 0 0 0 0 0 0

0 0 0 0 0 0 0

0 0 1 2 3 0 0

0 0 4 5 6 0 0

0 0 7 8 9 0 0

0 0 0 0 0 0 0

0 0 0 0 0 0 0

To visualize the computation, flip g both vertically and horizontally:

g^* = 9 8 7

6 5 4

3 2 1

The values for y are found by placing g* over each element of f, and then performing a dot product; that is, multiply element-for-element, and then add the products (sum of products). One way to crop out the middle 5x5 is to manually redefine y:

>> y = y(2:6,2:6);

;This command redefines y to be row 2 through 6 and column 2 through 6. The result is cropping the internal 5x5 image of the previous convolution so that the filtered version is now the same size as the original image, f.

y = 0 0 0 0 0

0 1 2 3 0

0 4 5 6 0

0 7 8 9 0

0 0 0 0 0

A better alternative for cropping out the central pixels is to use the same attribute in the conv2 Matlab function when determining y:

>> y=conv2(f,g,’same’); y = 0 0 0 0 0

0 1 2 3 0

0 4 5 6 0

0 7 8 9 0

0 0 0 0 0

Exercise 2.1-3

Given the 2D filter f and image x, give the output filtered image y = f * x.

f = -2 0 2 x = 0 0 0

-1 0 1 0 1 1

0 1 0 0 0 0

Give both the cropped and uncropped representations of the output.

Answers:

Uncropped: y = 0 0 0 0 0 Cropped: y = -2 -2 2

0 -2 -2 2 2 -1 -1 1

0 -1 -1 1 1 0 1 1

0 0 1 1 0

0 0 0 0 0