In the last article we saw that the Point Spread Function and the Modulation Transfer Function of a lens could be easily obtained numerically by applying Discrete Fourier Transforms to its generalized exit pupil function twice in sequence.^{[1]}

Obtaining the 2D DFTs is easy: simply feed MxN numbers representing the two dimensional complex image of the pupil function in its space to a fast fourier transform routine and, presto, it produces MxN numbers that represent the amplitude of the PSF on the sensing plane. Figure 1a shows a simple case where pupil function is a uniform disk representing the circular aperture of a perfect lens with MxN = 1024×1024. Figure 1b is the resulting PSF.

Simple *and* fast. Wonderful. Below is a slice through the center, the 513th row, zoomed in. Hmm…. What are the physical units on the axes of displayed data produced by the DFT?

Less easy – and the subject of this article as seen from a photographic perspective.

#### A Numerical Perspective

For this discussion we take a numerical view of the simplified diagram representing the exit pupil and sensing plane (Figure 3 below, not drawn to scale). Thus the exit pupil function refers to an array of MxN numbers indexed by integers u and v which range from 0 to M-1 and 0 to N-1 respectively. Similarly for the relative Point Spread Function on the sensing plane and Modulation Transfer Function . The only physical parameters on the diagram are shown in red: aperture diameter , effective focal length and sampling pitch and (not shown and which would however be the same as and resp. in a square-pixeled sensor).

#### The Physical Units of the Exit Pupil Function

To get the axes of the PSF we first have to determine the physical units of the exit pupil function, which fortunately we do have. Using the same example as in the last article, with the diameter of the disk at 5mm, we can easily put axes on Figure 1a above since I counted that exactly 9 diameters fit within a width and a height of the exit pupil image^{[2]}. That means that the pupil plane must correspond to a 45mmx45mm space.

If the plane is sampled times in the horizontal direction, with = N = 1024 per Figure 1, that means that linear physical sample spacing is 45mm/1024 = 0.0439mm on the physical pupil plane.

Denoting the number of times the diameter of the aperture fits horizontally within the pupil plane as , the horizontal axis of the pupil function in Figure 1a at index position would then correspond to physical value

(1)

and the spacing between samples to , all in the same physical units as since and are dimensionless. in this case represents integers 0,1,2… to n-1.

That’s with the origin top left, as shown in Figure 1. The convention is however to put the origin of the pupil function in the center of the image as shown in Figure 3. If we shift the image halfway down and to the right^{[1]} we need to let the units change accordingly by modifying and . Shifted then becomes

Of course remains unaffected by shifting. Assuming square pixels, is the same as and physical values for can be obtained similarly (in fact the axes are the same in the case of a square pixel and a square image, as in this example).

#### From Pupil Function to PSF Units

As described in the previous article, to generate the Point Spread Function on the sensing plane we take the Fourier Transform of the Pupil Function after performing coordinate substitution

(2)

with the wavelength of light and the effective focal length of the lens (the distance between parallel and planes in Figure 3 above).

Since our data is digitized with samples at regular intervals we can use the Discrete Fourier Transform to perform the conversion numerically. The DFT produces an output data set of the same size as the input MxN data set and at a physical spacing equal to the inverse of the input sample spacing. One over a spacing is sometimes referred to as a ‘frequency’. We can think of it like that if it helps us follow the math – though in this case the result of the Fourier Transform leads to a PSF, clearly a quantity in the spatial domain.

The sequence of values output by the DFT corresponds to an integer multiple of this frequency, from 0 to n-1 times, by definition. Here is the Discrete Fourier Transform expressed with the variables at hand^{[3]}

(3)

in this case represents integers 0,1,2… to N-1 and M-1 for and respectively. So with physical units of the horizontal axis of the PSF in Figure 1b at integer index position would correspond to physical value

To convert to units of on the sensor we simply multiply by according to (2) so that for instance the fourth pixel on the sensing plane now represents 3 in physical units.

Substituting effective f-number for in that fraction we obtain the following physical quantity at integer index position on the horizontal axis of the sensing plane

(4)

and , all in the same physical units as because and are dimensionless. That’s with the origin top left. If we want the PSF to be centered we shift the origin half an image both ways as before:

With the origin shifted to the middle of the image, the 513th sample in Figure 2 becomes the origin. It is then easy to see from its central slice that the Airy pattern hits first zero at x = 11, corresponding to 11/9*0.5*10 = 6.11 um in physical units (from the last article , um and ). Correct, because theory predicts that the Airy’s first ring radius should be 1.22 .

#### PSF to MTF Units

To get Modulation Transfer Function from Point Spread Function we take a second Discrete Fourier Transform of the data, as shown by the third line of code up top. This time there is no coordinate substitution so the physical units corresponding to the output of the DFT are the bona fide linear spatial frequencies present in the PSF. Since the physical units corresponding to integer index on the horizontal axis of the relative MTF are

(5)

and , all in cycles per the same units as . Sometimes it is useful to shift the origin of the frequency diagram to the center of the image. If so proceed as before to also shift .

#### Bringing Them All Together

In Octave or Matlab, the three lines of code up top are therefore complemented by the three lines of code below to generate axes for the three planes discussed in this post:

Plots of the pupil, image and frequency spectrum can therefore now be shown with axes in physical units:

and in the vertical direction can be obtained using the same procedure starting from .

#### On Q and N

You may have noticed that the only variables required to determine the physical units of exit pupil, PSF and MTF are , , and . A Q value of about one order of magnitude is mainly necessary in order to obtain relatively unaliased data from the discrete fast fourier transform routine but for continuous photographic thought experiments it can be considered to be equal to one. In that case ** the physical units of PSFs and MTFs are seen to depend only on a single figure, the product **. This realization can provide some neat insights into the basic properties of light on the sensing plane after it has gone through a lens.

#### Notes and References

_{1. The substance of the first line is represented by the Fast Fourier Transform of the Pupil Function. The shifts are there to ensure that the resulting PSF is centered. Thanks to Alan Robinson for this particular shifting sequence. The second line is a better and faster way to compute the magnitude of the intensity spectrum squared, thus producing the PSF. The third line takes the magnitude of the Discrete Fourier Transform of the PSF, which produces the MTF.
2. Simplified Figure 3 is not drawn to scale.
3. The definition of the Discrete Fourier Transform can be found here.
}