Linear Systems Analysis - Vision Systems - Lecture Notes, Study notes for Ophthalmology

21 pages
1000+Number of visits
Description
Vision Systems lecture handout. Keywords of this handout are: Linear Systems Analysis, Constraints Defining Linear Systems, Impulse Response Functions, One Dimensional Systems, Two Dimensional Systems, Impulse Response F...
20 points
this document
Preview3 pages / 21
Microsoft Word - Vision Notes.doc

57

IV. Linear systems analysis

Linear systems analysis refers to a set of mathematical techniques that can be used to analyze and describe input/output systems that satisfy certain linearity assumptions. The inputs and outputs of linear systems are represented by functions. These functions can be of any number of variables, but for present purposes there will only be one variable, time (t) or space (x), or two variables, space-space (x,y), or space-time (x,t). For example the input functions might represent the variation in light level on a single receptor over time, the spatial distribution of light entering the eye, or the spatial distribution of neural activity transmitted from the eye to the brain along the optic nerve. The output functions might represent the electrical response of the receptor over time, the spatial distribution of light falling on the retina, or the spatial distribution of neural activity of a group of cells receiving inputs from the optic nerve.

The letter L will be used to represent a linear system; i(t), i(x,y), etc. will represent input functions; and o(t), o(x,y), etc. will represent output functions. Thus, we can write in the one-dimensional case, o(t) = L i(t){ } (4.1) or in the two-dimensional case, o(x,y) = L i(x,y){ } (4.2) The curly brackets are used as a reminder that a linear system takes a whole function as input and gives a whole function as output. For example, the optical system of the eye can be regarded as a linear system that transforms the two-dimensional light distribution in an object plane, i(x,y), into a two-dimensional light distribution in the image plane, o(x, y).

Figure 4.1 illustrates graphically a one-dimensional input function, a linear system, and a one dimensional output function. To make the figure concrete the reader might imagine that the left side represents a temporal signal (rectangular pulse) turned on shortly after time 0, while the right side represents the response of the linear system to this signal.

A. Constraints defining linear systems

In order for a one-dimensional system to be linear, it must satisfy the following constraint, L ai1(t) + bi2 (t){ }= aL i1(t){ }+ bL i2(t){ } (4.3)

Docsity.com

58

where i1(t) and i2(t) are two arbitrary input functions and a and b are arbitrary constants (scalars). In words, the constraint is that the response of a linear system to the sum of two input function equals the sum of the responses to each input taken alone, and also scaling an input function by some amount scales the output function by exactly the same amount. In Figure 4.1B, another pulse has been added to the one in Figure 4.1A. If the system is linear then the response must be the sum of the individual responses (as shown). In Figure 4.1C, the input pulse in Figure 4.1A has been doubled in amplitude. If the system is linear the response should double in amplitude (as shown).

Figure 4.1

The constraint is essentially the same for two-dimensional stimuli, L ai1(x, y) + bi2 (x, y){ }= aL i1(x, y){ }+ bL i2 (x, y){ } (4.4) where i1(x, y) and i2(x, y) are two arbitrary input functions and a and b are arbitrary constants (scalars).

A somewhat more constrained type of system (the one that we will focus on) is the linear shift-invariant (LSI) system. Shift invariant means that shifting the input along the input dimensions shifts the output by a corresponding amount without changing the output shape. In other words, for a system to satisfy this constraint, a shift in the starting time or position of an input should not change the shape of the output. More formally, if

Docsity.com

59

o(t) = L i(t){ } (4.5) then o(t t0 ) = L i(t t0 ){ } (4.6) where t0 is an arbitrary shift along the t axis. In the two dimensional case, if, o(x,y) = L i(x,y){ } (4.7) then o(x x0, y y0 ) = L i(x x0, y y0 ){ } (4.8) where x0 and y0 are arbitrary shifts along the x-axis and y-axis, respectively. (Note that the linear system may perform, among other things, an arbitrary linear transformation of the axes.) Many systems that are not shift invariant globally can be approximated and analyzed as LSI systems in local regions.

In Figure 4.1D the pulse in Figure 4.1A has been shifted in time. If the system is shift invariant the output should be shifted without a change in shape (as shown).

Any system that satisfies the above constraints can be studied and characterized using the linear-systems techniques described below. However, the above constraints are very strong; thus, it is important to keep in mind that there are many systems that cannot be analyzed meaningfully with these techniques. Nonetheless, there a good number of situations where the constraints hold accurately and many others where they hold well enough for a linear-systems analysis to be useful.

B. Impulse-response functions and convolution

If a system is linear and shift invariant it can be completely characterized by its response to an instantaneous pulse of finite area (for one-dimensional systems) or of finite volume (for two-dimensional systems). This instantaneous pulse is known as an impulse function or delta function. The response to an impulse of unit area or volume is called the impulse response function. It is of course impossible to produce pulses of zero duration and an area or volume of one. However, for every physically realizable linear system, it is possible to select the pulse duration short enough so that the response is essentially identical to the impulse-response function. The use of impulse functions is a mathematical convenience.

Docsity.com

60

An impulse function is represented graphically as an arrow with height equal to its area or volume. The arrow on the left in Figure 4.2A represents an impulse function of unit area presented a location 0. The function labeled h(t) on the right represents an impulse-response function. The impulse-response function completely characterizes a LSI system because the response of the system to any input can be predicted from knowledge of the impulse-response function alone. This fact rests on the convolution formula, which will now be derived in an intuitive (not rigorous) fashion.

Figure 4.2

1. Convolution formula for one-dimensional systems

Consider an arbitrary input function, i(t), to a LSI system with impulse response function h(t). This arbitrary input is illustrated in Figure 4.2B. The first thing to note is that the input function can be broken down into the sum of short slices of width ∆α. One of these slices (the one at position α) has been shaded in the figure. As illustrated just below this slice, the effect of the slice on the system can be approximated by an impulse function positioned at α with a of height i(α) ∆α (the approximate area of the slice). The smaller is ∆α, the smaller the slices and the more accurate will be the approximation. What is the response to this slice? Because of shift invariance (Figure 4.1D), the response must have the same shape as what would have been obtained if the slice were centered at the origin (0); but it will be shifted over by α. Because of the scaling constraint (Figure 4.1c), the shape (if presented at the origin) must be that of the impulse- response function but scaled (multiplied) by i(α) ∆α (the area of the impulse). Thus, the response to the slice will be approximately h(t-α) i(α) ∆α, as illustrated in the right side of Figure 4.2B. This is the response to the slice at an arbitrary position α. Because of the

Docsity.com

61

summation constraint (Figure 4.1B) the response to the sum of the slices (the whole input) is the sum of the responses to each slice. Thus,

( ) ( ) ( )j j j

o t i h tα α α ∞

=−∞

= − ∆∑ (4.9) Intuitively, as ∆α gets smaller the slices are better approximated by impulse functions and the summation becomes more accurate. Thus, in the limit, as ∆α approaches 0.0,

( ) ( ) ( )o t i h t dα α α∞ −∞

= −∫ (4.10) This is the convolution formula and it can be used to compute the response to arbitrary inputs. The symbol * will be used to represent operation of one-dimensional convolution. Thus, the shorthand notation for equation (4.10) is

o(t) = h(t) * i(t)

(4.11)

2. Convolution formula for two-dimensional systems

Computing the response of a two-dimensional LSI system from its impulse response function is entirely analogous. Let the arbitrary input function be i(x, y) and the impulse response function be h(x, y), then the output is given by,

( ) ( ) ( ), , ,o x y i h x y d dα β α β α β∞ ∞ −∞ −∞

= − −∫ ∫ (4.12) The "intuitive" derivation of this formula is analogous to the one-dimensional case and is left for the reader. The symbol ** will be used to indicate the operation of two- dimensional convolution. Thus, equation (4.12) can be represented as

( ) ( ) ( ), , ,o x y h x y i x y= ∗∗ (4.13)

3. Measuring the impulse-response function

The beauty of a linear system is that once the impulse response function is known the response of the system to any input can be predicted. There are several common methods of measuring the impulse response function. The most direct method is simply to present an impulse and measure the system's response. As mentioned above, it is usually possible to approximate an impulse with a narrow pulse (perfect impulses are not

Docsity.com

62

physically realizable). A second method is to measure the response to edges (step- or edge-response functions). When the system is two-dimensional the responses to edges at several orientations must be measured. A third method that is often used with two- dimensional systems is to measure the response to lines of various orientations (line- response functions). Finally, a fourth and very useful method is to measure the response to sinusoidal functions. This last method will be described in more detail later in the discussion of Fourier analysis.

For all but the first method, the responses to the inputs must be mathematically transformed to obtain the impulse response function. These other methods are used because they tend to yield more accurate (less noisy) measurements of the impulse response function.

C. Linear weighting functions and cross correlation

One very common operation used in image processing and in modeling biological vision is the cross-correlation of signals with a weighting function. The operation of cross correlation consists of multiplying the signal by the weighting function and summing (integrating) the result, for each possible location of the weighting function. In image-processing applications, the weighting function might be a spatial filter shape (like a difference of Gaussians); in biological-modeling applications, the weighting function might represent the shape of a spatial receptive field or the shape of a neural temporal integrator.

1. Cross-correlation formula for one-dimensionalsystems

Figure 4.3 illustrates the operation of cross-correlation for the one-dimensional case. To derive the cross-correlation formula, let i(t) represents an arbitrary input signal, w(t) represents the weighting function (when it is positioned at the origin), and o(t) represents the output function. To compute the output at a particular time, t, the weighting function is shifted over to position t and multiplied by the input signal. Let α be an arbitrary point along the t axis. The product at this point is i(α) w(α-t). Summing across all values of α we have:

( ) ( ) ( )o t i w t dα α α∞ −∞

= −∫ (4.14) The symbol ⊗ is used to represent the operation of cross correlation, thus the short-hand notation for equation (4.14) is

o(t) = i(t) ⊗ h(t)

(4.15) It can be shown that any system that produces its output by applying a single weighting function is a linear shift invariant system. (This can be shown by checking that equation

Docsity.com

63

(4.14) satisfies the properties of a LSI system.) Indeed, equation (4.14) is seen to be very similar to equation (4.10). Note that equation (4.14) can be written as

( ) ( ) ( )( ) ( ) ( )o t i w t d i t w tα α α∞ −∞

= − − = ∗ −∫ (4.16) Thus, the impulse response function of a cross-correlation system is just the system weighting function with the sign of the argument reversed; h(t) = w(-t). For symmetrical weighting functions, the impulse response and weighting functions are identical.

Figure 4.3

2. Cross-correlation formula for two-dimensional systems

The formula for two-dimensional cross correlation is entirely analogous:

( ) ( ) ( ), , ,o x y i w x y d dα β α β α β∞ ∞ −∞ −∞

= − −∫ ∫ (4.17) The symbol ⊗⊗ is used to represent the two-dimensional cross correlation. Also the impulse response function of a 2D cross-correlation system is the weighting function with the sign of the arguments reversed; h(x, y) = w(-x,-y). Thus,

o(x,y) = i(x,y) ⊗⊗ w(x,y) = i(x,y) ** w(-x,-y) (4.18)

Docsity.com

64

3. Measurement of the weighting function

Because the weighting function of a linear system is just the impulse-response

function with the sign of the arguments reversed, the weighting function can be found by measuring the impulse response function or vice versa. The only thing new worth noting is that weighting function is most directly measured by recording the output of the system at a fixed point while presenting input impulses (edges, lines or sinewaves) at different positions. Sometimes this is more convenient than measuring the impulse response function (i.e., the output at different points to an input presented at just one point).

D. The Fourier transform

The two major steps in applying a linear systems analysis are (1) measuring the impulse-response function (or weighting function) and (2) computing the predicted outputs for different inputs. The methods of Fourier analysis and synthesis often prove useful for carrying out both steps. A brief overview of the Fourier methods will now be given before presenting some of the mathematical details.

The fundamental concept in Fourier analysis and synthesis is that any function can be

decomposed into a sum of sinusoidal functions of differing frequencies, amplitudes, and phases (positions). This decomposition is called the Fourier transform. Formally, the Fourier transform is a complex valued function, which can be regarded as a pair of real- valued functions: the amplitude spectrum (sinewave amplitude as a function of frequency) and the phase spectrum (sinewave phase as a function of frequency). In other words, the Fourier transform of a function is essentially a pair of functions describing the sinewaves that when added together would give back the function.

The inverse Fourier transform adds up the all the sinewaves described by the

amplitude and phase spectra. In other words, the inverse Fourier transform takes the amplitude and phase spectra as input and returns the original function.

Fourier analysis and synthesis are demonstrated in Figure 4.4. The straight-edged

function at the bottom is a square wave and it (like any other function) can be decomposed into sine-wave functions. If the repetition frequency of the square wave (i.e., the number of squares per unit distance or time) is ω then the sine-wave frequencies that are needed to reconstruct the square wave are ω, 3ω, 5ω, 7ω, … and so on. The first component is the fundamental and the rest of the components are harmonics of the fundamental. (The harmonics of a frequency are the frequencies that are integer multiples of that frequency. By this definition the fundamental is also the first harmonic.) A square-wave contains only the odd numbered harmonics. In other words, the amplitude and phase spectra of a square wave are zero everywhere except at the frequency of the square wave and at the odd harmonics of that frequency. The first several sine-wave components are shown at the top of the figure along with the sum of

Docsity.com

65

these components. As can be seen, the sum does approximate the square wave. Adding more components would give a better approximation, but an infinite number would be required for an exact synthesis.

Figure 4.4

Fourier analysis and synthesis of two-dimensional functions, such as visual image, requires using two-dimensional sinewaves, like the one illustrates in Figure 4.5A. Two dimensional sine-wave patterns that modulate in luminance (or radiance) are called sine- wave gratings. As shown in Figure 4.5A, these two-dimensional sine waves have both a vertical and a horizontal spatial frequency. In other words, the luminance modulates with one frequency in the vertical direction and another frequency in the other direction. Vertically oriented gratings have a spatial frequency of zero in the vertical direction, horizontally oriented gratings a frequency of zero in the horizontal direction, and all other gratings have nonzero frequencies in both directions. (A slice of any orientation through a two-dimensional grating has a sinusoidal profile.) Thus, for two-dimensional functions the amplitude and phase spectra are functions of both vertical and horizontal frequency. Figure 4.5 B demonstrates Fourier analysis and synthesis of a vertical square-wave grating using vertical sine-wave gratings. The waveform at the bottom of the figure is the sum of the sine-wave gratings above. The Fourier representation of something complex like a face will contain components of all different orientations (i.e., all different pairs of horizontal and vertical frequency).

The Fourier transform of the impulse-response function of a linear system is called as the transfer function (TF). The amplitude and phase spectra contained in the transfer function are called the amplitude transfer function (ATF) and the phase transfer function (PTF), respectively. The impulse-response function can be obtained from the transfer function by applying an inverse Fourier transform. Thus, the impulse-response function

Docsity.com

66

can be measured by measuring the transfer function, and then applying an inverse Fourier transform. Measuring the transfer function often gives a much more accurate estimate of the impulse response function than any other method. The direct method for measuring the transfer function is based on the following facts of Fourier theory: (a) A sine-wave input to an LSI system produces a sine-wave output that can differ from the input only in amplitude and phase (it must of the same frequency). (b) The amplitude spectrum of the output of an LSI system (Ao) is the product of the amplitude spectrum of the input (Ai) and the amplitude spectrum of the impulse-response function (Ah, which is the ATF):

Ao = Ah Ai

(4.19) (c) The phase spectrum of the output (Po) is the sum of the phase spectrum of the input (Pi) and the phase spectrum of the impulse-response function (Ph, which is the PTF):

Po = Ph + Pi

(4.20)

These properties imply that one can measure the transfer function of a linear system by presenting sine-wave inputs of all different frequencies. If the amplitudes and phases of these sine-wave inputs and outputs are recorded at each frequency, then the amplitude transfer function can be obtained by dividing the output amplitudes by the input amplitudes:

Ah = Ao / Ai

(4.21) and the phase transfer function can be obtained by subtracting the input phases from output phases:

Ph = Po - Pi

(4.22)

Properties (b) and (c) are also useful for computing the output of the linear system once the impulse response function is known. Specifically, one can take the Fourier transform of the input, and then use the transfer function and equations (4.19) and (4.20) to obtain the Fourier transform of the output. A final inverse Fourier transformation yields the output signal. This circuitous route through the Fourier domain is usually faster on a digital computer than directly calculating the convolutions.

Docsity.com

67

Figure 4.5

A

B

Docsity.com

68

A few mathematical details concerning Fourier transforms of real-valued functions

will now be given. Many of the properties described below also hold for complex-valued functions, but we will not need to deal much with complex functions here. Also, the intuitive concepts are easier to grasp for real functions. More complete treatments are available in Bracewell (1978) and Gaskill (1978).

1. One-dimensional systems

The Fourier transform, F(u), of a real-valued function, f(t), is defined by the following equation:

( ) ( ) 2j utF u f t e dtπ∞ − −∞

= ∫ (4.23) where u is frequency (e.g., cycles per second) and j is the imaginary constant -1. The inverse Fourier transform is given by,

( ) ( ) 2j utf t F u e duπ∞ −∞

= ∫ (4.24) This can be verified by substituting equation (4.18) into equation (4.19) (see Bracewell, 1978 or Gaskill, 1978).

The term e-j2πut is a complex sine wave. It is called this because of the following equality,

e-j2πut = cos(2πut) - j sin(2πut)

(4.25) Equation (4.25) is derived by expanding both sides into Taylor's series and comparing terms. [Complex numbers are defined as pairs of numbers (e.g., (a, b)) where the first number is the real part of the complex number and the second is the imaginary part. Complex numbers are also written in the form a + jb. This latter form is helpful in doing algebra with complex numbers. The basic properties and algebraic rules for complex numbers can be found in Gaskill and in most college algebra textbooks.]

Substituting equation (4.25) into equation (4.23) gives,

( ) ( ) ( ) ( ) ( )cos 2 sin 2F u f t ut dt j f t ut dtπ π∞ ∞ −∞ −∞

= −∫ ∫ (4.26)

Docsity.com

69

The real part is called the cosine transform (C(u)), the imaginary part the sine transform (S(u)):

( ) ( ) ( )cos 2C u f t ut dtπ∞ −∞

= ∫ (4.27)

( ) ( ) ( )sin 2S u f t ut dtπ∞ −∞

= ∫ (4.28) and thus,

F(u) = C(u) - jS(u)

(4.29) C(u) is always symmetric about the origin (u=0) and S(u) is always anti-symmetric about the origin.

This shows that the Fourier transform of a real-valued function consists of two real- valued transforms (in fact, the complex numbers play a relatively minor role). In words, the Fourier transform of a real function is obtained by summing (integrating) the product of the function and a unit amplitude cosine function, and by summing the product of the function and a unit amplitude sine function. This calculation is done for each frequency (u), and thus the Fourier transform is a function of u.

Substituting equations (4.25) and (4.29) into the formula for the inverse Fourier transform (equation 4.24) gives some insight into what the Fourier transform does:2

( ) ( ) ( ) ( ) ( )cos 2 sin 2f t C u ut dt S u ut dtπ π∞ ∞ −∞ −∞

= +∫ ∫ (4.30) Equation (4.30) says that f(t) is equal to the sum of cosine functions and sine functions with amplitudes C(u) and S(u), respectively. In other words, C(u) and S(u) describe the sine and cosine waves that when added together reproduce the function.

We can also express C(u) and S(u) in polar coordinates using the inverse polar coordinate transformations:

( ) ( ) ( )( )cosC u A u P u= (4.31)

2 The cross product terms that appear after the substitutions, integrate to zero because f(t) is a real valued function and thus C(u) and cos(2πut) are even symmetric about the frequency origin and S(u) and sin(2πut) are odd symmetric about the origin. When an even function is multiplied by an odd function the result is odd, and odd functions must integrate to zero.

Docsity.com

70

and,

( ) ( ) ( )( )sinS u A u P u= (4.32) Substituting these into equation (4.30) and simplifying shows that

( ) ( ) ( )( )cos 2f t A u ut P u dtπ∞ −∞

= +∫ (4.33)

Equation (4.33) says that f(t) is also equal to the sum of cosine functions with amplitudes A(u) and phases P(u). In other words, the Fourier transform of a function gives the cosine waves (of various amplitudes and phases) that when added together reproduce the function. A(u) is called the amplitude spectrum, and P(u) the phase spectrum. A(u) is always symmetric about the origin and P(u) is always anti-symmetric about the origin. Figure 4.6 shows a one-dimensional function and its amplitude spectrum.

Most computer subroutine for doing Fourier transforms give C(u) and S(u) as outputs. The amplitude and phase spectra can be obtained (if desired) using the polar coordinate transformations:

A(u) = C(u)2 + S(u)2

(4.34) and,

P(u) = arctan(S(u)

C(u) )

(4.35)

Figure 4.6

Docsity.com

71

Figure 4.7a

Docsity.com

72

Figure 4.7b

The Fourier transform can also be expressed in terms of the phase and amplitude spectra by substituting the inverse polar coordinate transformations (4.31 and 4.32) into equation (4.29):

F(u) = A(u) ejP(u)

(4.36) A table of properties of the one-dimensional Fourier transform (from Gaskill, 1978) is given in Figure 4.7. The only one that will receive special mention now is the Fourier convolution property: the Fourier transform of a convolution of two functions [ i(t)*h(t) ] equals the product of the Fourier transforms of the individual functions [ I(u)H(u) ]. This property is easily proved (see Gaskill). Given the convolution property, equation (4.36) implies that,

Docsity.com

73

Ao(u) = Ah(u) Ai(u)

(4.37) and,

P0 (u) = Ph (u) + Pi (u)

(4.38) where Ao, Ai and Ah are the amplitude spectra of o(t) (the result of the convolution), i(t) and h(t), respectively, and Po, Pi and Ph are the corresponding phase spectra. As mentioned earlier, these properties are useful for measuring the impulse-response function of a linear system and for computing its output to different inputs (c.f., equations 4.21 - 4.22).

2. Two-dimensional systems

The definition of the Fourier transform for two-dimensional (2D) functions is analogous to that for one-dimensional signals:

( ) ( ) ( )2, , j ux vyF u v f x y e dxdyπ∞ − + −∞

= ∫ (4.39) where u is frequency in one direction (e.g., cycles per horizontal degree), v is frequency in the other direction (e.g., cycles per vertical degree), and j is the imaginary constant 1− . The inverse two-dimensional Fourier transform is given by,

( ) ( ) ( )2, , j ux vyf x y F u v e dudvπ∞ + −∞

= ∫ (4.40) As in the one-dimensional case, the two-dimensional Fourier transform can be expressed in terms of the cosine and sine transforms,

F(u,v) = C(u,v) - j S(u,v)

(4.41) or in terms of the amplitude and phase spectra,

F(u,v) = A(u,v) ejP(u,v)

(4.42)

Docsity.com

74

C(u, v) and A(u, v) are both symmetric about the origin (in the two-dimensional plane) and S(u, v) and P(u, v) are both anti-symmetric about the origin.3 Thus, the Fourier transform is completely specified by its values on one half of the two-dimensional (u, v) plane. Figure 4.8 shows a two-dimensional function and its amplitude spectrum.

A table of properties of the two-dimensional Fourier transform (from Gaskill, 1978) is given in Figure 4.9. Comparison with Figure 4.6 shows that the properties of the 2D transform are very similar to those of the 1D transform. Importantly, the convolution property also holds in the 2D case: the Fourier transform of a convolution of two functions [i(x, y)**h(x, y)] equals the product of the Fourier transforms of the individual functions [I(u, v) H(u, v)]. Thus,

( ) ( ) ( )0 , , ,h iA u v A u v A u v= (4.43)

and

( ) ( ) ( )0 , , ,h iP u v P u v P u v= + (4.44) where Ao, Ai and Ah are the amplitude spectra of o(x, y) (the result of the convolution), i(x, y) and h(x, y), respectively, and Po, Pi and Ph are the corresponding phase spectra. As mentioned earlier, these properties are useful for measuring the impulse-response function of a 2D linear system and for computing its output to different inputs (c.f., equations 4.14 - 4.17).

Figure 4.8

3 If a two dimensional function f(x, y) is symmetric about the origin then f(x, y) = f(-x,-y); if it is antisymmetric then f(x, y) = -f(-x,-y).

Docsity.com

75

Figure 4.9 A

Docsity.com

76

Figure 4.9a (continued)

Docsity.com

77

Figure 4.9b

Docsity.com