PhD Thesis: Chapter 2

 

"Phase-Only Optical Information Processing"

University of Edinburgh, D.J.Potter, 1992.

Index   Chapter 1  2  3  4  5  6  7  8  9 

Chapter 2
Fourier Series Phase Object Spectra

Introduction

In chapter one, it was noted how one may describe a general phase object with optical path variations described by f(x) as the argument in the exponential e^i f(x). Zernike performed a simple Taylor expansion of the exponential which allowed him to obtain a relatively simple expression for the phase object spectrum. The validating qualification being that powers of f(x) greater than unity be negligible relative to the linear term i.e. f(x) be very small. This is known as the `weak phase' approximation.

In this chapter, an alternative perspective with which to view the whole subject of phase object spectra is introduced. The underlying principle at work in the formation of a phase object spectrum, it will be shown, is one of convolution. By developing a mathematical framework based on convolution, several interesting observations emerge as to both the range of validity of the `weak phase' approximation and the assumptions of image intensity linearity with object phase after phase contrast filtering. These form the subject of chapter three for which this chapter specifies the algorithm used in the calculation of the phase object spectrum.

1  The Nature of the Problem

For simplicity the problem is considered in one dimension in all that follows. Let the object plane of an optical processor contain a general phase object which results in a complex transmission function of
where a Taylor expansion of the exponential is utilised and f(x) is assumed to be small. Upon Fourier Transformation the light distribution on the back focal plane of the transform lens has a field given by
where F(n) is the Fourier Transform of f(x) and d denotes the Dirac delta-function. Many image processing operations seek to obtain a visual representation of the underlying phase structure f(x), and some (phase-contrast imaging, for example) ideally seek an image intensity which is linearly proportional to f(x).

Although strictly a non-linear transformation, utilising the Taylor expansion of e^i f(x) shows it to be an approximately linear process if the argument of the exponential (the object phase) is very small. Where the small phase approximation is valid, it is allowable to describe a spectrum transfer function whose sole effect is to pass the spectrum of f(x) unchanged in amplitude but [(p)/2] out of phase with the strong light field at the zero spatial frequency position. The transfer function allows quantification of the resulting image spatial frequency distribution G(n) with that of the desired image field spectrum F(n) (See, for instance [16], [17,page 456], [18,page 315]).

Consider, however, what would happen if f(x) were not small. An appropriate description of the phase object would be forced to realise the truly non-linear nature of the exponential in equation 2.1 above and thus include higher powers of f(x) in the Taylor expansion so that

resulting in a spectrum described by
where * denotes a convolution. To define a transfer function in this case would be much more difficult, indeed this simple expansion leads straight to the fundamental problem with phase object spectra - that it is simply not possible to define a single transfer function for the spectrum of an arbitrary function f(x). From equation 2.5 above it will be seen that this is because the actual spectrum in the Fourier plane G(n) is a function of multiple convolutions of the desired spectrum F(n) and thus is unique to the particular function f(x). Without a spectrum transfer function it is no longer possible to specify the nature of the image spatial frequency distribution precisely.

The one-to-one mapping of spatial frequency with a unique position in the Fourier plane that occurs if f(x) were truly an amplitude transmittance object (or a very weak phase object) is no longer a valid association. In deeper analysis equation 2.5 shows the introduction of light at frequencies other than those specified by the desired spectrum F(n) so that a further optical transform will interpret this light as belonging to spatial frequencies present in f(x). It is reasonable, after performing a phase contrast operation, to expect an image with an intensity distribution which now only approximates f(x) in cases where f(x) is no longer `small'.

1.1  Fourier Series Complex Objects

It has been common practice to study aspects of the imaging of an optical system by using a Fourier Series object. The merit of this method is demonstrated in chapter one where the analysis due to Zernike [5] was presented as a framework within which to analyse various phase visualisation techniques. It was shown most clearly, for example, that the Schlieren image is akin to the spatial derivative of the optical path variations of the object, f(x). The ease with which this result is obtained highlights the fact that a series approach can give considerable insight to the processes involved.

Of considerable importance in applying this method is the work of Ichioka et al [19]. Starting with an expression derived by Hopkins [20] for the intensity distribution of an arbitrary object illuminated by a partially coherent light source, Ichioka performs a detailed study into the imaging of a single frequency sinusoidal complex object. The image intensity is expressed as a Fourier series, the coefficients of which are related to the degree of coherence of the illumination and the complex transmission of the object. Shortly after this work a follow up paper (Ichioka, Suzuki [21]) demonstrated that their analysis could be extended to cover the imaging of a general periodic complex object. By specifying the general properties of a trapezoidal-like periodic function (pulse height - `C', DC bias - `A', etc.), Fourier coefficients were determined as functions of these parameters. Thus the effect of the partially coherent optical system on a whole range of similar objects could be carried out by performing the analysis on just one Fourier series. Selective alteration of `C', `A' etc. in the resulting expression for the image intensity effectively changed the object transmission function.

Ichioka et al limited their analysis to cases where the phase nature of a complex object was identical in structure to the amplitude transmittance of the object (Figure 2.1). This situation models the common phenomenon known to occur in photographic transparencies, where a relief image of height proportional to the amplitude transmittance of the transparency is found.

In a partially coherent system it was found that an abrupt change in either amplitude or phase of a complex object is found to have a strong influence on the appearance of the intensity image of the object. Also, for the class of complex object described in the last paragraph, if the phase structure of the object is always small, phase changes become of lesser importance in changing the image structure as the amplitude contrast of the object increases. This result was found to be independent of the degree of coherence of the illumination and provides confirmation of the intuitive limiting case where the phase variations become negligible in comparison with amplitude changes, and we have almost a pure amplitude object.

1.2  Pure Phase Objects

Of more relevance to this chapter is the work carried out on the imaging of pure phase objects in a partially coherent optical system. Ichioka et al concluded that in near coherent or coherent optical systems (such as the 6-f optical information processing bench) sharp boundaries of phase in a phase object are highly influential in the image structure. In particular a sharp boundary appears as either a peak or dip of the local intensity, surrounded by characteristic `ringing' phenomena associated with coherent imaging of sharp boundaries. (It is assumed throughout that a [(p)/2] phase contrast filter is used to allow formation of a visible image of the phase object.) The small phase expansion of Zernike (equation 2.2), shows the spectrum of a weak phase object bears a close resemblance to an amplitude version of f(x). Using this fact as a link between the worlds of amplitude and phase objects, it is therefore plausible to suggest that phenomena occurring in the coherent imaging of amplitude objects, such as that described above, should be expected to occur in the phase-contrast method of imaging pure phase objects.

2  Phase Object Spectra

Ichioka et al came to their conclusions from the study of a range of similar objects. Fundamental to their calculations was the equation for image intensity. If the object space has a transmittance described by
where n₀ is the fundamental spatial frequency, a_j is the Fourier coefficient of spatial frequency jn₀ and x denotes position in the object space. It is shown [,page 922] that the image intensity may be expressed in the form
where x^¢ denotes the image space co-ordinate and T(j,k) is a function known as the transmission cross-coefficient, which accounts for the coherence of the illumination.

For the special type of object used in their modelling (object phase proportional to object amplitude transmittance for a rigorously defined object function) the a_j are known precisely. However, for the case of an arbitrarily shaped phase object the a_j are much more difficult to calculate. This shall be expanded upon in section 2.2.1 but for now it has hopefully been demonstrated that

1. The representation of a complex object by a Fourier series is a valid technique often yielding both insight into a problem and valuable quantitative results.
2. It is of practical importance, as in the example here, to pursue methods of calculating the Fourier coefficients a_j in equation 2.7 describing the complex transmission of the object space.

2.1  Spectrum Computation - Taylor Series

In order that a comparison be made between methods, it shall prove beneficial to illustrate the difficulties in calculating a phase object spectrum with the following example.

Recall that the `large object' Taylor expansion for a pure phase object is given by

and results in a spectrum of form
We can utilise the Fourier series approach in the following manner. Let the phase variation f(x) be represented as
where the C_m are complex. The function f(x) is required to be REAL as it describes the spatial nature of the phase retardance, so that
Then our Taylor expansion becomes
The Fourier transform of this expansion is then

Already the equation has become cumbersome by including only three powers of f(x) in the object space expansion. Let us calculate the spectrum of the simplest phase object allowable in this method

where `n_j' is an integer, at the zero'th, fundamental and first harmonic spatial frequencies. For this object we have

The summation for the zero spatial frequency becomes

and insertion of the specific value of the coefficients for this example gives

The summation for the fundamental spatial frequency at n = +n_j is

which becomes

In a similar fashion the series for the first harmonic is found to be

The main drawback of this method is the intensive calculation required to find even a few terms in the series describing the spectrum. From the equations above one may observe that if the spectrum G(n) is to include terms in ([(a_j)/2])^M one must expand the exponential series to include powers of f(x)^M, which results in `M' products of series in equation 2.13. This, in turn, requires the determination of M variable sets m, n, p etc. which obey the relation

where M=3 in this example.

Consider the more realistic example where C_m extend from -L to +L rather than having an infinite number of terms. To obtain a spectrum of highest power in [(a_j)/2] of `M' it is required to find all a_i obeying the relation

where

It is therefore necessary to search through (2L)^M combinations to include the M'th power of [(a_j)/2]. Furthermore, for the single frequency case taken in the example above, it appears that the most significant term in the spectrum series for the K'th harmonic frequency has a lowest power of ([(a_j)/2])^K+1. Thus

1. An enormous increase in computation is required to increase the accuracy of the series representation of the spectrum at any spatial frequency n.
2. An equivalent increase in computation is required to find even the first term of the series as n increases.

The function (2L)^M is graphed in figure 2.2 for the case of a phase object described by only sixteen Fourier coefficients.

To recap on what has been learned from this example, the calculation using a single spatial frequency shows

1. The generation of light at frequencies other than at the fundamental.
2. The most significant term in the series representation of the spectrum at the K'th harmonic is proportionate to ([(a_j)/2])^K+1. The fundamental frequency light also obeys this relation if the value K=0 is assigned to it.
3. For small phase objects, (a_j small) the results are in agreement with a Taylor expansion including only the first power of [(a_j)/2]. For a simple sinusoidal phase object, if the Taylor expansion accurately represents the object space complex transmission by including powers of [(a_j)/2] up to K, say, the spectral orders at frequencies higher than n = Kn_j are negligible. Here, `n_j' is the fundamental spatial frequency.

3   The Bessel Function Method

It has been shown that much complexity in spectrum calculation arises from the fundamental convolution nature of the problem. In this section a straightforward expansion of the complex object transmission is introduced allowing a cleaner approach to spectrum calculation. The results of the previous section are incorporated into an alternative procedure of greatly reduced complexity.

In place of the Taylor expansion, the Jacobi-Anger expansion [22] will be shown to be of fundamentally greater use. The expansion may be written

where J_m(a_j) is the m'th order Bessel function of argument a_j. For the purpose of this thesis the argument of the exponential represents a spatial frequency component of the object phase retardance. Thus the above equation can be written
Upon Fourier transformation we immediately have the complete spectrum of a single spatial frequency in the phase object without resort to approximation. The frequency spectrum is given by
Each sinusoid of the phase retardance f(x), equation 2.10, results in an infinite spectrum of equally spaced d-functions, this array being termed a Dirac comb or more frequently just a comb [17, pages 60-63]. The separation of d-functions on a comb is proportional to the spatial frequency which generated that comb. Note that the m'th d-function has an amplitude described by the m'th order Bessel function, so that the spectral order of the comb is synonymous with the Bessel function order describing the amplitude there. For the m=±1 comb orders, observe that although the amplitude is changed by J_±1(a_j) the phase information of these orders is identical to that which would describe an pure amplitude object of transmission
This is a most significant observation since the phase transfer function of a spectrum (rather than the amplitude transfer function) plays the greater part in determination of the image characteristics [23], [24]. It will be shown that these orders are primarily responsible for the image intensity being a linear function of object phase in the weak phase approximation.

The form of J_m(a_j) for m up to 3 is shown in figure 2.3. One representation [22] of the Bessel function set is by the infinite power series

The power series representations for J₀, J₁ and J₂ are
so that the lowest power of [(a_j)/2] in J_m(a_j) is ([(a_j)/2])^m.

Comparison of the equations for the spectrum derived from the Taylor expansion in the last section for the zero spatial frequency, fundamental and first harmonic terms (equations 2.18, 2.20 & 2.21) with the above expressions reveals that, apart from a complex premultiplier of `i', they are identical as expected. In effect, the Taylor series expansion of the spectrum contains terms describing the amplitude at each comb location ( J_m(a_j) ), and the information on how the combs belonging to each spatial frequency are convolved with one another. The principle advantage of the Bessel series expansion is that it may be thought of as incorporating a great deal of convolution under the banner of just one function, J_m(a_j). A phase object comprising of many spatial frequencies can be described by

where the fundamental spatial frequency n_j of f(x) is taken as unity for simplicity. The Fourier Transform of a product is the convolution of the spectra of each factor [14,page 108], so that the problem of calculating the phase object spectrum then reduces to one of finding a method of convolving many d-function combs.

3.1  Terminology

It will prove convenient to introduce some terminology which greatly simplifies lengthy discussion. It is assumed in the following discussion that the object Fourier components have amplitudes a_j lying in a range bounded by a value a_max such that the first order Bessel function just begins to show a non-linear response. From figure 2.4 one might assign a_max to be 0.5 as a reasonable limit.

For any particular comb the first d-functions either side of the comb origin are denoted in this thesis as the primary comb orders ( m=±1 in equation 2.27). These orders have an amplitude proportional to J_± 1(a_j) where a_j is the object spatial frequency responsible for that comb. The phase of the j'th object spatial frequency F_j is unaltered at the primary comb order also. Those spectral orders described by the higher order Bessel functions ( m ³ 2 ) are non-linear functions of both a_j and F_j. These orders are, among other things, centres of convolution for every other comb in the expression for the spectrum of an N spatial frequency phase object. Due to the fact that their presence is almost unnoticed for small a_j though they are generally undesirable for larger a_j, these orders shall be termed ghost orders. Figure 2.5 illustrates the terminology introduced here.

Only at the primary orders of each comb do we have a one-to-one mapping of the spatial frequencies of f(x) with position in the frequency plane. More will be said on this towards the end of section 2.4.1 but for now note that this statement has some additional qualifications.

1. The j`th spatial frequency component of f(x) with amplitude a_j is represented in the frequency plane by the approximately linear function J₁(a_j) at the positions n = ±j.
2. Once convolution takes place with the combs from every other object spatial frequency, the light field at n = ±j has additional terms to J₁(a_j). Convolution effects thus serve to corrupt the one-to-one mapping partially present in each comb.

4  Convolution - a mathematical framework

Having introduced the idea of Bessel function combs in the last section, this section will specify a simple technique whereby the spectrum of a phase object described by N Fourier components may readily be determined.

4.1  Two Frequency Case

From equation 2.33, an N-frequency phase object may be described by
Replacing each exponential factor with the corresponding Bessel function comb and performing the convolution leads to a spectrum of form
It might seem that nothing has been gained over the Taylor expansion method if complexity were the deciding factor. To calculate the spectrum at any frequency one must find all m_i satisfying
which is one linear equation in N unknowns. An alternative formulation of the problem is clearly required.

The purpose of this section is to introduce just such a method, whereby the spectrum is built up in stages. The example of a two frequency phase object is chosen to illustrate this new technique both for the simplicity with which the terminology is introduced and because it will prove to be of central importance. In the process of convolution, one specific frequency location n = g is chosen as a point of observation at which to calculate the spectrum. A suitable label with which to refer to any particular comb is the d-function spacing of that comb, so that in this analysis, the spectral comb resulting from the a₁ Fourier coefficient is denoted the unit comb and the comb resulting from the a₂ coefficient is denoted as the 2-comb.

The summation indices m_i in equation 2.35 denote the d-functions of each comb, as illustrated in figure 2.5, where each arrow represents a d-function location. In the process of convolution, the 2-comb origin sits on each d-function of the unit comb in turn and the complex amplitude at n = g is built up term by term. The summation index m₁ thus denotes which d-function of the unit comb the 2-comb is centred upon, so that m_N denotes which d-function of the 2-comb then sits at the point of observation g.

In the example shown in figure 2.6, g = +1 and the 2-comb sits on top of the negative primary order of the unit comb. The positive primary order of the 2-comb then sits on on point of observation. This particular configuration is therefore described as having m₁=-1 and m₂=+1. Table 2.1 lists several comb configurations also resulting in an additional contribution to the complex amplitude at g = +1.