1

Computing the 2-D Discrete Fourier Transform or Sweet-talking MATLAB into Making Cool Pictures Marshall. C. Overcast, MOIL MSU, and Paul. W. Nugent, ORSL MSU

Abstract—This paper covers the methods used to produce accurate and decipherable two dimensional Fourier Transforms using the MATLAB software package. This software implements a two dimensional Fast Fourier Transform (FFT) in its fft2() script. This paper covers the pertinent assumptions that are made by this algorithm. Once the 2-D Fourier Transform of an image is computed, the methods of computing the Point Spread Function and Optical Transfer Function are also discussed. The guidelines presented in this paper assume that the reader is familiar with the integral form of both the one dimensional and two dimensional Fourier transforms, and familiarity with viewing the frequency space outputs of the transforms.

I. INTRODUCTION

T

He Fourier transform was initially based on the concepts of Jean Baptiste Joseph Fourier who in 1822 published the paper “Théorie analytique de la chaleur “. Within this paper he develops the idea that all discontinuous functions could be expressed as a series of sinusoidal waves. His formulation was flawed; however, this idea was a breakthrough and lead to the development of the Fourier series. The Fourier transform of today is capable of expressing a signal as a composition of sinusoidal function. With the development of the computer during the mid 20th century there came a need for a quick method of determining the discrete Fourier transform of a signal. Initially these methods were very processor and memory intense. With the computer technology of the mid 20th century, this meant very large computer and lots of time (24 hrs for a 128 x 128 image!). To decrease the time it took to process these discrete transforms many versions of the Fast Fourier Transform (FFT) were created. The mostly widely used algorithm was developed in 1965 by J.W. Cooley and John Tukey name the Cooley-Tukey algorithm. This algorithm reduced the FFT computation time to O(nlog(n)) from a O(n2). The version of the FFT implemented in MATLAB is largely based on the M. C. Overcast is with Micro Optics and Imaging Lab at MSU, Bozeman, MT 59715 USA (e-mail: [email protected]). P. W. Nugent is with the Optical Remote Sensing Lab at MSU, Bozeman, MT 59715 USA (e-mail: [email protected]).

Cooley-Tukey algorithm with other optimizations. “The execution time for fft depends on the length of the transform. It is fastest for powers of two. It is almost as fast for lengths that have only small prime factors. It is typically several times slower for lengths that are prime or which have large prime factors.”[1]

II. IMPLEMENTATION OF THE FFT2 IN MATLAB The MATLAB software environment provides a powerful set of tools for performing signal analysis and other applications. This tool set comes with a learning curve that can lead to some non-intuitive results when one first attempts an operation. This is very much the case for the implementation of the fft2 MATLAB function. This following section provides a guideline of the basic operation of the function and shows how we sweet-talk the function into producing the outputs we desire. The first of these problems arises with the special layout of the input and the output of these functions. Many of us have become familiar with the frequency space results of the one dimensional Fourier transform. We expect that from most signals there will be a relatively large ‘DC’ component located at the origin and the transform will be mirrored about the origin. This holds true for the two dimensional case where we find the ‘DC’ component located at the origin of the two dimensional plane and the frequency space function is mirrored about this point. This however is not the result of directly using the fft2 function without any pre or post processing. First, we need to create a simple function to be transformed. For this example, a simple 2-D rect function will be used. In order to place our rect in exactly in the center of our matrix we need to define both the rect and the matrix with odd dimensions. We can define our rect using the following MATLAB commands.

2

rect = zeros(51); rect(25:35,25:35)=ones(11); surf(rect)

input and the output of the fft2 function. The following section shows how this can be done for the previously defined rect function. First the rect will be shifted in image space. rectshifted=fftshift(rect); surf(rect)

Figure 1: The mesh of a two dimensional rect function 11x11 units wide. This function uses an odd number for its side and its image to ensure that it is exactly centered. Now the fft2 command can applied to this rect. rect_fft = fft2(rect); These complex functions are commonly viewed as magnitude plots. To view the frequency space function the following command should be used. surf(abs(rect_fft))

Figure 3: The rect function after it has been shifted using the fftshift command. If this image was replicated in both the x and y axis, it would still remain continuous at the boundaries. Now that we have shifted the image, we apply the two dimensional Fourier Transform to it. rectfft=fftshift(rectshifted); surf(rectfft)

Figure 2: The direct fft2 of the rect function in Figure 1. The DC components of the frequency space plot have been placed in the corners of the matrix. This result is not what has come to be expected from experience with Fourier Transforms. We expect that the low frequency components should be located near the center of the output matrix and not the edges. MATLAB supplies functions that both pre and post process the image both in image space and frequency space so that the expected results are obtained. The MATLAB function fftshift can be used to shift both the

Figure 4: The rect function after being Fourier Transformed. Note that the DC components are still in the four corners of the image.

3

Now after a post processing shift, again using the fftshift command, the expected 2-D sinc pattern can be seen. shiftedrectfft=fftshift(rectfft); surf(shiftedrectfft)

It may be desired, by Tibetan monks post Doctorial Mathematicians, to view the phase of the transformed image. If so desired this may be viewed using the following commands. phase=angle(atan(imag(shiftedrectfft)/real(shiftedrectfft))); surf(phase)

Figure 5: After applying the fftshift command to the Fourier Transform, the image looks a bit more familiar.

Figure 7: This is the phase of the rect function. Note this information goes through many 2π phase shifts. These shifts give the phase a stair step pattern where the phase for a rect should be nearly linear [3].

Now it may be desired at times to see not only the magnitude of the Fourier Transform, but also the finer detail that arises from the power spectrum of the image. To do this, the power spectrum of the frequencies can be displayed using the following commands.

So now that the pre and post processing of an image and the information delivered by the fft2 function is understood, several common 2-D functions and their Fourier Transforms can be explored. The first function is a cosinusoidal wave defined by the following MATLAB commands.

power=10*log(abs(shiftedrectfft).^2+1); surf(power)

2dcos = cos(2πfxx + 2πfyy);

(Note the +1 is to remove errors due to a log(0) possibly occuring)

Figure 6: This is the power spectrum of the rect function plotted in DB instead of linearly. From this the finer structure of the transform is visible.

Figure 8: A cosinusoidal wave propagating at 45 degrees. For this propagation angle the values fx and fy in the function are equal.

4

The Fourier Transform of 2dcos defined above is shown below.

After crunching numbers for a few seconds, MATLAB returns the following Fourier Transform.

Figure 9: The Fourier Transform of the cosinusodial wave. Note that there are 2 delta functions representing the frequency of the cosinusodial function. Also note that these two delta functions are turned 45 degrees from axes.

Figure 11: The Fourier Transform of the circ function shown above. Notice the non-airiness of the transform at the edges.

Another common function is the circular aperture that is commonly found in optics labs. By exploring its Fourier Transform some insight into its diffraction pattern can be gained. The circ used below is defined by using the fspecial(‘disk’,radius) command. This command softens the edges of the function, alleviating the discrete nature the circular would otherwise demonstrate. Discretely defined circles tend to have edges that are jagged. These jagged edges introduce frequency components that would not be present in a real circular aperture. These frequencies act as noise in the pattern and corrupt the frequency space diagram. To avoid this problem the ‘soft’ edged circ should be used.

Now, this does not look very much like an Airy pattern. Portions in the center appear to exhibit Airy like qualities, but the tops and edges have some other junk going on. This is due to the assumptions made by the Cooley-Tukey algorithm to decrease the processing time. The algorithm assumes that the image is repeating in both the X and Y directions. To bring about a more interpretable result this assumption must be avoided. The easiest solution is to simply pad the image with zeros increasing the period of the repeating pattern. The circ shown below is the same size as the previous one, but has had zeros padded around it increasing the image size by a factor of four.

Figure 10: A simple circ function which is representative of a circular aperature. This function is designed with "soft" edges using the fspecial('disk',radius) command.

Figure 12: This circ is the same size as the circ in Figure 10. The image has been padded with 0's to increase the period of the repeating pattern assumed in the MATLAB fft2 algorithms.

5

Now the Fourier Transform of the padded circ can be taken as before.

Figure 13: The Fourier Transform of the circ with 0 padding viewed over the same region in frequency space. Note that compared with the other fourier transoform that this has much more of an Airy pattern to it. The frequency space image now resembles an airy pattern. The quality of this Airy pattern has not been determined, nor the real extent of the zero padding necessary to accomplish an Airy pattern of given error. These tests would be interesting though none the less and were not within the scope of this project due to time constraints. With the zero padding technique other circularly symmetric images and their Fourier Transform can now be explored. The image below is a circularly symmetric cosinusoidal wave. We have defined this wave by 2

2 0.5

r = cos((2πf (x + y ) )

Now taking the Fourier Transform of this wave is shown below.

Figure 15: The Fourier Transform of the cosinusoidal wave is shown above. Note that this resembles a delta function repeated at a given radius distributed through all angles. III. THE POINT SPREAD FUNCTION AND THE OPTICAL TRANSFER FUNCTION Now that the 2-D Fourier Transform of a given image can be computed with reasonable accuracy this transform can be used to compute other properties of the optical system. The point spread function (PSF) is one such property. The PSF defines the propagation of electromagnetic radiation from a point source. It is defined in spherical coordinates for a Lambert-type radiator[2]. To compute the PSF of a square aperture, we use the following MATLAB commands. sqft=(fftshift(fft2(fftshift(sq1)))); psf = sqft.*conj(sqft); surf(psf);

Figure 14: A cosinusoidal wave with zero padding is defined above. This could represent a circularly symmetric grating material or something else. Figure 16: The PSF of the rect function shown in Figure 1. Another useful optical property of a system that is related to the PSF is the Optical Transfer Function (OTF). The Optical transfer function represents image contrast to object contrast

6

ratio of an optical system when viewed in frequency space, taking into account the phase shift between positions occupied by the actual and ideal image [2]. The OTF can be calculated by taking the inverse Fourier Transform of the PSF. To calculate the OTF of the previously defined rectangular aperture the following MATLAB commands can be used.

VI. REFERENCES [1] [2] [3] [4]

otf = fftshift(ifft2(fftshift(psf))); surf(otf);

[5]

[6]

[7] [8]

Figure 17: The OTF of the rect function defined in Figure 1. This closely resembles the convolution of the rect function with itself. IV. CONCLUSIONS We have demonstrated in this paper that an accurate 2-D Fourier Transform of a given image can be determined using MATLAB’s fft2 command. Pre and Post processing of the image matrix through the use of fftshift makes the interpretation of the fft2 results easier. In the case of circularly symmetric functions, zero padding may be required bring about the desired results. This is due to the fft2’s assumption that the image is infinitely periodic in both x and y. From these Fourier Transforms it is trivial to calculate the Point Spread Function (PSF) & Optical Transfer Function (OTF) of an image as well. The interpretation of the PSF and the OTF is left for a later more in-depth discussion. For now it is know that these functions are useful tools in analyzing an optical system.

V. ACKNOWLEDGMENTS The authors gratefully acknowledge the contributions of Dr. J. A. Shaw for his help in troubleshooting these MATLAB functions and aid in interpreting their results.

MATLAB Help Files “MATLAB Function Reference fft”, MATLAB version 7.0.1, Matworks, September, 2004 Wilson, G. 1995. Fourier Series and Optical Transform Techniques in Contemporary Optics. New York: John Wiley & Sons, Inc. Shaw, Joseph. Inquisition by authors, 25 Jan 2006, Bozeman, MT. Mentally recorded. Mathworks. (2006). MATLAB (Version 7.0.4) [Computer software]. Natick: Mathworks. Wikipedia contributors, “Joesph Fourier" Wikipedia: The Free Encyclopedia, http://en.wikipedia.org/wiki/Joseph_Fourier (accessed January 24, 2006) Wikipedia contributors, Cooley-Tukey FFT Algorithm," Wikipedia: The Free Encyclopedia, http://en.wikipedia.org/wiki/CooleyTukey_FFT_algorithm (accessed January 24, 2006) Goodman, J. 2005. Fourier Optics. Englewood: Roberts & Company Phillips, C. & Parr, J. & Riskin, E. 2003. Signals, Systems, and Transforms. Upper Saddle River: Prentice Hal

Computing the 2-D Discrete Fourier Transform or Sweet ...

MT 59715 USA (e-mail: [email protected]). Cooley-Tukey algorithm with ... The first of these problems arises with the special layout of the input and the ...

2MB Sizes 0 Downloads 198 Views

Recommend Documents

3D Fourier based discrete Radon transform
Mar 21, 2002 - School of Computer Science, Tel Aviv University, Tel Aviv 69978, Israel ... where x = [x,y,z]T, α = [sinθ cosφ,sinθ sinφ,cosθ]T, and δ is Dirac's ...

3D discrete X-ray transform
The continuous X-ray transform of a 3D function f(x,y,z), denoted by Pf , is defined by the set of all line integrals of f . For a line L, defined by a unit vector θ and a ...

Fourier transform solutions to Maxwell's equation.
The aim of this set of notes is to explore the same ideas to the forced wave equations for the four vector potentials of the Lorentz gauge Maxwell equation.

Warped Discrete Cosine Transform Cepstrum
MFCC filter bank is composed of triangular filters spaced on a logarithmic scale. ..... Region: North Midland) to form the train and test set for our experiments.

3D discrete X-ray transform - ScienceDirect.com
Available online 5 August 2004. Communicated by the Editors. Abstract. The analysis of 3D discrete volumetric data becomes increasingly important as ... 3D analysis and visualization applications are expected to be especially relevant in ...

Uniform Discrete Curvelet Transform
The effective support of these curvelet functions are highly anisotropic at fine scale, ... The set of windows is then used in Section VI to define a multi-resolution ...

Fractional Fourier Transform Based Auditory Feature for ...
where a is the normalization constant, b is a bandwidth related parameter, fc is ..... evaluation. [Online]. Available: http://www. nist.gov/speech/tests/lang/index.htm.

Broad-Band Two-Dimensional Fourier Transform Ion ...
w, and w2 that correspond to the mass-to-charge ratios of the ions that participate as reactants and products, respectively. ... The first rf pulse PI in the sequence of eq ..... was incremented in 120 steps of 1 ps, with 1K data points recorded in f

Fast Fourier Transform Based Numerical Methods for ...
Analyzing contact stress is of a significant importance to the design of mechanical ...... A virtual ground rough surface can be formed through periodically extend-.

Image Compression Using the Discrete Cosine Transform
NASA Ames Research Center. Abstract. The discrete ... The discrete cosine transform of a list of n real numbers s(x), x = 0, ..., n-1, is the list of length n given by:.

Image Compression and the Discrete Cosine Transform ...
We are now ready to perform the Discrete Cosine Transform, which is accomplished by matrix multiplication. D : TMT r. 5. In Equation (5) matrix M is first multiplied on the left by the DCT matrix T from the previous section; this transforms the rows.

THE FOURIER-STIELTJES AND FOURIER ALGEBRAS ...
while Cc(X) is the space of functions in C(X) with compact support. The space of complex, bounded, regular Borel measures on X is denoted by. M(X).