Главная
Computational Fourier Optics : a MATLAB tutorial (SPIE Tutorial Texts Vol. TT89)
Computational Fourier Optics : a MATLAB tutorial (SPIE Tutorial Texts Vol. TT89)
David George Voelz
0 /
0
Насколько вам понравилась эта книга?
Какого качества скаченный файл?
Скачайте книгу, чтобы оценить ее качество
Какого качества скаченные файлы?
Computational Fourier Optics is a text that shows the reader in a tutorial form how to implement Fourier optical theory and analytic methods on the computer. A primary objective is to give students of Fourier optics the capability of programming their own basic wave optic beam propagations and imaging simulations. The book will also be of interest to professional engineers and physicists learning Fourier optics simulation techniqueseither as a selfstudy text or a text for a short course. For more advanced study, the latter chapters and appendices provide methods and examples for modeling beams and pupil functions with more complicated structure, aberrations, and partial coherence. For a student in a course on Fourier optics, this book is a concise, accessible, and practical companion to any of several excellent textbooks on Fourier optical theory.
Категории:
Год:
2010
Издательство:
SPIE Press
Язык:
english
Страницы:
249
ISBN 10:
0819482048
ISBN 13:
9780819482044
Серия:
SPIE Tutorial Texts TT89
Файл:
PDF, 4,96 MB
Ваши теги:
Скачать (pdf, 4,96 MB)
 Открыть в браузере
 Checking other formats...
 Конвертировать в EPUB
 Конвертировать в FB2
 Конвертировать в MOBI
 Конвертировать в TXT
 Конвертировать в RTF
 Конвертированный файл может отличаться от оригинала. При возможности лучше скачивать файл в оригинальном формате.
 Пожалуйста, сначала войдите в свой аккаунт

Нужна помощь? Пожалуйста, ознакомьтесь с инструкцией как отправить книгу на Kindle
В течение 15 минут файл будет доставлен на Ваш email.
В течение 15 минут файл будет доставлен на Ваш kindle.
Примечание: вам необходимо верифицировать каждую книгу, которую вы отправляете на Kindle. Проверьте свой почтовый ящик на наличие письма с подтверждением от Amazon Kindle Support.
Примечание: вам необходимо верифицировать каждую книгу, которую вы отправляете на Kindle. Проверьте свой почтовый ящик на наличие письма с подтверждением от Amazon Kindle Support.
Возможно Вас заинтересует Powered by Rec2Me
Ключевые слова
fourier^{179}
irradiance^{171}
propagation^{168}
simulation^{167}
fig^{161}
optical^{154}
sampling^{147}
exp^{146}
diffraction^{135}
imaging^{133}
axis^{130}
analytic^{121}
lens^{117}
fft^{115}
sample^{112}
transform^{112}
frequency^{111}
coherence^{110}
array^{110}
optics^{108}
coherent^{108}
plot^{107}
rect^{98}
fraunhofer^{96}
spatial^{95}
pupil^{95}
fresnel^{94}
magnitude^{91}
wavefront^{90}
samples^{88}
psf^{82}
sinc^{80}
lambda^{79}
beam^{77}
matlab^{75}
sampled^{74}
bandwidth^{69}
periodic^{66}
fourier transform^{65}
convolution^{63}
vector^{60}
coordinates^{59}
fourier optics^{57}
width^{57}
chirp^{57}
wavelength^{53}
discrete^{53}
sample interval^{50}
incoherent^{49}
temporal^{49}
grating^{49}
impulse response^{48}
spectral^{48}
cos^{47}
script^{47}
aperture^{47}
propagator^{47}
criterion^{45}
aberration^{44}
xlabel^{44}
Связанные Подборки
0 comments
Вы можете оставить отзыв о книге и поделиться своим опытом. Другим читателям будет интересно узнать ваше мнение о прочитанных книгах. Независимо от того, пришлась ли вам книга по душе или нет, если вы честно и подробно расскажете об этом, люди смогут найти для себя новые книги, которые их заинтересуют.
1

2

Tutorial Texts Series Modeling the Imaging Chain of Digital Cameras, Robert D. Fiete, Vol. TT92 Cells Illuminated: In Vivo Optical Imaging, Lubov Brovko, Vol. TT91 Polarization of Light with Applications in Optical Fibers, Arun Kumar, Ajoy Ghatak, Vol. TT90 Computational Fourier Optics: A MATLAB® Tutorial, David Voelz, Vol. TT89 Optical Design of Microscopes, George Seward, Vol. TT88 Analysis and Evaluation of Sampled Imaging Systems, Richard H. Vollmerhausen, Donald A. Reago, Ronald Driggers, Vol. TT87 Nanotechnology: A Crash Course, Raúl J. MartinPalma and Akhlesh Lakhtakia, Vol. TT86 Direct Detection LADAR Systems, Richard Richmond, Stephen Cain, Vol. TT85 Optical Design: Applying the Fundamentals, Max J. Riedl, Vol. TT84 Infrared Optics and Zoom Lenses, Second Edition, Allen Mann, Vol. TT83 Optical Engineering Fundamentals, Second Edition, Bruce H. Walker, Vol. TT82 Fundamentals of Polarimetric Remote Sensing, John Schott, Vol. TT81 The Design of Plastic Optical Systems, Michael P. Schaub, Vol. TT80 Fundamentals of Photonics, Chandra Roychoudhuri, Vol. TT79 Radiation Thermometry: Fundamentals and Applications in the Petrochemical Industry, Peter Saunders, Vol. TT78 Matrix Methods for Optical Layout, Gerhard Kloos, Vol. TT77 Fundamentals of Infrared Detector Materials, Michael A. Kinch, Vol. TT76 Practical Applications of Infrared Thermal Sensing and Imaging Equipment, Third Edition, Herbert Kaplan, Vol. TT75 Bioluminescence for Food and Environmental Microbiological Safety, Lubov Brovko, Vol. TT74 Introduction to Image Stabilization, Scott W. Teare, Sergio R. Restaino, Vol. TT73 Logicbased Nonlinear Image Processing, Stephen Marshall, Vol. TT72 The Physics and Engineering of Solid State Lasers, Yehoshua Kalisky, Vol. TT71 Thermal Infrared Characterization of Ground Targets and Backgrounds, Second Edition, Pieter A. Jacobs, Vol. TT70 Introduction to Confocal Fluorescence Microscopy, Michiel Müller, Vol. TT69 Artificial Neural Networks: An Introduction; , Kevin L. Priddy and Paul E. Keller, Vol. TT68 Basics of Code Division Multiple Access (CDMA), Raghuveer Rao and Sohail Dianat, Vol. TT67 Optical Imaging in Projection Microlithography, Alfred KwokKit Wong, Vol. TT66 Metrics for HighQuality Specular Surfaces, Lionel R. Baker, Vol. TT65 Field Mathematics for Electromagnetics, Photonics, and Materials Science, Bernard Maxum, Vol. TT64 HighFidelity Medical Imaging Displays, Aldo Badano, Michael J. Flynn, and Jerzy Kanicki, Vol. TT63 Diffractive Optics–Design, Fabrication, and Test, Donald C. O’Shea, Thomas J. Suleski, Alan D. Kathman, and Dennis W. Prather, Vol. TT62 FourierTransform Spectroscopy Instrumentation Engineering, Vidi Saptari, Vol. TT61 The Power and EnergyHandling Capability of Optical Materials, Components, and Systems, Roger M. Wood, Vol. TT60 Handson Morphological Image Processing, Edward R. Dougherty, Roberto A. Lotufo, Vol. TT59 Integrated Optomechanical Analysis, Keith B. Doyle, Victor L. Genberg, Gregory J. Michels, Vol. TT58 ThinFilm Design: Modulated Thickness and Other Stopband Design Methods, Bruce Perilloux, Vol. TT57 Optische Grundlagen für Infrarotsysteme, Max J. Riedl, Vol. TT56 An Engineering Introduction to Biotechnology, J. Patrick Fitch, Vol. TT55 Image Performance in CRT Displays, Kenneth Compton, Vol. TT54 Introduction to Laser DiodePumped Solid State Lasers, Richard Scheps, Vol. TT53 (For a complete list of Tutorial Texts, see http://spie.org/x651.xml.) Tutorial Texts in Optical Engineering Volume TT89 PRESS Bellingham, Washington USA Voelz, David George, 1959Computational fourier optics : a MATLAB tutorial / David G. Voelz. p. cm. Includes bibliographical references and index. ISBN 9780819482044 1. Image processingDigital techniquesMathematics. 2. Fourier transform optics. 3. MATLAB. I. Title. TA1637.V674 2010 621.3601'515723dc22 2010050505 Published by SPIE P.O. Box 10 Bellingham, Washington 982270010 USA Phone: +1 360.676.3290 Fax: +1 360.647.1445 Email: Books@spie.org Web: http://spie.org Copyright © 2011 Society of PhotoOptical Instrumentation Engineers (SPIE) All rights reserved. No part of this publication may be reproduced or distributed in any form or by any means without written permission of the publisher. The content of this book reflects the work and thoughts of the author(s). Every effort has been made to publish reliable and accurate information herein, but the publisher is not responsible for the validity of the information or for any outcomes resulting from reliance thereon. Printed in the United States of America. To my dad, my best friend who always wanted to write a book. Introduction to the Series Since its inception in 1989, the Tutorial Texts (TT) series has grown to cover many diverse fields of science and engineering. The initial idea for the series was to make material presented in SPIE short courses available to those who could not attend and to provide a reference text for those who could. Thus, many of the texts in this series are generated by augmenting course notes with descriptive text that further illuminates the subject. In this way, the TT becomes an excellent standalone reference that finds a much wider audience than only short course attendees. Tutorial Texts have grown in popularity and in the scope of material covered since 1989. They no longer necessarily stem from short courses; rather, they are often generated independently by experts in the field. They are popular because they provide a ready reference to those wishing to learn about emerging technologies or the latest information within their field. The topics within the series have grown from the initial areas of geometrical optics, optical detectors, and image processing to include the emerging fields of nanotechnology, biomedical optics, fiber optics, and laser technologies. Authors contributing to the TT series are instructed to provide introductory material so that those new to the field may use the book as a starting point to get a basic grasp of the material. It is hoped that some readers may develop sufficient interest to take a short course by the author or pursue further research in more advanced books to delve deeper into the subject. The books in this series are distinguished from other technical monographs and textbooks in the way in which the material is presented. In keeping with the tutorial nature of the series, there is an emphasis on the use of graphical and illustrative material to better elucidate basic and advanced concepts. There is also heavy use of tabular reference data and numerous examples to further explain the concepts presented. The publishing time for the books is kept to a minimum so that the books will be as timely and uptodate as possible. Furthermore, these introductory books are competitively priced compared to more traditional books on the same subject. When a proposal for a text is received, each proposal is evaluated to determine the relevance of the proposed topic. This initial reviewing process has been very helpful to authors in identifying, early in the writing process, the need for additional material or other changes in approach that would serve to strengthen the text. Once a manuscript is completed, it is peer reviewed to ensure that chapters communicate accurately the essential ingredients of the science and technologies under discussion. It is my goal to maintain the style and quality of books in the series and to further expand the topic areas to include new emerging fields as they become of interest to our reading audience. James A. Harrington Rutgers University Contents Preface .................................................................................................. xiii Chapter 1 Analytic Fourier Theory Review .......................................... 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 A Little History and Purpose ....................................................................1 The Realm of Computational Fourier Optics............................................2 Fourier Transform Definitions and Existence ...........................................3 Theorems and Separability .......................................................................3 Basic Functions and Transforms ...............................................................5 Linear and SpaceInvariant Systems .........................................................7 Exercises .................................................................................................10 References ...............................................................................................12 Chapter 2 Sampled Functions and the Discrete Fourier Transform 13 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 Sampling and the Shannon–Nyquist Sampling Theorem .......................13 Effective Bandwidth ...............................................................................15 Discrete Fourier Transform from the Continuous Transform .................18 Coordinates, Indexing, Centering, and Shifting ......................................20 Periodic Extension ..................................................................................21 Periodic Convolution ..............................................................................24 Exercises .................................................................................................26 References ...............................................................................................27 Chapter 3 MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms ................................................................... 29 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 Defining Functions .................................................................................29 Creating Vectors .....................................................................................32 Shift for FFT ...........................................................................................34 Computing the FFT and Displaying Results ...........................................36 Comparison with Analytic Results .........................................................38 Convolution Example .............................................................................39 Two Dimensions .....................................................................................41 Miscellaneous Hints ................................................................................43 Exercises .................................................................................................45 ix x Contents Chapter 4 Scalar Diffraction and Propagation Solutions ................ 47 4.1 4.2 4.3 4.4 Scalar Diffraction ....................................................................................47 Monochromatic Fields and Irradiance ....................................................48 Optical Path Length and Field Phase Representation .............................50 Analytic Diffraction Solutions ................................................................51 4.4.1 Rayleigh–Sommerfeld solution I...................................................51 4.4.2 Fresnel approximation ...................................................................53 4.4.3 Fraunhofer approximation .............................................................55 4.5 Fraunhofer Diffraction Example .............................................................56 4.6 Exercises .................................................................................................59 4.7 References ...............................................................................................61 Chapter 5 Propagation Simulation ..................................................... 63 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 Fresnel Transfer Function (TF) Propagator ............................................63 Fresnel Impulse Response (IR) Propagator ............................................64 Square Beam Example ............................................................................66 Fresnel Propagation Sampling ................................................................69 5.4.1 Square beam example results and artifacts ....................................69 5.4.2 Sampling regimes and criteria .......................................................72 5.4.3 Criteria applied to square beam example ......................................74 5.4.4 Propagator accuracy ......................................................................75 5.4.5 Sampling decisions ........................................................................77 5.4.6 Splitstep simulation, windowing, and expanding grids................78 Fraunhofer Propagation ..........................................................................79 Coding Efficiency ...................................................................................83 Exercises .................................................................................................83 References ...............................................................................................86 Chapter 6 Transmittance Functions, Lenses, and Gratings ............. 89 6.1 6.2 6.3 6.4 Tilt...........................................................................................................89 Focus .......................................................................................................93 Lens.........................................................................................................96 Gratings and Periodic Functions .............................................................98 6.4.1 Cosine magnitude example ...........................................................99 6.4.2 Squarewave magnitude example ................................................102 6.4.3 Onedimensional model ..............................................................105 6.4.4 Periodic model.............................................................................106 6.5 Exercises ...............................................................................................108 6.6 References .............................................................................................111 Chapter 7 Imaging and DiffractionLimited Imaging Simulation ... 113 7.1 Geometrical Imaging Concepts.............................................................113 7.2 Coherent Imaging .................................................................................116 Contents xi 7.2.1 Coherent imaging theory .............................................................116 7.2.2 Coherent transfer function examples ...........................................117 7.2.3 Diffractionlimited coherent imaging simulation ........................119 7.2.4 Rough object ...............................................................................124 7.3 Incoherent Imaging ...............................................................................127 7.3.1 Incoherent imaging theory...........................................................127 7.3.2 Optical transfer function examples ..............................................128 7.3.3 Diffractionlimited incoherent imaging simulation .....................129 7.4 Exercises ...............................................................................................132 7.5 References .............................................................................................139 Chapter 8 Wavefront Aberrations ..................................................... 141 8.1 Wavefront Optical Path Difference ......................................................141 8.2 Seidel Polynomials ...............................................................................142 8.2.1 Definition and primary aberrations .............................................142 8.2.2 MATLAB function ......................................................................144 8.3 Pupil and Transfer Functions ................................................................146 8.3.1 Pupil function ..............................................................................146 8.3.2 Imaging transfer functions...........................................................147 8.4 Image Quality .......................................................................................147 8.4.1 Point spread function ...................................................................147 8.4.2 Modulation transfer function .......................................................148 8.5 Lens Example—PSF and MTF .............................................................148 8.6 Wavefront Sampling .............................................................................153 8.7 Superposition Imaging Example ...........................................................157 8.7.1 Image plane PSF map ..................................................................157 8.7.2 Image simulation .........................................................................160 8.7.3 Practical image simulation ..........................................................163 8.8 Exercises ...............................................................................................163 8.9 References .............................................................................................168 Chapter 9 Partial Coherence Simulation .......................................... 169 9.1 Partial Temporal Coherence .................................................................170 9.1.1 Quasimonochromatic light .........................................................170 9.1.2 Partial temporal coherence simulation approach .........................172 9.1.3 Partial temporal coherence example ............................................173 9.2 Partial Spatial Coherence ......................................................................177 9.2.1 Stochastic transmission screen ....................................................177 9.2.2 Partial spatial coherence simulation approach ............................178 9.2.3 Partial spatial coherence example ...............................................182 9.3 Reducibility, Number of Spectral Components, and Phase Screens .....186 9.4 Exercises ...............................................................................................187 9.5 References .............................................................................................189 xii Contents Appendix A Fresnel Propagator Chirp Sampling ............................ 191 A.1 Fresnel Transfer Function Sampling ....................................................191 A.1.1 Oversampled transfer function .................................................192 A.1.2 Critically sampled transfer function .........................................194 A.1.3 Undersampled transfer function ...............................................194 A.2 Fresnel Impulse Response Function Sampling ....................................195 A.2.1 Undersampled impulse response ..............................................196 A.2.2 Critically sampled impulse response ........................................196 A.2.3 Oversampled impulse response ................................................197 A.3 Summary ..............................................................................................198 A.4 References ............................................................................................198 Appendix B Fresnel TwoStep Propagator ...................................... 199 B.1 Approach ..............................................................................................199 B.2 Sampling Considerations .....................................................................202 B.2.1 Similar side lengths....................................................................203 B.2.2 Significantly different side lengths ............................................203 B.2.3 Comments and recommendations ..............................................204 B.3 MATLAB Code ...................................................................................204 B.4 References ............................................................................................205 Appendix C MATLAB Function Listings .......................................... 207 C.1 C.2 C.3 C.4 C.5 C.6 Circle ....................................................................................................207 Jinc .......................................................................................................207 Rectangle ..............................................................................................207 Triangle ................................................................................................208 Unit Sample “Comb” ...........................................................................208 Unit Sample “Delta” ............................................................................208 Appendix D Exercise Answers and Results .................................... 209 D.1 D.2 D.3 D.4 D.5 D.6 D.7 D.8 D.9 Chapter 1 .............................................................................................209 Chapter 2 .............................................................................................210 Chapter 3 .............................................................................................211 Chapter 4 .............................................................................................212 Chapter 5 .............................................................................................214 Chapter 6 .............................................................................................217 Chapter 7 .............................................................................................220 Chapter 8 .............................................................................................223 Chapter 9 .............................................................................................225 Index ..................................................................................................... 229 Preface This book began as a collection of notes and computer examples prepared for a firstyear graduate course on Fourier optics. In teaching Fourier optics over a number of years, I found that I developed a better conceptual understanding of the analytic material after setting up examples for the class on the computer. The examples required careful consideration of the sample coordinates, amplitude scaling, practical dimensions, display settings, sampling conditions, and a number of other issues. It wasn’t long before I started designing computer exercises for the students to do—figuring that if it helped me, it would probably help them. In addition, applying the theory to produce a display of a beam pattern or a blurry image of some object seemed to bring the application of Fourier optics to life for many students. At the same time, the research being performed by my group at New Mexico State University involved wave optics simulation of laser beam propagation through atmospheric turbulence. The synergy of the teaching and research activities led to the idea of a book on computer methods and Fourier optics. I did some research and found a scattering of material on numerical Fourier optics, but no book with the content I envisioned. So with that, the project began. Computational Fourier Optics is a text that shows the reader in a tutorial form how to implement Fourier optical theory and analytic methods on the computer. A primary objective is to give students of Fourier optics the capability of programming their own basic wave optic beam propagations and imaging simulations. The book will also be of interest to professional engineers and physicists learning Fourier optics simulation techniques—either as a selfstudy text or a text for a short course. For more advanced study, the latter chapters and appendices provide methods and examples for modeling beams and pupil functions with more complicated structure, aberrations, and partial coherence. For a student in a course on Fourier optics, I envision this book as a companion to any of several excellent textbooks on Fourier optical theory. I felt a companion book should be concise, accessible, and practical—so those are also goals for this text. The book begins in Chapter 1 with a short review of the Fourier optical results that are central to wave optics simulation development. The review is intended to be a quick, consolidated reference. In Chapter 2 the discrete Fourier transform (DFT) is developed, of which the fast Fourier transform (FFT) version is a primary tool for simulations. FFT scaling aspects, index formatting, and other differences from the analytic xiii xiv Preface transform are introduced. These differences later come to play in the scaling and interpretation of the simulation results. The handson tutorial part of the book begins in Chapter 3 where stepbystep examples are presented that involve programming functions, vectors, equations, and taking transforms in MATLAB®. Students with a range of backgrounds—electrical engineers, astronomers, physicists—take my Fourier optics course. The nonengineers often have never used MATLAB, so the idea of combining a MATLAB tutorial with a computational Fourier optics tutorial was natural and led to Chapter 3. The MATLAB programming environment is optimized for vector and matrix operations; therefore, it is a good tool for Fourier optics simulation, which generally involves at least two dimensions. Furthermore, MATLAB has a heritage in this subject since several optical propagation codes, such as the AOTools and WaveProp toolboxes, are written in MATLAB. The material in this chapter has been tested by students in my Fourier optics course, and even those without any MATLAB experience have found they could get up and going quickly with the tutorial. Chapter 4 is a quick review and summary of scalar diffraction and optical propagation theory. The expressions presented in Chapter 4 are taken into the computer domain in Chapter 5. Implementations of the Fresnel and Fraunhofer diffraction expressions are described with stepbystep coding instructions. The methods are demonstrated for an illuminated aperture. Attention is paid to sampling issues that can be the bane of wave optics propagation simulations. Chapter 6 covers techniques that add further application to the diffraction simulations. Methods are described for applying tilt and focus to an optical wavefront, and lenses and diffraction gratings are considered. A review of coherent and incoherent imaging theory and modeling techniques applied to diffractionlimited imaging examples are presented in Chapter 7. Imaging simulation is extended in Chapter 8 to the more practical circumstance involving wavefront aberrations. Chapter 9 provides a short review of coherence theory and demonstrates approaches for simulating partial temporal and partial spatial coherent illumination. Exercises at the end of each chapter (with answers in the back of the book) give the reader a chance to work with both theory and computer implementations. The appendices cover: (a) further sampling details for Fresnel diffraction; (b) a twostep diffraction propagation technique that allows arbitrary grid scaling between the source and observation planes; (c) listings of basic MATLAB functions developed in the text; and (d) answers to the exercises. Please visit http://www.ece.nmsu.edu/~davvoelz/cfo/ for updates, errata, files and other resources. This book would have never happened were it not for a sabbatical leave in 2008 in the Upper Peninsula of Michigan. I owe Mike Roggemann at Michigan Tech a big debt for all of his care and feeding of a displaced New Mexican. My discussions with him, on and off the lake, helped shape much of the content of Preface xv this book. I also thank the faculty and staff at Michigan Tech for all their support. Go Huskies! As this project was getting underway, Jason Schmidt kindly sent a first draft of his book Numerical Simulation of Optical Wave Propagation with Examples in MATLAB®. I tried to avoid studying it too closely as I wanted to put my own spin on related material. But I had to peek from time to time to see what he had to say on certain matters. His book was a valuable resource. Xifeng Xiao at New Mexico State University deserves credit for pioneering much of the partial coherence material. She also combed through all the chapters, working examples and checking equations. Our discussions over the years on numerical simulation are deeply imbedded in this book. It has been a great pleasure to work with her. The students of a succession of Fourier optics courses since 2003 at New Mexico State University have been, often unknowingly, a constant source of insight and inspiration for this book. Their reactions and feedback to the material helped change many things for the better and encouraged me to keep going. For all those spurofthemoment questions and sudden inquiries of howdoesthatwork, I thank my colleagues at the Klipsch School of ECE at New Mexico State University, especially Deva Borah, Laura Boucheron, Chuck Creusere, Philip DeLeon and Mike Giles  a good group of folks. Finally I cannot thank my wife, Judi, enough for supporting this project in every way, including proofreading the manuscript. Our children, Alex, Katie and Brian, have had to deal with an absent dad while I worked on this book, so I thank them for their patience. My family is my support and I couldn’t do what I do without them! David Voelz December 2010 Chapter 1 Analytic Fourier Theory Review 1.1 A Little History and Purpose The branch of optical science known today as “Fourier optics” had its genesis in the 1940s through the 1960s with the application of new telecommunications and circuit design analysis techniques in optical diffraction theory.1 In 1968 this upstart discipline was given a permanent foothold with the publication of Introduction to Fourier Optics, by Joseph W. Goodman, a seminal textbook that explained and united the fundamental concepts, and which continues to add significantly to the application of Fourier optics in subsequent editions.2 Fourier optics is now the cornerstone for the analysis of diffraction, coherence, and imaging, as well as specialized topics such as wavefront control, propagation through random media, and holography. The study of Fourier optics today leads naturally toward the computer for at least two reasons: (1) diffraction integral expressions are difficult to solve analytically for all but a few of the simplest aperture functions, and (2) the fast Fourier transform (FFT) algorithm combined with the linear systems framework of Fourier optics provides an extremely efficient computational approach for solving wave optics problems. Certainly, the computer can be applied directly in finding exceedingly accurate solutions to diffraction problems using numerical integration techniques.3 However, this book is really about the FFT and how to apply it to a variety of Fourier optics problems. The computer coding steps mirror the analytic concepts and the FFT’s speed makes it possible to perform thousands of optical propagation or imaging simulations in a reasonable amount of time. In fact, the methods explored in this book form the basis for wave (or physical) optics simulation tools that are widely used in industry. But, of course, there’s no free lunch (…if there were, perhaps we could be eating while studying Fourier optics...). It turns out the FFT is an accomplice to various numerical artifacts. We do our best in this book to expose these issues and provide constraints to help minimize the damage. This is also a tutorial text with stepbystep instructions, not only for coding Fourier optics problems, but also for MATLAB, our software application of choice. So, if you are new to MATLAB, don’t worry! Chapter 3 starts at the 1 2 Chapter 1 beginning (“Open MATLAB”) and leads you through the basics of working with the FFT. By the end of the book you will be programming diffraction problems involving partially coherent light—at least that’s the goal! Exercises at the end of the chapters give you room to tinker with the programs and stretch out with your own code. It is assumed the reader has some familiarity with Fourier optics. Presenting the topic from the ground up is too much material to cover and would obscure our purpose. However, the analytic theory required is presented in summary form throughout the text. The notation and form closely follow Goodman’s presentation in Introduction to Fourier Optics.2 For further details and explanations of the analytic foundations of Fourier theory and Fourier optics the reader is encouraged to consult Goodman’s book as well as the many other excellent references that exist on the topic.4–7 1.2 The Realm of Computational Fourier Optics In this book, the variables, vectors, and arrays in the computer code are defined as much as possible in terms of physical quantities. For example, the coordinates of samples in an array that models a spatial plane are defined in units of meters. Integers for indexing arrays show up only when they can’t be avoided. This approach allows a clear connection between the physical world being modeled and the computer code. MATLAB’s vectorized structure is also suited to this approach. Thus, programming examples presented in the book involve specific aperture sizes, wavelengths, and distances. Although some examples are simply academic, others are something one might encounter in the real world. However, the reader will soon notice an emergent theme: the finite size of the sample array in the computer limits the range of parameters that can be considered. We might consider this difficulty in light of the optical designer’s dilemma: When does one transition between a geometrical optics prediction of system performance and a wave optics prediction? The difference between these predictors is that geometrical optics assumes rectilinear (straightline) propagation of the rays of light and ignores diffractive spreading due to the wave nature of light. The usual answer for the dilemma is that for small departures from perfection (near the “diffraction limit”) a wave optics description is needed. For large departures a geometrical ray optics description, which has more flexible implementation options, is adequate.8,9 So, although analytic Fourier optics theory is quite general, the finite array size tends to limit the computer modeling to the “nearperfection” situations. Typically, this means small divergence angles for optical beam propagation, small simulated image area, and so forth. For practical applications, this is the same realm as the wave optics performance prediction for optical system design. The remainder of this chapter is a summary of the fundamental Fourier transform definitions, theorems, basic functions, and transform pairs. A review of linear systems theory is also included. So, let’s go! Analytic Fourier Theory Review 3 1.3 Fourier Transform Definitions and Existence Fourier optics problems often involve two spatial dimensions. The analytic Fourier transform of a function g of two variables x and y is given by G f X , fY g ( x, y ) exp j 2π f X x fY y dxdy , (1.1) where G(fX , fY) is the transform result and fX and fY are independent frequency variables associated with x and y, respectively. This operation is often described in a shorthand manner as g ( x, y ) G ( f X , f Y ) . Similarly, the analytic inverse Fourier transform is given by g x, y G( f X , fY ) exp j 2π f X x fY y df X dfY . (1.2) The shorthand notation for this operation is 1 G( f X , f Y ) g ( x, y ). For the Fourier transform to be realizable in a mathematical sense, g(x,y) must satisfy certain sufficient conditions. These conditions are commonly listed as: (a) g must be absolutely integrable over the infinite range of x and y; (b) g must have only a finite number of discontinuities; and (c) g must have no infinite discontinuities. Goodman2 illustrates that in a number of important cases, one or more of these conditions can be weakened, and a generalized transform approach using idealized mathematical functions can be employed to find useful transform representations. Some generalized transform results of interest include 1 ( f X , fY ) , cos 2πf 0 x 12 ( f X f 0 , fY ) 12 ( f X f 0 , fY ) , where is the Dirac delta function. 1.4 Theorems and Separability The theorems listed in Table 1.1 find considerable application in Fourier analysis. In Table 1.1, A, B, a, and b are scalar constants. An important property of certain functions is separability. A twodimensional (2D) function is separable if it can be written as the product of two functions of a single variable, such as 4 Chapter 1 g S x, y g X x g Y y . (1.3) Separability reduces the Fourier transform of a 2D function to the product of two onedimensional (1D) transforms or g S x, y g X x g Y y . (1.4) Table 1.1 Fourier transform theorems. Theorem Expression Linearity Ag x, y Bhx, y Ag x, y Bh x, y Similarity x y g , ab G af X , bf Y a b Shift g ( x a, y b) G f X , fY exp j 2π f X a fY b Parseval’s (Rayleigh’s) g ( x, y) Convolution Autocorrelation 2 dxdy G ( f X , f Y ) df X df Y 2 g ( , )h( x , y )d d G ( f g ( , ) g ( x, y)d d G( f , f ) g ( x, y ) G ( , )G ( f , f )d d 2 X 2 Y X Crosscorrelation , fY ) H ( f X , f Y ) X Y g ( , )h ( x, y)d d G ( f X , fY ) H ( f X , fY ) g ( x, y ) h ( x, y ) G ( , ) H ( f X , f Y ) dd Fourier integral 1 g ( x, y ) 1 g ( x, y ) g ( x, y ) Successive transform g ( x, y ) g ( x, y ) Central ordinate g ( x, y ) f X 0 G (0,0) g ( x, y )dxdy fY 0 1 G( f X , f Y ) x 0 g (0,0) G ( f X , f Y )df X df Y Note: A, B, a, and b are scalar constants y 0 Analytic Fourier Theory Review 5 1.5 Basic Functions and Transforms Several basic functions, or combinations thereof, are used to describe various physical or analytic structures encountered in optics, such as a circle function to describe a circular aperture. Thus, these functions and their Fourier transform pairs are of considerable utility. The definitions in Table 1.2 are adopted. Functions of one variable are illustrated in Fig. 1.1. These can be combined as products to represent separable 2D functions. The circle function is a symmetric 2D function where a single radial variable r = (x2+y2)1/2 is often used. A shorthand name is not defined for the Gaussian, but this function appears often. The form we use is convenient for Fourier analysis. The circle and a 2D Gaussian function are plotted in Fig. 1.2 for illustration. Table 1.2 Basic functions. Function Definition Rectangle 1 1, x 2 1 1 rect( x) , x 2 2 0, otherwise Sinc sinc( x) Triangle 1 x , x 1, ( x) otherwise. 0, Comb comb( x ) sin(πx) πx ( x n) n Gaussian Circle exp πx 2 circ x2 y 2 1 x2 y 2 , 1 2 1 1 x2 y 2 , 2 2 0 otherwise. 6 Chapter 1 Figure 1.1 Basic 1D functions. Figure 1.2 Examples of 2D functions. Analytic Fourier Theory Review 7 If optical structures and apertures are modeled with basic functions, then corresponding Fourier transforms can aid in finding diffraction solutions or image results. The basic functions and their Fourier transforms are presented in Table 1.3. J1 is a Bessel function of the first kind, order 1, and appears in the transform of the circle function. The transform of the circle is illustrated in Fig. 1.3. In Table 1.3, the last row gives a pair of “chirp” functions that will become quite familiar in the following chapters. 1.6 Linear and SpaceInvariant Systems The power of Fourier methods to analyze the response of a physical system to an input is significantly enhanced if the system can be modeled as linear and shift(or space) invariant. There are many aspects of optical systems that can be modeled in this way. In general, the operation of a system on a twovariable input Table 1.3 Basic functions and their transforms. Function x rect a Transform a sincaf X x sinc a a rectaf X x a a sinc 2 af X x comb a a combaf X x2 exp π 2 a a exp πa 2 f X 2 x2 y2 circ a a x2 y 2 exp π 2 2 b a ab exp π a 2 f X2 b 2 fY2 x2 y 2 exp jπ 2 2 b a j ab exp jπ a 2 f X2 b 2 fY2 2 J1 2πa f X2 fY2 a f X2 fY2 Note: J1 is a Bessel function of the first kind, order 1. 8 Chapter 1 J1 2π f X2 fY2 / f X2 fY2 Figure 1.3 Circle function transform; peak value at fX = 0, fY = 0 is . function g1 to produce an output function g2 can be described by g 2 x2 , y 2 S g1 x1 , y1 , (1.5) where S indicates the operation performed by the system. The “test” for linearity is the following: S Ag A x1 , y1 Bg B ( x1 , y1 ) ASg A x1 , y1 B S g B ( x1 , y1 ) , (1.6) where A and B are scalar constants. For a sum of input functions—for example, gA and gB in Eq. (1.6) with constant multipliers—the output is a sum of the individual responses. If the input can be “decomposed” into a sum of “elementary functions,” then the output of a linear system can be determined if the response to the elementary functions is known. Linearity leads to the following expression, known as a superposition integral: g 2 x2 , y2 g , h x , y ; , d d . 1 2 2 (1.7) The function h is the impulse response of the system, and the integrals indicate that the output of the system is a superposition—or sum—of an infinite set of impulse responses that multiply the input function. The impulse response is modeled by hx 2 , y 2 ; , S x1 , y1 . (1.8) Analytic Fourier Theory Review 9 A linear system is completely characterized by its responses to impulse functions, but to use this property in practice the responses must be known for all locations in the input plane (x1, y1). Linearity represents one level of simplification. Further simplification is afforded by the property of space invariance, where in its most basic form we write g 2 x 2 , y 2 S g 1 x1 , y1 , (1.9) and, therefore, the impulse response simplifies to h x 2 , y 2 ; , h x 2 , y 2 . (1.10) This impulse response does not depend on the absolute position in the input plane or the output plane. It only depends on the relative separation of the input and output points as if they were to appear in a common x–y plane. An interpretation of this situation is that an impulse anywhere in the input plane creates a corresponding response in the output plane that changes position with the input but always has the same relative form. The superposition integral now becomes a convolution integral: g 2 x2 , y2 g , h x 1 2 , y2 d d . (1.11) In shorthand notation with the convolution operator , Eq. (1.11) is written as g 2 x, y g 1 x , y h x, y , (1.12) where the subscripts on the x and y variables are no longer necessary. Taking the Fourier transform of each function in Eq. (1.12) and applying the convolution theorem yields G 2 f X , f Y G1 f X , f Y H f X , f Y , (1.13) where H(fX, fY) is the Fourier transform of the impulse response h(x, y) and is known as the transfer function. Two valuable features of linear spaceinvariant systems are apparent from Eq. (1.13): 1. Rather than directly tackling the convolution integral, a more computationally appealing route can be taken of transforming the input to the Fourier domain, multiplying by the transfer function, and inverse transforming to find the result. 2. The transfer function model is analogous to frequencyfiltering operations that are found in electric circuit theory, digital signal 10 Chapter 1 processing, and many other disciplines involving signal analysis. Insights from these areas can often be applied to Fourier methods associated with optical systems. 1.7 Exercises 1.1 Sketch the following functions: (a) x rect , 2 (b) rectx 2 , (c) x 2 , 2 (d) exp 3πx2 , (e) x x comb 4 x rect 12 , (f) circ x 22 y 2 circ x 22 y 2 . 1.2 Using known transform pairs and theorems, find the Fourier transforms of the following: (a) y x rect rect , 2w 2w (b) x x0 y rect rect , 2w 2w (c) x2 y2 exp w2 (d) x2 y2 x2 y2 , circ circ w1 w2 (e) circ , x d / 2 2 y 2 w circ x d / 22 y 2 w . 1.3 Perform the following convolutions by applying the convolution theorem: Analytic Fourier Theory Review 11 (a) y y x x rect rect rect rect , 2w 2w 2w 2w (b) x2 y 2 x2 y 2 exp π exp π , 32 42 (c) x x sinc sinc y sinc sinc y . 2 4 1.4 Find the autocorrelations of the following: (a) y x rect rect , 2w 2w (b) x2 y 2 exp π . w2 1.5 Apply the central ordinate theorem to find g ( x, y) dxdy for the following, and compare the results with simple area calculations: (a) y x g ( x, y) rect rect , 2w 4w (b) x2 y2 g ( x, y) circ 3 . 1.6 Demonstrate whether the following operations are linear and/or space invariant, where A, B, a, and b are scalar constants: (a) S g ( x, y) Ag( x, y) , (b) S g ( x, y) Ag( x, y) +B, (c) S g ( x, y) Ag ( x, y) , (d) S g ( x, y) xgx, y , (e) 2 Ave g ( x, y ) 1 ab x a b y 2 2 x a b y 2 2 g ( , ) d d . 12 Chapter 1 1.8 References 1. W. T. Rhodes, “History and evolution of the teaching of Fourier optics,” Proc. SPIE, 3572, 50–56 (1999). [doi:10.1117/12.358418]. 2. J. W. Goodman, Introduction to Fourier Optics, 3rd Ed., Roberts & Company, Greenwood Village, CO (2005). 3. R. Barakat, “The calculation of integrals encountered in optical diffraction theory,” in The Computer in Optical Research, Topics in Applied Physics, Vol. 41, B. R. Frieden (ed.), Springer, Berlin (1980). 4. J. D. Gaskill, Linear Systems, Fourier Transforms, and Optics, WileyInterscience, New York (1978). 5. O. K. Ersoy, Diffraction, Fourier Optics and Imaging, WileyInterscience, New York (2006). 6. C. A. Bennett, Principles of Physical Optics, Wiley, Hoboken, NJ (2008). 7. E. G. Steward, Fourier Optics: An Introduction, 2nd Ed., Dover, New York (2004). 8. W. J. Smith, Modern Optical Engineering, 4th Ed., McGrawHill Professional, New York (2007). 9. J. M. Geary, Introduction to Lens Design with Practical ZEMAX® Examples, WillmannBell, Richmond, VA (2002). Chapter 2 Sampled Functions and the Discrete Fourier Transform When implementing Fourier optics simulations on the computer it is necessary to represent functions by discrete arrays of sampled values and apply transform and processing methods designed for these discrete signals. To come as close as possible to simulating continuous space, it would be great to model the physical elements with a gazillion samples. However, computer memory and execution time limitations won’t allow this. Thus, devising practical Fourier optics simulations becomes an act of balancing acceptable sampling artifacts and available computer resources. This chapter begins to address this matter with discussions of the sampling of continuous functions, the Shannon–Nyquist sampling theorem, and the concept of effective bandwidth. The remainder of the chapter concerns the discrete Fourier transform (DFT), the workhorse tool for computational Fourier analysis. We look at its relationship with the analytic transform, describe implementation details, and discuss how DFT results differ from the analog world. 2.1 Sampling and the Shannon–Nyquist Sampling Theorem Consider the twodimensional (2D) analytic function g(x,y) and suppose it is sampled in a uniform manner (Fig. 2.1) in the x and y directions, which is indicated by g ( x, y ) g ( mx, ny ) , (2.1) where the sample interval is x in the x direction and y in the y direction, and m and n are integervalued indices of the samples. The respective sample rates are 1/x and 1/y. In practice, the sampled space is finite and, assuming it is composed of M N samples in the x and y directions, respectively, m and n are often defined with the following values: m M M , , 1 , 2 2 13 n N N ,, 1 . 2 2 (2.2) 14 Chapter 2 x y y y x x LY LX (a) (b) Figure 2.1 Twodimensional function: (a) analytic and (b) sampled versions. This is a standard index arrangement where M and N are assumed to be even. Even numbers of samples are used in this book for reasons associated with discrete Fourier transform efficiency and sample arrangement (see Section 2.4). A finite physical area (e.g., units of m2) is spanned by the sampled space, and this is given by LX LY, where LX is the length along the x side of the sampled space and LY is the length along the y side (Fig. 2.1). LX and LY are referred to as the side lengths. They represent physical distances and are related to the sampling parameters by L X Mx , LY Ny . (2.3) An obvious sampling concern is whether all the significant values of g(x,y) “fit” within the physical area defined by LX LY. The support of g(x,y) refers to the span of the significant values. This concept is illustrated in Fig. 2.2(a) for one axis. If DX is the support in the x direction and DY is the support in the y direction, then for the significant values of g(x,y) to be contained within the array requires DX LX , DY LY . (2.4) Another concern is whether the sample intervals are small enough to preserve the features of g(x,y). For functions that are bandlimited, where the spectral content of the signal is limited to a finite range of frequencies, a continuous function can be recovered exactly from the samples if the sample interval is smaller than a specific value. The Shannon–Nyquist sampling theorem, extended to two dimensions, states this requirement as1 x 1 , 2B X y 1 , 2 BY (2.5) Sampled Functions and the Discrete Fourier Transform 15 G(fX, fY) g(x, y) x ~ BX ~ DX (a) fX (b) Figure 2.2 Illustration of the (a) support DX and (b) bandwidth BX along the x axis of g(x, y). Bandwidth is commonly defined as a halfwidth measure and is illustrated here with a profile of G(fX, fY), the Fourier transform magnitude of g(x, y). where BX is the bandwidth of the spectrum of the continuous function along the x direction and BY is the bandwidth along the y direction. Bandwidth is illustrated in Fig. 2.2(b). Violating Eq. (2.5) results in aliasing, in which undersampled highfrequency components in the signal are interpreted erroneously as lowfrequency content. This issue is considered further in Section 2.5. A related parameter is the Nyquist frequency given by f NX 1 , 2 x f NY 1 , 2 y (2.6) which is half the sample rate and corresponds to the maximum spatial frequency that can be adequately represented given the interval x or y. 2.2 Effective Bandwidth Practical functions such as those defined in Chapter 1 are not bandlimited. In fact, any function with finite support, like the rectangle or circle functions, cannot be bandlimited. Often these functions have an effective bandwidth that encompasses the most significant frequency values. Even though the criteria posed by the Shannon–Nyquist theorem may not be completely satisfied, a small enough sample interval can be found to provide an acceptable representation of the analytic function where the effects of aliasing are small. For example, consider a 2D square signal with halfwidth of w: y x f ( x, y ) rect . rect 2w 2w The analytic Fourier transform yields the spectrum (2.7) 16 Chapter 2 F f X , fY 4 w2 sinc 2 wf X sinc 2 wfY . (2.8) One approach to defining the effective bandwidth is to find the spectral width (radius) that contains a high percentage of the total power in the spectrum. Applying Parseval’s theorem, the total spectral power of Eq. (2.8) is 4w PT 2 2 sinc 2 2 wf X sinc 2 2 wfY df X dfY rect 2 x 2 y 2 rect dxdy 4w 2w 2w . (2.9) A practical criterion for the effective bandwidth B is to include 98% of the total spectral power. Converting to polar coordinates to allow a radial bandwidth value B to be considered leads to 1 PT 2π B 0 0 4w 2 2 sinc 2 2 w cos sinc 2 2 w sin d d 0.98 ,(2.10) where fX = cos and fY = sin . The integrals on the left side can be evaluated numerically for various values of B until Eq. (2.10) is satisfied. With this approach the effective bandwidth is found to be B 5 . w (2.11) Figure 2.3 illustrates the portion of the spectrum that encompasses 98% of the spectral power. Substituting Eq. (2.11) into Eq. (2.5) for the bandwidth gives x w , 10 (2.12) which says at least 10 samples across the halfwidth of the rect function (20 across the full width) are required to retain the effective bandwidth indicated in Eq. (2.11). It is important to realize that the part of the analytic spectrum that lies beyond the Nyquist frequency does not simply disappear. Even though small in power, it can introduce noticeable aliased frequency content that is erroneous. Table 2.1 shows effective bandwidth values for square, circle, and Gaussian functions computed in the same way. A larger effective bandwidth can be used if there is a need to include more of the spectral power. Sampled Functions and the Discrete Fourier Transform 17 5/w 5/w 0 fy 0 5/w 5/w fx (a) 5/w 5/w 0 fy 0 5/w 5/w fx (b) Figure 2.3 (a) Magnitude of the Fourier spectrum and (b) power spectrum of g(x,y) = rect(x/w) rect(y/w), comprising 98% of the total spectral power. 18 Chapter 2 2.3 Discrete Fourier Transform from the Continuous Transform The DFT, usually in the form of its highly efficient offspring—the fast Fourier transform (FFT)—is a fundamental tool for modeling Fourier optics problems on the computer. The objective for this section is to develop the DFT, a discrete implementation of the continuous Fourier transform. The derivation is helpful for understanding the scaling of the spatial sample coordinates and frequency coordinates, as well as the constants multiplying a discrete transform result. This scaling is an important part of modeling a physical optics problem. Only the aspects of the DFT and FFT that are critical to the simulation approaches covered in this book are highlighted, so there are many more details to be discovered (or rediscovered) in other resources.2–4 The analytic Fourier transform of a function g of two variables x and y is repeated here for reference G f X , fY g ( x, y ) exp j 2π f X x fY y dxdy . (2.13) First, assume g(x,y) is sampled as indicated in Eqs. (2.1) and (2.2). To simplify some of the notation, the following substitution can be used where the actual sample intervals are not explicitly shown: g ( m x , n y ) g~ ( m , n ) . (2.14) Next, the integrals in Eq. (2.13) can be approximated using a Riemann sum: dxdy N / 2 1 M / 2 1 x y . (2.15) n N / 2 m M / 2 Table 2.1 Effective bandwidth for 98% power. Function Effective Bandwidth B y x rect rect 2w 2w 5 w x2 y2 circ w 5 w x2 y 2 exp π 2 w 0.79 w Sampled Functions and the Discrete Fourier Transform 19 Because the DFT operation is performed generically on a discrete array of values without specific sample interval information, the multipliers xy in Eq. (2.15) are not included in the DFT definition. However, these multipliers need to be applied subsequently to the DFT operation for appropriate scaling of a physical problem. The convention for the frequency domain is to divide this continuous “space,” indicated by fX and fY, into M and N evenly spaced coordinate values. This involves the following substitutions: fX p , M x where p M M , , 1; 2 2 q , N y where q N N , , 1; 2 2 fY (2.16) where p and q are integers that index multiples of the frequency sample intervals f X 1 1 , Mx L X and f Y 1 1 . Ny LY (2.17) In fact, p and q take on the same values as m and n, respectively, since the spatial and frequency arrays have the same number of elements. Note that the maximum absolute values of the frequency coordinates in Eq. (2.16) are the Nyquist frequencies 1 /( 2x) = fNX and 1 /(2y ) = fNY. Incorporating Eq. (2.16) into the complex exponential kernel of Eq. (2.13) yields p q exp j 2( f X x fY y) exp j 2 mx ny N y M x pm qn exp j 2 . (2.18) M N Finally, substituting Eqs. (2.14)–(2.18) into Eq. (2.13), we arrive at the following form of the DFT: G p, q ~ pm qn , g (m, n)exp j 2π N M m M /2 n N / 2 M /2 1 N /2 1 (2.19) where G ( p, q ) represents the DFT of g~ (m, n) . The inverse discrete Fourier transform (DFT−1) is derived in a similar way and is written as 20 Chapter 2 g m, n 1 MN pm qn . G ( p, q)exp j 2π N M p N /2 q M /2 M /2 1 N /2 1 (2.20) The appearance of the 1/MN multiplier in Eq. (2.20) requires some discussion. The factor is equal to the product xyfXfY [see Eq. (2.17)], which comes about when numerically evaluating a forward, and then an inverse Fourier integral in succession. This factor allows the DFT followed by the DFT−1 to return the original function values, which is consistent with the Fourier integral theorem. The application of 1/MN varies in different software tools. For example, MATLAB implements the transform and inverse transform based on the definitions in Eqs. (2.19) and (2.20), but some applications apply a (MN) −½ factor to both the forward and inverse transforms. In some modeling situations we will need to account for this factor. The forward and inverse DFTs are not usually accomplished with a direct execution of Eqs. (2.19) and (2.20), but rather they are accomplished with the computationally efficient FFT and FFT−1 algorithms. These algorithms implement a scheme that is not of specific importance here other than to say that the result is consistent with Eqs. (2.19) and (2.20). FFT algorithms are most efficient when M and N are a power of 2, although computation times can be nearly as fast for other values. A practical issue for FFT implementation concerns the arrangement and indexing of data values in an array. This issue is now discussed below. 2.4 Coordinates, Indexing, Centering, and Shifting Uniform sampling and square grids, where y = x, N = M, and LX = LY = L, are often used in practice. This will be the case for all the examples presented in this book; so, to simplify the presentation often only one set of variables is discussed. Considering Eqs. (2.1) and (2.2), the coordinates of the samples along one dimension can be described as L L x : x : x , 2 2 (2.21) where the above notation is borrowed from MATLAB and indicates that the coordinates range from –L/2 to L/2−x in steps of x. The yaxis coordinates are defined similarly. Assuming a FFT relationship between the spatial and spectral domains, then from Eqs. (2.16) and (2.17) the following is derived: 1 1 1 1 f X , : : 2x L L 2x (2.22) Sampled Functions and the Discrete Fourier Transform 21 which indicates that the spatial frequency coordinates range from 1/(2x) to 1/(2x)1/L in steps of fX = 1/L. Again, the fY coordinates are similar. The integer index variables m and n, as well as p and q introduced in Eqs. (2.2) and (2.16) span negative as well as positive values. However, software applications use positive integer values for vector or array (matrix) indexing. In the case of MATLAB, indexing for a onedimensional (1D) vector begins at (1). For display purposes it is convenient to “center” the function of interest in the vector, which means the zero coordinate will correspond to the (M/2+1) index. However, the convention for the 1D FFT algorithms is that the data value placed in the first index position corresponds to the zero coordinate. Thus, a “shift” of the centered vector values is needed before an FFT operation. Figure 2.4 illustrates arrangements of values and indices of a 1D sampled rect function and its spectrum. The arrangement in Figs. 2.4(a) and (b) is consistent with the analytic development in Section 2.1, where indices span negative and positive values. Figures 2.4(c) and (d) illustrate the centered arrangement in the computer with positive indices, and Figs. 2.4(e) and (f) show the shifted arrangement that is necessary prior to an FFT or FFT−1. Conversion between the centered and shifted arrangements can be done easily with the MATLAB command fftshift. Reversing the shift to get back to the centered arrangement is done with ifftshift. For a 2D array, MATLAB indexing begins at (1,1). A centered function has the zerocoordinate [x,y] = [0,0] value located at index (N/2+1, M/2+1). Prior to the 2D FFT, a shift is needed to place the zerocoordinate value at index (1,1). Figure 2.5 illustrates the centered and shifted arrangements for 2D arrays. In Fig. 2.5, row order is toptobottom, column order is lefttoright. The required shift is actually a swapping of array quadrants. Again, fftshift and ifftshift can be applied to move between the centered and shifted arrangements. A confusing detail for 2D array data is that MATLAB uses a row–column indexing scheme, where (i, j) indicates the ith row and the jth column. This is, in a sense, a reverse of standard Cartesian coordinate notation, where x (horizontal axis or “column”) is listed first and y (vertical axis or “row”) is listed second. Thus, the jindices correspond to xcoordinates and iindices correspond to ycoordinates. This explains the (j, i) listing (N/2+1, M/2+1) paired with [x,y] values in Fig. 2.5. Fortunately, MATLAB’s vector operation notation is developed for the Cartesian coordinate system; so, this issue is mostly transparent as far programming is concerned. It is only when the actual integer index values are used in codes that this array arrangement becomes an issue. 2.5 Periodic Extension Roughly speaking, all of the Fourier transform theorems listed in Table 1.1 can be applied in the discrete domain. For example, a convolution can be performed by computing the FFTs of two discrete functions, multiplying the results (pointwise) and computing the FFT−1. However, discrete transform results differ from analytic results in a way that is characterized by a concept known as periodic 22 Chapter 2 Analytic fX x Index: coords: M 2 0 L 2 0 M 2 1 2 x M 1 2 L x 2 M 1 2 1 1 2x L 0 0 (b) (a) Centered fX x Index: 1 coords: L 2 M 1 2 0 M 1 L x 2 1 2 x M 1 2 0 M 1 1 2x L (d) (c) Shifted fX x Index: 1 M 2 coords: 0 L 2 x M 2 1 M 1 M 2 L 2 x 0 1 1 2x L (e) M 2 1 M 1 1 2x L (f) Figure 2.4 Sampling arrangements for 1D spatial (left column) and frequency (right column) vectors: (a) and (b) analytic; (c) and (d) centered; and (e) and (f) shifted for FFT operations. extension. Here, we provide a short discussion of this property and an illustration. The topic is covered in more depth in other resources, such as the work by Brigham.4 Consider the 1D, analytic function f(x) and a sampled version given by x x f ( x) f ( x) comb rect . L x (2.23) Sampled Functions and the Discrete Fourier Transform (1, 1) [LX/2, LY/2] (1, 1) [0, 0] Centered 3 23 4 Shifted 1 2 (N/2+1, M/2+1) [LX/2, LY/2] (N/2+1, M/2+1) [0,0] 2 1 (a) 4 3 (b) Figure 2.5 Sampling arrangements for 2D spatial array: (a) centered and (b) shifted for the 2D FFT. The comb function “picks out” the sample values at intervals of x and the rect function sets the overall sampled space as L. Taking the analytic Fourier transform of Eq. (2.23) gives F ( f X ) F ( f X ) x combxf X L sincLf X . (2.24) This result is a continuous function where the analytic spectrum F(fX) is repeated at intervals of 1/x. The sinc convolution is a “smoothing” process. However, the FFT operation actually produces a sampled result that can be modeled by altering Eq. (2.24) with a sampling term FP ( f X ) F ( f X ) x combxf X L sincLf X L combLf X . (2.25) The new comb term sets the sample spacing in the frequency domain to be 1/L. This is consistent with Eq. (2.17). By inverse transforming Eq. (2.25), we find the function that corresponds to the spectrum FP ( f X ) : x x x f P ( x) f ( x) comb rect comb . L L x (2.26) So, the periodic extension concept can be described as follows: although we start with a sampled version of f (x) in the spatial domain, when the FFT is performed, it is as if we started with the periodic function f P ( x) and produced the periodic spectrum FP ( f X ) . To illustrate, an analytic rectangle function is shown in Fig. 2.6 (solid line) along with a sampled version (dots). The sampled version is contained in a vector of length 20, where L = 20 and x = 1. The periodic form of the function, which 24 Chapter 2 extends (virtually) beyond the original span of the sample vector, is also indicated (dashed line). Figure 2.6(b) shows the magnitude of the analytic spectrum of the rectangle (solid), the FFT result (dots), and the periodic spectrum (dashed). Figure 2.6(b) illustrates the FFT samples follow the periodic spectrum. The most obvious difference between the analytic and sample spectra in this case is slightly larger sample values in the magnitude lobes at higher frequencies. This effect results from aliasing of undersampled frequencies in the rectangle spectrum. The periodic extension concept is an instructive way to define this effect. In practice, by sampling a function to preserve the effective bandwidth— for example, the 98% level—only a small amount of aliasing is introduced. 2.6 Periodic Convolution Another issue discussed in the context of periodic extension is periodic convolution. If the FFT is applied to perform the convolution f(x)h(x), the result is actually a scaled form of fP(x)hP(x), where the periodic form fP(x) is Figure 2.6 Periodic extension illustration: (a) rect function—analytic (solid), periodic extension (dash) and sampled (dot); and (b) rect spectral magnitude—analytic (solid), periodic extension (dash) and FFT samples (dot). Sampled Functions and the Discrete Fourier Transform 25 convolved with the periodic form hP(x). This situation is illustrated in Fig. 2.7 for rectangle and right triangle functions. Sampled versions of both functions are contained in vectors of 20 samples [Figs. 2.7(a) and (b)]. Conceptually, a convolution requires one function to be reversed, or “flipped,” relative to the axis. The reversed function is translated across the second function, and the overlap area is recorded. Figure 2.7(c) shows the right triangle at a particular position in translation across the rectangle. The overlap area for an analytic convolution is the small gray triangle. However, because of periodic extension, a strip on the left side of the triangle is also overlapping a rectangle copy. This area is included in the discrete convolution result. As the righttriangle continues translating to the right, the erroneous area drops away, but the edges of the convolution result [dots in Fig. 2.7(d)] follow the periodic result. In order for the periodic convolution to match the analytic convolution, the combined support of the two functions being convolved needs to be less than the array side length, or D1 D2 L . (2.27) For the example in Fig. 2.7, D1 = 9 and D2 = 15, but L = 20; so, Eq. (2.27) is violated, and the FFTderived result does not match the analytic result. 2.7 Exercises 2.1 For a sample interval of x = 10 m and side length L = 5 mm, what is the sample number M? What is the Nyquist frequency? What is the frequency sample interval? What is the range of coordinates in the spatial domain? What is the range of the coordinates in the frequency domain? 2.2 Consider the following: (a) x2 y2 , w = 1 mm, g ( x, y ) circ w (b) x2 y2 g ( x, y ) exp w2 , w = 1 mm. For each function determine the following: (1) the effective bandwidth; (2) the maximum sample interval x necessary to satisfy the sampling theorem given the effective bandwidth; and (3) assuming 256 samples (linear dimension), the maximum side length that can be modeled. 26 Chapter 2 Figure 2.7 Periodic convolution illustration: (a) sampled rect function; (b) sampled righttriangle function; (c) shift and overlap illustration—periodic functions (dash) and sampled (dot) functions; and (d) convolution results—analytic (solid), periodic (dash) and FFT (dot). Sampled Functions and the Discrete Fourier Transform 27 2.3 What is the support of the following function along one dimension if the support is defined by where the function value drops to 1% of its peak? x2 y 2 . g x, y exp π w2 2.4 What is the bandwidth along one axis for the following? (a) y x sinc sinc , w w (b) y x sinc 2 sinc 2 . w w 2.5 Consider the following two functions: x y g 1 x, y , d 2d x2 y2 . g 2 x, y circ d (a) What is the minimum side length required to accommodate a convolution of these two functions? (b) What is the minimum side length required to accommodate the autocorrelation of g2(x, y)? 2.8 References 1. J. W. Goodman, Introduction to Fourier Optics, 3rd Ed., Roberts & Company, Greenwood Village, CO (2005). 2. R. M. Gray and J. W. Goodman, Fourier Transforms: An Introduction for Engineers, Kluwer Academic, Boston (1995). 3. A. V. Oppenheim and R. W. Schafer, DiscreteTime Signal Processing, 3rd Ed., PrenticeHall, Upper Saddle River, NJ (2009). 4. E. Brigham, The Fast Fourier Transform and its Applications, PrenticeHall, Upper Saddle River, NJ (1988). Chapter 3 MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms By following the examples presented in this chapter, you will gain some familiarity with MATLAB, learn how to implement a function, compute the discrete Fourier transform, and compare the result with analytic theory. For instructional purposes, most of the examples are onedimensional (1D) problems, but two dimensions are introduced at the end of the chapter. MATLAB is a mathematics and graphics software application with its own interpreted language that is widely used for simulation and modeling in science and engineering disciplines. It is optimized for vector and matrix operations and, therefore, is a good tool for Fourier optics simulations, which generally involve at least two dimensions. The examples in this chapter and throughout the book use a basic set of MATLAB features, in part to keep the material at a tutorial level but also because the details of the programming steps are more obvious. As you become familiar with the software, you may find more convenient and efficient ways to implement the programming. MATLAB Version 7.1 is used in this book. 3.1 Defining Functions Open MATLAB. The windows that are commonly displayed include the “Current Directory,” “Command Window,” and “Command History.” These windows are often grouped together as part of the main window that comprises the MATLAB “Desktop” (Fig. 3.1). The Current Directory shows the folder in which your work will be stored and MATLABrelated files that are in that folder. Code can be entered in the Command Window where it is executed a line at a time as it is entered. Numerical output, such as the value of a variable, can also be displayed in the Command Window. The Command History shows a compact listing of the commands that have been entered in the Command Window. Another important window that can be called up is the “Editor Window,” where 29 30 Chapter 3 Figure 3.1 The MATLAB Version 7.1 desktop screen shows the Current Directory window (upper left), the Command History window (lower left), the Editor (“docked” in the upper right of the Desktop by clicking the downcurving arrow in the Editor window), and the Command Window (lower right). code can be entered and saved as a file (an Mfile with extension “.m”). Two different types of Mfiles can be created: a script or a function. A script is a collection of code lines that can be executed as a program. When a script is executed, numerical output shows up in the Command Window. Function Mfiles can be called from the Command Window or from a script to do a particular task and return a result. The examples presented here are described as being entered in the Editor Window, which allows your script to be saved in an Mfile; but, the code can also be entered directly and executed in the Command Window. The first thing to do is make a new folder for your work. On the toolbar above the Current Directory, click on the “New Folder” icon and a folder should appear with a cursor positioned to enter the folder name. Type in a name, such as “Fourier Optics,” and click anywhere on the folder line or hit Enter. Then double click on this folder in the Current Directory window to open it as your current storage area. The first programming step will be to create a rectangle (rect) function. Once created, this piece of code can be called from a script to generate a vector containing a sampled rectangle function. Call up the Editor Window by clicking on the New Mfile icon (a documentlike icon), which is on the far left of the MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 31 MATLAB main Desktop toolbar. A cursor next to the line number 1 should appear in the window. Type in the following: 1 2 3 4 5 function[out]=rect(x); % rectangle function out=abs(x)<=1/2; end To do this, simply type out each line, and hit the Enter key at the end of each line. MATLAB uses a simple text format, so special function keys are not required for entering the code. The “%” character indicates that the text that follows is a comment and not to be interpreted for execution. Some of the text in the display will be colored. Comment text is displayed in green. As always, documentation of your code is important and a few comments are included in the examples presented throughout this book; but, in general, they are kept to a minimum to keep the presentation concise. Comments can also be added on the same line following the code. Bluecolored text indicates builtin MATLAB functions. The function command sets up this Mfile to be a function named “rect” with an input vector “x” and an output vector “out.” MATLAB is most proficient with vector and matrix operations, so the code tends to involve variables that are vectors and arrays as opposed to single parameters. The code in line 4 uses a vectorized logical test feature in MATLAB. The “abs” command takes the absolute value of each x element and “<=” applies the “less than or equal” test to each element. If the test is “true” (less than or equal to 1/2), then 1 is returned for that element. If the test is “false,” a 0 is returned. The vector out has the same number of elements as x and will contain the sampled rectangle function. The semicolon at the end of a line suppresses the output in the Command Window as the code is executing. In the rect code, the samples at the edges of the rectangle are not allowed to take on a value of 1/2 as in the analytic rectangle function definition (see Section 1.5). Doing so can be interpreted as a “slope” at the edge rather than a sharp transition, so we choose to work with a single point transition. The coding we use for rect has the characteristic that the full width of the function is always created with an odd number of samples. If a situation comes about where a rectangle with an even number of samples is requested, one more sample than expected will be returned. Some of the examples in the book use a width or radius value that is not a round number—this is simply to be consistent with the odd number of samples that will return from calling rect. When you are done typing the code, click on “Save” (floppy disk icon) on the Editor Window toolbar and save the function with the name “rect.” It is recommended that the file name be the same as the function name. An Mfile with the name “rect.m” should appear in the Current Directory window. Close 32 Chapter 3 the Editor Window. The function “rect” can now be called by any other script that is operating within the “Fourier Optics” folder. 3.2 Creating Vectors Suppose we want to program a simulation of a 1D spatial rectangle function with a half width of 0.055 m (full width of 0.11 m). This rectangle function will be created in a vector that corresponds to a physical length (side length L) of 2 m. If the number of samples in the vector is 200, then the sample interval x (or “dx” in the code) is 0.01 m. Open a new Editor Window by going to the main MATLAB desktop toolbar and clicking on New Mfile. In the new window, enter the following lines: 1 2 3 4 5 6 % fft_example  Chapter 3 fft example w=0.055; L=2; M=200; dx=L/M; %rectangle halfwidth (m) %vector side length (m) %number of samples %sample interval (m) Here, the parameters for the rectangle physical size, the vector side length, the number of samples, and sample interval are defined. Add the following code: 7 8 x=L/2:dx:L/2dx; f=rect(x/(2*w)); %coordinate vector %signal vector The “*” indicates multiplication. In this code, the coordinate vector x is created with values ranging from 1 (= –L/2) to 0.99 (= L/2dx) in steps of 0.01 (dx). The colons separate the range and step size. The subtraction of one step off the high limit results in a vector of 200 samples. The signal vector “f” is created with the help of the previously defined rect function and the use of the coordinate vector x. Since 2*w is a scalar, the interpretation is that each sample in x is divided by 2*w and the resulting vector is input to the rect function. To display what you have created in a plot of the f values against the x values, add the following lines: 9 10 figure(1) plot(x,f); %plot f vs x The “figure(1)” command opens a window labeled “Figure 1” and “plot” makes a plot. Before executing this script, you need to save it to a disk. Click on Save on the Editor Window toolbar and save this Mfile with the name fft_example. Don’t use any spaces in the name, as that can be interpreted as a request to link to another function; to execute, click on Run in the Editor Window toolbar (a triangle icon or for older MATLAB versions an arrow beside a document icon). If there are no errors in your code, the plot of f should appear MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 33 1 0.8 0.6 0.4 0.2 0 1 0.5 0 0.5 1 Figure 3.2 A rect function example plot. in the Figure 1 window. If something is wrong with the code, you will see a message in the Command Window. The plot should look like Fig. 3.2. There are many options for adjusting the plot display, which include changing the limits on the axes, adding labels, and changing graph lines and symbols. One approach is to use the “Edit Plot” (arrow icon) editor on the toolbar that appears at the top of each figure window. In this tutorial we will not go over the use of this editor and the reader is referred to the MATLAB “Help” (the “?” icon on the MATLAB Desktop window) to learn more. The plot editor is a way to interactively work with the display. However, the changes made to the plot are not recorded in the code; so, the next time the script is executed, the plot will be displayed as it was before editing. Commands can be included in the script to adjust the plot. For our example, go back and edit the last plot command and add “axis” and “xlabel” command lines: 11 12 13 14 figure(2) plot(x,f,'o'); %plot f vs x axis([0.2 0.2 0 1.5]); xlabel('x(m)'); The added option in the plot command ('o') places the marker o at each sample in the plot. The axis command sets limits on the values plotted on x and y axes and the xlabel command labels the xaxis. The figure(2) command will open a second window for this new plot. Without this command, the new plot would overwrite the first plot in the Figure 1 window. Note that the plotmodifying statements come after the plot call. Click on Run to execute the code (this also automatically saves the latest version of the code) and the plot appearing in the Figure 2 window (Fig. 3.3) more clearly shows the rectangle function. Note the distance across the top of the 34 Chapter 3 rectangle corresponds to 0.1 m (11 samples) whereas the base of the rectangle corresponds to 0.12 m (13 samples). The correct interpretation for this sampled rectangle function is a width of 0.11 m, which lies between the two measures. In fact, because the rect function coding generates an odd number of samples, if a half width of 0.05 m were used rather than 0.055 m, the resulting signal vector would be identical to that shown in Fig. 3.3, which means there would be a small disparity in the intended width and the digital representation. 3.3 Shift for FFT Disregarding the coordinate vector x for a moment, in terms of the 200 samples in f, a rectangle function with a width of 11 samples has been created that is centered in the middle of the vector f. This can be visualized by adding the following code: 15 16 17 18 figure(3) plot(f,'o'); axis([80 120 0 1.5]); xlabel('index'); With the x argument removed from the plot command, the vector f is plotted as a function of the vector index values. The resulting plot (Fig. 3.4) shows the center sample of the rectangle function at index position 101. Centering the rectangle in the middle of the vector allows for easy viewing, but as described in Section 2.4, the FFT algorithm expects the zerocoordinate value to be in the first index location. To shift the values for the FFT operation, the “fftshift” command can be applied, which takes the samples of the first half of the vector and swaps them with the samples in the second half. Add the following the code to the program and run the script to get the plot in Fig. 3.5. 1.5 1 0.5 0 0.2 0.1 0 x(m) 0.1 0.2 Figure 3.3 Expanded view of sampled rect function. MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 35 1.5 1 0.5 0 80 90 100 index 110 120 Figure 3.4 Sampled rect function versus index value. 19 20 21 22 23 f0=fftshift(f); %shift f figure(4) plot(f0); axis([0 200 0 1.5]); xlabel('index'); A close look at Fig 3.5 should convince you that the center sample of the rectangle is at index location 1, and five samples are found to the “right” and the remaining five to the “left,” where the left group is placed at the end of the vector. 1.5 1 0.5 0 0 50 100 index 150 200 Figure 3.5 Shifted rect function versus index value. 36 Chapter 3 A few comments about shifting for the FFT: first, it is more straightforward to use an even number of samples M for the vector and center the function of interest at an index of M/2+1. The fftshift function will then shift the M/2+1 sample to index 1. In our example M = 200 and the definition of the x vector causes the rectangle function to be centered at the M/2+1 position (101). More care needs to be taken when shifting vectors that contain an odd number of samples. For example, if the fftshift function is used with an odd number of samples, then the ifftshift function should be used to undo the shift. For an even number of samples, the fftshift function works both forward and backward. A second comment is that without the shift operation the FFT algorithm generates a transform for a function that is translated from the zero position, which means a linear phase term will be present in the result (shift theorem!). 3.4 Computing the FFT and Displaying Results To compute the FFT of the vector “f0,” add the following below the last piece of code: 24 25 26 27 28 29 30 31 32 33 F0=fft(f0)*dx; %FFT and scale figure(5) plot(abs(F0)); %plot magnitude title('magnitude'); xlabel('index'); figure(6) plot(angle(F0)); title('phase'); xlabel('index'); %plot phase Here, the 1D FFT algorithm in MATLAB is used. A capital letter is used for the frequency domain vector. Multiplying the result by the sample spacing dx is necessary to correctly approximate the analytic Fourier transform integral. Since each sample in F0 contains two pieces of information (the real part and the imaginary part); or alternatively, the magnitude and the phase, two plots can be used to display this result. The “abs” command takes the absolute value (magnitude) of the samples in F0 and the “angle” command extracts the phase of the samples, with output values ranging from − to . Plot titles and xaxis labels have been included in this code. Run the script and the plots should look like those in Fig. 3.6. MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 37 Figure 3.6 (a) Magnitude and (b) phase plots for FFT result versus index. The sinc function nature of the magnitude result in Fig. 3.6 is obvious where the peak of the sinc is centered at index 1. However, note that the samples are not located exactly where the “zeros” of the magnitude would occur. Thus, the valleys shown in the curve do not necessarily appear to reach zero. The phase plot may look curious, but there are essentially three phase values in the plot: 0, −, and . However, − and represent the same value in a modulo 2 format, so the sharp transitions in the plot between − and are not particularly important, and the phase is effectively constant in these places. These sharp 2 transitions occur because of slight numerical differences between samples. With this interpretation, you should understand that the important phase transition is from 0 to (or −), which occurs about every 18 points. Comparing the phase plot with the magnitude plot, it is apparent that the transition occurs at the magnitude zeros. Furthermore, a phase value of is equivalent to placing a minus sign on the magnitude value. Combining all of this information, the magnitude and phase plots of Fig. 3.6 represent a realvalued sinc function where the values in the main lobe are positive, the values in the first lobe are negative, the second lobe values are positive, and so on. In this case the realvalued sinc function could have simply been displayed on one plot; but, in general, Fourier transform results are complex. Once again for display reasons, it is helpful to center the FFT result in the vector. In addition, the spatial frequency coordinates need to be determined. Add the following to your script (cutting and pasting from the earlier code can help accomplish this quickly): 34 35 36 37 38 39 F=fftshift(F0); %center F fx=1/(2*dx):1/L:1/(2*dx)(1/L); %freq cords figure(7) plot(fx,abs(F)); %plot magnitude title('magnitude'); 38 Chapter 3 40 41 42 43 44 45 xlabel('fx (cyc/m)'); figure(8) plot(fx,angle(F)); %plot phase title('phase'); xlabel('fx (cyc/m)'); The frequency coordinates in fx are created as defined in Section 2.4. Running the script generates the plots shown in Fig. 3.7, where the sampled Fourier magnitude and phase are centered, and the frequency axis is scaled appropriately. 3.5 Comparison with Analytic Results When possible, it is good practice to test a new piece of code by applying simple inputs where the result can be compared with an analytic result. This helps diagnose problems and lets you build on previous code with confidence. For the example given in Section 3.4, an analytic result is easily found. For the expression x f ( x) rect , 2w the Fourier transform is F ( f X ) 2w sinc2wf X . To compare the discrete result with the analytic result, add the following code to the “fft_example” script: 46 47 48 49 50 51 52 53 54 55 56 57 58 F_an=2*w*sinc(2*w*fx); %analytic result figure(9) plot(fx,abs(F),fx,abs(F_an),':'); %plot magnitude title('magnitude') legend('discrete','analytic') xlabel('fx (cyc/m)') figure(10) plot(fx,angle(F),fx,angle(F_an),':'); %plot phase title('phase') legend('discrete','analytic') xlabel('fx (cyc/m)') Here, the frequency coordinates fx are used as input to the analytic solution to create the vector F_an. Note that the MATLAB builtin sinc function is used since it conforms to our definition (see Table 1.2). The plot command MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 39 Figure 3.7 (a) Magnitude and (b) phase plots for centered FFT result. arguments are set up to graph the discrete and analytic results on the same plot, and the legends are added to the display. The resulting plots (Fig. 3.8) compare the FFT and analytic results. The magnitude results are nearly identical, but the FFT result has slightly higher values than the analytic curve at the edges. This is a result of the periodic extension property of the FFT (Section 2.5). The phase plots differ only in some transitions between and − for the digital result, which are of no consequence. 3.6 Convolution Example A convolution can be performed using the FFT and applying the Fourier convolution theorem. The example presented here involves the convolution of two Gaussian functions of different widths. In the Editor Window, select New Mfile and enter the following code that defines and generates sample values for the two functions fa and fb: magnitude phase 4 discrete 0.1 discrete analy tic 0.06 0 0.04 2 0.02 0 50 analy tic 2 0.08 0 fx (cyc/m) (a) 50 4 50 0 fx (cyc/m) (b) Figure 3.8 Comparison of (a) magnitude and (b) phase of FFT and analytic results. 50 40 Chapter 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 % conv_example  Convolution: two Gaussian functions wa=0.3; wb=0.2; L=2; M=200; dx=L/M; %Gaussian 1 width [exp(pi) radius](m) %Gaussian 2 width [exp(pi) radius](m) %side length (meters) %number of samples %sample interval (m) x=[L/2:dx:L/2dx]; fa=exp(pi*(x.^2)/wa^2); fb=exp(pi*(x.^2)/wb^2); %x coordinates %Gaussian a %Gaussian b figure(1) plot(x,fa,x,fb,''); title ('functions'); xlabel('x (m)'); Select the name “conv_example” for this Mfile. Note the command for squaring each value in the x vector requires a period before the “^” symbol. Without the period, MATLAB will attempt a vector rather than a single element operation. Running this code produces the plot shown in Fig. 3.9(a). Now add the following to compute and plot the convolution result: 16 17 18 19 20 21 22 23 24 Fa=fft(fa); Fb=fft(fb); F0=Fa.*Fb; f0=ifft(F0)*dx; f=fftshift(f0); %transform fa %transform fb %multiply pointwise %inverse transform and scale %center result figure(2) plot(x,f); title('convolution') xlabel('x (m)') functions convolution 1 0.15 0.8 0.6 0.1 0.4 0.05 0.2 0 1 0.5 0 x (m) (a) 0.5 1 0 1 0.5 0 x (m) 0.5 (b) Figure 3.9 (a) Two Gaussian functions and their (b) convolution result. 1 MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 41 The FFTs of the Gaussian vectors are first computed. The functions fa and fb need not be shifted prior to the FFTs since the convolution only depends on their relative positions. The transform results are next multiplied “pointwise” orelement by element. This operation is indicated with a period placed before the product operator (.*). Without the period a vector/matrix multiplication would be attempted. Running the code produces the result shown in Fig. 3.9(b). As introduced in Section 2.6, a critical consideration for a convolution computed in this way is the periodic extension property of the FFT. The criterion is that the sum of the function’s supports should be less than the vector length. From Fig. 3.9(a), the estimated support of the significant values in the first Gaussian is Da 0.9 m and for the second Db 0.6, thus the sum of these two is less than the vector length L = 2 m. Exercise 3.2 gives you a chance to check if the curve in Fig. 3.9(b) is indeed a good representation of the analytic convolution. In this book convolutions are coded directly by applying the convolution theorem, but MATLAB has the builtin function conv and for two dimensions, conv2. To strictly avoid artifacts due to the periodic convolution, these functions “zeropad” the initial vectors, placing them in doublesized vectors or arrays of zeros, and then performing the convolution. For computational speed and efficiency we tend to work with fixed array sizes and pay heed to the support of the signals, as was done in the example above. 3.7 Two Dimensions Physical optics problems typically involve at least two spatial dimensions. Start a New Mfile, named fft2_example, and enter the following code to generate a twodimensional (2D), sampled rectangle function: 1 2 3 4 5 6 7 8 9 10 11 % fft2_example  2D FFT example wx=0.1; %rect x halfwidth (m) wy=0.05; %rect y halfwidth (m) L=2; %side length x&y (m) M=200; %samples/side length dx=L/M; %sample interval (m) x=L/2:dx:L/2dx; %x coordinates y=x; %y coordinates [X,Y]=meshgrid(x,y); %X and Y grid coords g=rect(X/(2*wx)).*rect(Y/(2*wy)); %signal In this example, a 0.2 0.1 m rectangle is modeled in a 200 200 element array corresponding to a physical size of 2 2 m. The side length L is the physical length along one edge of the array where it is assumed that the x and y dimensions are the same. Two identical sample coordinate vectors x and y are defined for the two dimensions. The “meshgrid” command generates the coordinate arrays X and Y (capital letters) where the rows of X are copies of the 42 Chapter 3 vector x and the columns of Y are copies of y. X and Y are used to produce the sampled version of the 2D rect function in the array g. By using X and Y, the coded version of g appears much like an analytic expression. An image is a common way to visualize optics simulation results. To display the contents of g as an image, add the following: 12 13 14 15 16 17 figure(1) imagesc(x,y,g); %image display colormap('gray'); %linear gray display map axis square; %square figure axis xy %y positive up xlabel('x (m)'); ylabel('y (m)'); The command “imagesc” scales the array data to the full display range and presents the image and “colormap” provides a pseudocoloring of the image. The “gray” scale is used here, but check Help for other colormaps. The “axis square” command presents the figure border as a square rather than the default rectangular shape. This is helpful in this case where the side length is the same in the x and y directions. The first row of a conventional image file corresponds to the top of the picture. MATLAB image display functions assign yaxis values to increase from top to bottom. The axis xy command arranges the y axis to be displayed with increasing values from bottom to top. Running the script produces the image shown in Fig. 3.10(a). It is also often helpful to view 1D “slices,” or profiles, of 2D results. Add the following code to generate a 1D profile of the xaxis through the center of the array [Fig. 3.10(b)], where the syntax (M/2+1,:)selects the M/2+1 row: 18 19 figure(2) plot(x,g(M/2+1,:)); %x slice profile 1 0.8 y (m) 0.5 0.6 0 0.4 0.5 1 1 0.2 0.5 0 x (m) (a) 0.5 0 1 0.5 0 x (m) 0.5 (b) Figure 3.10 (a) Image and (b) x profile of 2D rectangle example. 1 MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 20 21 43 xlabel('x (m)'); axis([1,1,0,1.2]); To perform a 2D FFT of the rectangle function, add the following: 22 23 24 25 26 27 g0=fftshift(g); G0=fft2(g0)*dx^2; G=fftshift(G0); %shift %2D fft and dxdy scaling %center fx=1/(2*dx):1/L:1/(2*dx)(1/L);%x freq coords fy=fx; %y freq coords The fftshift command operates on a 2D array in the manner described in Section 2.4 and the fft2 command is needed for a 2D transform. Enter the following code to display the magnitude of the transform results as a surface plot along with a profile slice through the center (Fig. 3.11): 28 29 30 31 32 33 34 35 36 37 38 figure(3) surf(fx,fy,abs(G)) %display transform magnitude camlight left; lighting phong colormap('gray') shading interp ylabel('fy (cyc/m)'); xlabel('fx (cyc/m)'); figure(4) plot(fx,abs(G(M/2+1,:))); %fx slice profile title('magnitude'); xlabel('fx (cyc/m)'); The “surf” plotting command is introduced to illustrate another display method for a 2D array. The lighting, shading, and colormap commands can be used to change the display. The magnitude could also be displayed in other ways, such as an image or a contour plot. The phase of the result is not shown here, but it could also be displayed in a variety of ways. 3.8 Miscellaneous Hints A few other helpful MATLAB hints: Clear all: MATLAB retains the values of variables and arrays that are created when a script is executed. Typing a variable name in the Command Window and hitting enter displays the current value, which can be useful for analyzing your code. However, sometimes this gets confusing as code is being edited. The clear all command clears the variable memory. Close all: This command will close all the current figure windows. Display complex result: Sometimes a plot or image will not display and the error warning says the input is complex. After operations, like taking 44 Chapter 3 (a) magnitude 0.025 0.02 0.015 0.01 0.005 0 50 0 fx (cyc/m) 50 (b) Figure 3.11 (a) Surface plot display and (b) profile of 2D FFT magnitude result. the FFT, machine precision can leave some residual complex component in a vector or array that should be real. The real or abs functions can be applied to the array to allow the plot to display. Image contrast: The display function imagesc is used extensively throughout this book. For printing purposes the images are displayed in grayscale. However, it is easier to see lowvalue features with different colormaps—so try some other maps. To stretch the contrast of a grayscale image to more easily see dim features, a quick trick is to display the nth root of the image values. For example, nthroot(g,3) takes the third root of g. The higher the root, the more the contrast is stretched. Just be sure to remember that you are looking at a peakscaled, contraststretched image. Vector operations: Pay special attention to vectorized operations indicated with the period—for example, an elementbyelement multiplication indicated by A.*B. The most common programming error is to forget the period when a vectorized operation is needed. This mistake generally does MATLAB Programming of Functions, Vectors, Arrays, and Fourier Transforms 45 not create an error message, but the result will be unexpected. Make a habit of checking the vector operations when things are going wrong. 2D transform: Be sure and use fft2 and ifft2 for two dimensions—not fft and ifft. The 1D functions operating on the 2D array tend to give a repeated 1D result (it looks “striped”)—which is not good! 3.9 Exercises 3.1 Triangle function Mfile: (a) Create a triangle function in an Mfile. Try some lines like: T=1abs(x); ask=abs(x)<=1; out=T.*mask; (b) In a script, create a sampled triangle function using the following specifications: triangle base half width = 0.1 m, vector length = 2 m, and number of samples M = 200. (c) Plot the sampled function. (d) Compute the FFT. (e) Find the analytic Fourier transform of the function in (b). (f) Plot the FFT and analytic Fourier transform results together (both magnitude and phase). 3.2 Code the example for the convolution of the Gaussian functions presented in Section 3.6. Find the analytic convolution of these functions and compare this result with the discrete result in a plot. 3.3 Circle function Mfile: (a) Create a circle function in an Mfile. (b) Generate a sampled circle function in a 2D array with the following parameters: circle radius = 0.015 m, array side length = 0.2 m, and number of samples (one dimension) M = 200. (c) Display the sampled function as an image. (d) Take the FFT of the array and display the magnitude of the transform in surface and profile plots. Chapter 4 Scalar Diffraction and Propagation Solutions Perhaps the most fundamental task associated with Fourier optics is describing the evolution of an optical field as it propagates from one location to another. The phenomenon of diffraction underlies the behavior of propagating waves. Extensive theory developed for diffraction provides the basis for modeling optical propagation on the computer. This chapter is essentially a summary of scalar diffraction theory with a listing of the expressions commonly used today to describe