Efficient matrix approach to optical wave propagation and Linear Canonical Transforms

Sami A. Shakir; David L. Fried; Edwin A. Pease; Terry J. Brennan; Thomas M. Dolash

doi:10.1364/OE.23.026853

1. Introduction

The problem of scalar optical wave propagation is generally approached in one of three forms. One is the differential form, as represented by the Helmholtz wave equation [1–6]. A second form is the integral representation, which has been formulated in different forms, such as the Rayleigh-Sommerfeld formulation [1–3] and the Fresnel-Kirchhoff formulation. These integral formulations are based on a Green's function solution of the Helmholtz wave equation [1,2]. This solution is also called the impulse response and leads to a convolution formulation of the wave propagation problem. This convolution formulation displays a rich set of properties associated with general linear systems [1,2]. A third approach to wave propagation is in terms of the direct solution of Maxwell's equations [3,4]. This formulation is most convenient for problems dealing with systems having certain boundary conditions that may lead to wave confinement by the boundaries. These systems are typically called resonators and waveguides. In this work, we only focus on the solution of the scalar Fresnel diffraction integral, which is an approximation of the integral formulation of the Helmholtz wave equation, under the assumption of propagation within a narrow angle. This is commonly called the paraxial approximation. The results discussed here are actually more general than the Fresnel diffraction integral and are representative of a class of linear transforms called, the Linear Canonical Transforms (LCT) [7–9].

In what follows, we will focus mostly on the computational aspects of this efficient matrix algorithm (EMA). This matrix formalism can also form a useful analytical formalism reminiscent of the quantum matrix formalism. We will show that if one uses one of the efficient matrix multiplication algorithms, which is standard on some computational platforms such as MATLAB and the BLAS C + + library, then EMA can be competitive with wave propagation methods based on the Fast Fourier Transform (FFT) [10,11]. Due to the integral nature of EMA, it has transparent boundary conditions, which negates the need for very wide computational windows. Furthermore, the inherent flexibility of choosing a different number of points for the input and output planes can give EMA a significant speed advantage over FFT-based methods for cases where only a small number of solution points are of interest in the output plane.

2. An efficient matrix approach (EMA)

It is well known [5] that the free-space Fresnel diffraction integral can be cast in a form, called generalized Huygens-Fresnel diffraction integral [5], that applies to any composite paraxial optical system as long as the individual components, such as lenses and Gaussian apertures [5], can be represented by a paraxial optical ABCD-matrix [2,3,5]. This integral can be expressed in a slightly different form that represents a broad class of linear transforms called Linear Canonical Transform (LCT) [6–9],

U (x_{2}, y_{2}) = \frac{1}{i b} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} u (x_{1}, y_{1}) e^{i \frac{π}{b} [(a x_{1}^{2} - 2 x_{1} x_{2} + d x_{2}^{2}) + (a y_{1}^{2} - 2 y_{1} y_{2} + d y_{2}^{2})]} d x_{1} d y_{1}

where the transverse coordinates in the input (or source) plane at lateral location z = 0, and the output (or observation) plane, at lateral location z, are

{x_{1}, y_{1}}

and

{x_{2}, y_{2}}

, respectively, as shown in Fig. (1). The field in the source (input) plane is

u (x_{1}, y_{1})

and in the output plane is

U (x_{2}, y_{2})

, and the parameters {a, b, and d}, are the elements of a 2x2

" a b c d "

matrix that satisfies the relationship:

a d - b c = 1

. A specific choice of these elements reduces Eq. (1) to one of the well known linear integral transforms such as the Fourier, Fractional Fourier [9], free-space Fresnel diffraction integral, ABCD-paraxial optical system Fresnel diffraction integral, Laplace transform, Fractional Laplace transform, and Gauss-Weirestrass transform [6–9]. The particular configurations of the elements of the abcd matrix and their corresponding integral transform are summarized in Table 1.

Fig. 1 A graphical representation of a system being influenced by a Linear Canonical Transform. A source field, u(x₁,y₁) is transformed into an observation field U(x₂,y₂), in the observation plane a distance z from the source plane.

Download Full Size | PDF

Table 1. The abcd matrix configuration and corresponding Linear Canonic Transforms

View Table

For both numerical computation and analysis purposes, we construct a representation of Eq. (1) in terms of a matrix transformation as,

U = H^{T} u H

We will outline the derivation of Eq. (2) shortly. The corresponding inverse transform can be expressed as,

u = H^{*} U H^{*}^{T}

where “*” represent the complex conjugate, “T” represent the matrix transpose operation, and u represents the discrete form of the field source which is a matrix of dimensions

N_{1} x N_{1}

while U is the discrete form of the output fields and is a matrix of dimensions

N_{2} x N_{2}

, as depicted in Fig. (1). Without the loss of generality, we have assumed that the x and y transverse grid points of a given plane are similar. It is useful to note that the grid spacing do not have to be uniform, an observation that may be useful for problems where one wants to focus on certain regions but not others. The matrix H is the discrete form of the impulse response function (IR) of the linear system. It is also called the quadratic phase kernel to highlight its oscillatory nature. The impulse response matrix H is of dimension

N_{1} x N_{2}

and is defined by,

H_{j 1, j 2} = W_{j1} \frac{1}{\sqrt{i b}} e^{i \frac{π}{b} (a x_{j 1}^{2} - 2 x_{j 1} x_{j 2} + d x_{j 2}^{2})}

where the weight coefficient W_j1 represent the coefficients employed in discrete integration algorithms. For an equally spaced mesh of grid points with grid separation of Δ in both the x and y direction; the extended Trapezoidal Rule [10] algorithm requires,

W_{j} = (0.5, 1, 1, ....0.5) Δ

, and for the extended Simpson's Rule [10] algorithm, we have,

W_{j} = (1, 4, 2, 4, 2....4, 1) Δ / 3

. In effect, Eq. (2) transforms a discrete field distribution of size

N_{1} x N_{1}

at a source plane at lateral position z = 0 into a discrete field distribution of size

N_{2} x N_{2}

at an observation plane at lateral position z. The proof of Eq. (2) is based on the recognition that the impulse response function of the Linear Canonical Transform of Eq. (1) is separable in spatial coordinates {x, y}. Thus, for the free-space Fresnel diffraction integral,

\begin{array}{l} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} u (x_{1}, y_{1}) e^{i \frac{π}{λ z} [{(x_{n} - x_{1})}^{2} + {(y_{m} - y_{1})}^{2}]} d x_{1} d y_{1} \approx \sum_{j = 1}^{N} W_{j} e^{i \frac{π}{λ z} {(y_{m} - y_{j})}^{2}} \sum_{k = 1}^{N} W_{k} e^{i \frac{π}{λ z} {(x_{n} - x_{k})}^{2}} u_{j, k} \\ = \sum_{j = 1}^{N} \sum_{k = 1}^{N} H_{m, j} u_{j, k} H_{k, n} = {[H^{T} u H]}_{m, k} \end{array}

where

x_{k} = k Δ

, and

y_{j} = j Δ

; while k and j are integers,

j = 1, 2, ... N_{1}

and

k = 1, 2, ... N_{2}

. If one expands Eq. (2) using the definition of H in Eq. (4), then it becomes immediately clear the equivalence of Eqs. (5) and (2). The same method can be used to show the inverse form of Eq. (3). Note the freedom in choosing the number of grid points in either the source plane and/or the observation plane so that they don't have to be equal or uniform. We will call this approach, the Efficient Matrix Algorithm (EMA), and discuss the issues surrounding the implementation of this matrix formalism. In what follows, we focus on the free-space Fresnel diffraction integral as our main problem. The conclusions seem to be general and applicable to other transforms represented by Eq. (1), however, we have not yet studied the issues presented here, as applied to LCT systems in general.

3. The Mean-Kernel Method

The Fresnel diffraction integral is basically a convolution integral that cast the wave propagation problem as the linear superposition of the radiation field of point sources comprising an extended source. This is a mathematical expression of Huygens's principle, within the paraxial approximation. An alternative way is to view the source as a collection of infinite extent plane waves that are linearly superimposed at an observation plane. To extract the plane wave representation of the source, one has to resort to Fourier transforms to convert the source field from point sources to plane waves. The presence of a fast transform such as the Fast Fourier Transform (FFT) algorithm [10,11], helped make this approach widely used. However, transforming a spatial field distribution into a Fourier transform and back again, may introduce artifacts caused by aliasing [13,14], which is a byproduct of replacing the actual (non-periodic) problem by a periodic one. In practice, the effect of aliasing resembles reflection of the field from the boundaries. A popular solution is to introduce an artificial absorber around the boundaries, called “Guard Bands”, to reduce the field near the boundaries. This, hopefully, reduces the aliased field to a level that is negligible. The introduction of an artificial component may alter the problem sufficiently, that the results are not an accurate representation of the actual problem. Another computational issue is the possibility of impulse response function becoming highly oscillatory, which requires closely spaced grid points to sample the function properly, otherwise, high spatial frequency components might be missed and result in errors. These problems necessitate careful planning and monitoring of the propagation problem.

Since the EMA wave propagation approach does not rely on transforming the field into a Fourier transform, there is no requirement for the field to be zero at the boundaries. Sampling of the product of the source field and the impulse response function controls the accuracy level, as discussed in section-5. The impulse response function may pose a computational challenge when it becomes very oscillatory. This may take place when z becomes relatively small or/and the transverse observation position (x₂,y₂) is large, so as to result in a large quadratic phase factor of Eq. (4). This behavior is also commonly expressed in terms of the Fresnel Number, $F_{N}$ , which is a dimensionless number defined [1,2,5] as $F_{N} = a^{2} / λ z$ , where a is the source beam’s radius, and λ is the wavelength and z is the propagation distance. The Fresnel number establishes a coarse criterion to define the near field and far field cases. Essentially, if the Fresnel number is small - less than 1 - the beam is said to be in the “far field” [13,14]. If the Fresnel number is larger than 1, the beam is said to be “near field”. In effect, the phase term of the impulse response function in the Fresnel diffraction integral, will oscillate very rapidly when the Fresnel Number is large. This means that the sampling grid size must be fine enough to resolve the impulse response function accurately. This may cause the arrays sizes to be too large for acceptable computational speed. Fortunately, the nature of EMA wave propagation offers a simple and effective solution to this problem. Instead of sampling the source field and the impulse response function at integration grid points, we replace each integration segment by the product of an effective average field over the segment, and retain the impulse response function to be integrated over a continuous distribution of point sources over that segment. In effect, we are replacing a set of point sources by a “mean” source cell whose length is the same as the grid separation, Δ. A similar idea was used by Milonni and Paxton [15] to model carbon monoxide laser modes in an unstable resonator. To quantify this approach, we note that the contribution to the field at an observation point (x₂,y₂) is separable into an x-component and y-component factor for each two-dimensional cell, so that $Δ U = Δ U_{x} Δ U_{y}$ . The x-component is

Δ U_{x} = \sqrt{1 / i λ z} e^{i π z / λ} u_{o} (x_{1}, y_{1}) \int_{x_{1} - 0.5 Δ}^{x_{1} + 0.5 Δ} e^{i \frac{π}{λ z} {(x - x_{2})}^{2}} d x

and a similar equation for the y-contribution as,

Δ U_{y} = \sqrt{1 / i λ z} e^{i π z / λ} u_{o} (x_{1}, y_{1}) \int_{y_{1} - 0.5 Δ}^{y_{1} + 0.5 Δ} e^{i \frac{π}{λ z} {(y - y_{2})}^{2}} d y

The impulse response matrix of Eq. (4) becomes,

H (x_{1}, x_{2}) = \sqrt{\frac{e^{i k z}}{i λ z}} \int_{x_{1}_{} - 0.5 Δ}^{x_{1}_{} + 0.5 Δ} e^{i \frac{π}{λ z} {(x_{2}_{} - x)}^{2}} d x

Introducing a variable,

ξ = x - x_{1}

, and recalling the Fresnel approximation which implies that the grid size

Δ < < \sqrt{λ z}

, it becomes evident that Eq. (8) essentially represents the Fourier transform of a rect-function [1–3]. Equation (8) to a good approximation is,

H (x_{1}, x_{2}) = \frac{e^{i π z / λ} Δ}{\sqrt{i λ z}} e^{\frac{i π}{λ z} {(x_{2} - x_{1})}^{2}} [\frac{\sin (π (x_{2} - x_{1}) Δ / λ z}{π (x_{2} - x_{1}) Δ / λ z}]

With a similar formula for Eq. (7) but with x₁ and x₂ replaced with y₁ and y₂, respectively. The “mean-kernel” representation culminating in Eq. (9) is more accurate than simple sampling of the field and impulse response kernel at integration grid points. Moreover, within the Fresnel approximation removes the difficulty of large Fresnel numbers when z is small. Notice that for small Fresnel numbers the sinc-function [1–3] of Eq. (9) reduces to unity.

In closing this section, we summarize by stating that the EMA propagation approach is simply to compute the impulse response matrix H using Eq. (9), then apply the simple matrix product involving the field u as, $H^{T} u H .$

4. Applications and examples

The implicit periodic nature of FFT based methods, leads to the requirement that the field be vanishingly small at the boundaries in the output plane. Otherwise, aliasing can cause artifacts in the field. The fingerprint of these artifacts is usually in the form of interference patterns where they are not expected, which may cause substantial errors in the results. Taking a large computational window to insure vanishingly small fields at the boundaries is the common approach. However, this may not be practical since the number of grid points and resulting computational load may be too large. One widely used solution is to add an artificial absorber, called a “Guard Band” or a filter in the frequency space in order to reduce the field near the boundaries and thus, reduce the artificial contribution due to aliasing [13,14]. This approach requires careful design; otherwise it may alter the nature of the problem. For the EMA approach of Eq. (2), there is no such requirement since no Fourier transform and periodicity is invoked. The basic requirement is that the field in the source plane be sampled properly. In effect, EMA has transparent boundary conditions (TBC) [16]. Normally, the quadratic phase kernel needs to be sampled properly too. This constitutes a major challenge when the quadratic phase term, or Fresnel number, is large. To illustrate some aspects of EMA, we show an example of a square flat-top laser beam of 1 micron in wavelength, and an aperture half-width of a = 5 cm. The beam is propagated to a distance of z = 100 m in first case of Fig. (2a).We used N = 50 grid points for the source plane and observation plane. The figure shows the irradiance profile as predicted by analytical formulas for the Fresnel diffraction of a square aperture [1] (black solid line) and EMA's corresponding prediction (red points), as implemented in Eq. (9) for the impulse response matrix. In Fig. (2b) we illustrate the TBC nature of EMA by propagating a Gaussian beam of 2.5 cm amplitude radius (1/e) a distance of 10 km. The solid black line shows the analytical solution for a Gaussian beam propagated a distance of 10km [3], while the red dots are the EMA predictions. For comparison, we also show the solution using an FFT propagator (green dots) [13,14]. The observation window was chosen to highlight the TBC nature of EMA. Figure (2b) shows clearly that while EMA is free from aliasing, the FFT based solution is showing strong aliasing. The accuracy of EMA's approach is evident even for relatively large Fresnel numbers. The Fresnel number in Fig. (2a) was 25 and in Fig. (2b) was 6.25. The absence of aliasing in Fig. (2b) clearly shows the TBC [16] nature of EMA.

Fig. 2 (a) The irradiance distribution of a flat-top square laser beam with a half-width of 5 cm propagated a distance of 100 m. (b) A Gaussian beam of 2.5 cm radius propagated a distance of 10 km. In both cases the laser wavelength is one micron.

Download Full Size | PDF

The next example, shown in Fig. (3), demonstrates multiple steps propagation in the presence of an atmospheric phase screen. In this scenario, a Gaussian beam of 1.5 cm waist radius (1/e²) is propagated a distance of 10 km. We simulate the atmosphere by placing a phase screen midway, at z = 5km. The atmospheric coherence length, r_o, is assumed to be 5 cm. We employed the screen generation methods described by Schmidt and we also used the FFT propagation MATLAB codes described in the same reference [14]. We first propagated the Gaussian beam a distance of 5 km to where the atmospheric phase screen is located (Fig. 3a) and multiplied the propagated field by the phase screen phase exponential, and propagated the modified field a distance of 5 km to the observation plane (Fig. 3b). In both Figs. (3a) and (3b), we compare the EMA propagated field irradiance line-out along the x-axis (green solid line) to the corresponding field irradiance using the FFT propagation code (blue solid line) and for reference we show a Gaussian profile (red line), with an effective radius that is root-mean-square of the sum of the source beam’s radius and the square of the divergence radius caused by turbulence, i.e.

Fig. 3 Propagation of a Gaussian beam through the atmosphere. The 1-micron wavelength beam waist radius is 1.5 cm, is first propagated a distance of 5 km (a), then an atmospheric equivalent phase screen (r₀ = 5cm) is applied and the modified beam is propagated another 5 km (b). The two-dimensional beam at 10 km is shown as EMA calculated (c), and FFT calculated (d).

Download Full Size | PDF

r_{e f f e c t i v e} = \sqrt{a_{}^{2} + {(λ z / r_{o})}^{2}}

The two-dimensional irradiance distribution as calculated by EMA (c) and FFT propagation (d), are also shown in the figure. Unlike the FFT propagation, which requires that the computational window in the observation plane should be wide enough to avoid aliasing, the EMA calculates the field for a smaller observation window, without showing any trace of aliasing. Of course, the computational window in the intermediate plane has to be wide enough to encompass the entire significant field, for any method.

5. Accuracy and error bounds

One of the advantages of the EMA approach is the capability to quantify the truncation errors for a given choice of parameters, such as the grid size and number of sample points. This error can be quantified by noting that Eq. (2) is essentially based on numerical integration and the errors are basically truncation errors. This assumes that the integration limits are large enough to encompass all regions where the field is significant. Hence, if the weight factor (W) in Eq. (4) represent the Trapezoidal Rule, then the truncation error for 2D problems is given by [10],

E r r o r (x_{2}, y_{2}) \leq \frac{L^{2} Δ^{2}}{12 λ z} \max_{} | \frac{\partial^{2} f (x_{1}, y_{1}; x_{2}, y_{2})}{\partial^{2} x_{1}} + \frac{\partial^{2} f (x_{1}, y_{1}; x_{2}, y_{2})}{\partial^{2} y_{1}} |

Where

f (x_{1}, y_{1}; x_{2}, y_{2}) = u (x_{1}, y_{1}) e^{i \frac{π}{λ z} [{(x_{2} - x_{1})}^{2} + {(y_{2} - y_{1})}^{2}]}

, and Δ is the grid increment size, and L is the total integration length, with the limits of integration being

(\pm L / 2)

.

Equation (11) shows how the truncation error is proportional to the square of the length of the computation window (L) and also proportional with the grid/cell area, $Δ^{2}$ . Furthermore, the formula shows how the truncation error term depends on the rate of change of the product of the field and the impulse response function, and not the field alone. In addition to giving insight into the dependency of the error term, it also provides a way to choose parameters to keep the error within required tolerances. It is interesting that EMA’s accuracy improves with smaller computation window, L. This is the opposite of the trend with FFT-based propagation methods. Coupling this property with the TBC property of EMA may be a significant advantage, especially in propagation problems where one is only interested in the field in a certain region, such as that of a sensor or a receiver’s aperture. In choosing the computational width L, one has to make sure that L encompasses the entire radiation region contributing to the observation window. The requirement that the grid size be small enough to properly sample the impulse response function leads to the requirement [13], $L \leq λ z / 2 Δ$ . The equality represents what is called critical sampling [13].

The ability to estimate the upper bound on the computational error for a given set of parameters before actually performing the propagation is an advantage, especially for situations where a known level of accuracy needs to be maintained.

6. Flexibility and computational efficiency

Another important advantage of EMA is the flexibility in choosing the number and spacing of sampling points differently in different planes. While the source plane sampling is restricted by sampling theory requirements, there are no such requirements for the observation plane. In fact the number of grid points in the observation plane can be as small as a single point. This level of flexibility offers a tremendous speed advantage that can be as large as two or three orders of magnitude for problems where one is interested in a single point or a small region in the observation plane. This situation is very common in situations where one is only interested in metrics, such as the Strehl, and the Power-In-The-Bucket (PIB) estimates. This is in contrast to FFT based methods where the number of points must be the same for input and output planes.

FFT based methods derive their speed advantage over other methods from the fact that its computational complexity is of the order, $O (N \log_{2} (N))$ , where N is the linear number of grid points. Since matrix multiplication of two matrices, each of size $N x N$ , involves computational complexity of order, $O (N^{3})$ , then FFT based methods have a huge speed advantage. While this is true for 1D problems, the situation is more complex for higher dimensions, especially if one uses a “fast” matrix multiplication algorithm¹². In fact, matrix multiplication can be of order $O (N^{β})$ where, $2.35 < β < 3$ , depending on the algorithm used¹². Hence, for M-dimensional (M = 1, 2) LCTs, the number of multiply operations increases linearly with the dimension M, as, $O (M N^{β})$ , while it increases exponential for the FFT based methods, as $O (N^{M} \log_{2} (N^{M}))$ . To compare the speeds more accurately than just order of magnitude estimates, we applied the EMA method and an FFT based wave propagation on several problems, ranging from a simple tilted Gaussian beam to an array propagation in the atmosphere to compare the speed and accuracy of the two approaches. We used an FFT based propagator, such as the one in [14], which involves three FFT operations and one element by element multiplication of two arrays. The ratio of the time taken for the FFT method to that of EMA, we call the relative speed of EMA with respect to the FFT based method (> 1 means EMA is faster), and is very roughly given by

S p e e d \approx c \frac{1 + 3 \log_{2} (N^{2})}{N^{β - 2}}

Where, we assumed that the matrix multiplication algorithm is using the Strassen algorithm [12] which results in β = 2.8. The proportionality factor c of Eq. (12) depends on the propagation method used and the computer architecture. For best fit to timing data, we set the factor c to a value of 2. To generate the speed comparison curves in Fig. (4), we used MATLAB to run the wave propagation code. MATLAB was a convenient platform because it employs the Strassen algorithm for matrix multiplication, and uses optimized FFT routines for FFT calls that do not require the number of points to be integer exponents of 2. The computation was repeated 100 times on a four-core PC machine based on Windows 7 and Linux operating systems. The averages show similar results. For reference, we also show the expected speed for the so called “naive multiply” (β = 3), and the Coppersmith-Winograd algorithm (β = 2.4) [12]. We emphasize that the results shown in Fig. (4), are meant to be a rough indicator of trends for normal applications of EMA and FFT based wave propagation.

Fig. 4 Speed ratio for EMA over FFT wave-propagator. Solid circles represent actual simulation using MATLAB. Red line is the formula of Eq. (12) for β = 2.8. Blue broken line is for β = 2.4, and green line for β = 3.

Download Full Size | PDF

As the relative speed curves of Fig. (4) show, the EMA method is faster than the FFT based wave propagator for grid point numbers N < 300, but looses that advantage as N becomes larger than 300. However, if one employs the “mean kernel” method introduced in section-3, then one may use a smaller number of sampling points without sacrificing accuracy, resulting in a significant advantage in speed. Furthermore, as pointed out earlier, the flexibility in choosing the number of sampling points in the observation plane, can provide a significant speed advantage too. Finally, there is the prospect of parallel processing which may favor the matrix multiplication process over FFT algorithms. We have not studied this question extensively, but our preliminary application on GPU processors seem to speed up EMA more than FFT based wave propagation, as intuitively expected. The algorithms used in this context were optimized for FFT and matrix multiplication.

7. Conclusion

The Fresnel diffraction integral form of optical wave propagation and the more general Linear Canonical Transforms (LCT) are cast into a matrix transformation form. Taking advantage of recent efficient matrix multiply algorithms, this approach promises an efficient computational and analytical tool that is competitive with FFT based methods but offers better behavior because it supports transparent boundary conditions. The method also offers the advantage of freedom in choosing the number of sampling points and computational window sizes differently for the input and output planes. This flexibility makes the method significantly faster than FFT based propagators when only a single point, as in Strehl metrics, or a limited number of points, as in power-in-the-bucket metrics, are needed in the output observation plane. The ability to quantify the maximum computational error may also be significant for problems where a certain level of accuracy is required.

Acknowledgment

Tau Technologies acknowledges the partial support of Air Force Research Laboratory (AFRL) at Kirtland Air Force Base.

References and links

1. J. W. Goodman, Introduction to Fourier Optics, 3rd Edition (Roberts and Company, 2004).

2. M. Born and E. Wolf, Principles of Optics (Pergamon Press, 1969).

3. B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics, 2nd ed. (John Wiley& Sons, 2007).

4. D. Yevick and B. Hermansson, “Efficient beam propagation techniques,” IEEE J. Quantum Electron. 26(1), 109–112 (1990). [CrossRef]

5. A. E. Siegman, Lasers (University Science Books, 1986).

6. K. B. Wolf, Integral Transforms in Science and Engineering (Plenum Press, 1979).

7. M. Moshinsky and C. Quesne, “Linear canonical transformations and their unitary representations,” J. Math. Phys. 12(8), 1772–1783 (1971). [CrossRef]

8. B. M. Hennelly and J. T. Sheridan, “Fast numerical algorithm for the linear canonical transform,” J. Opt. Soc. Am. A 22(5), 928–937 (2005). [CrossRef] [PubMed]

9. H. M. Ozaktas and D. Mendlovic, “Fractional Fourier optics,” J. Opt. Soc. Am. A 12(4), 743–751 (1995). [CrossRef]

10. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientic Computing, 3rd ed. (Cambridge University Press, 2007).

11. E. O. Brigham, The Fast Fourier Transform, (Prentice Hall, 1974).

12. Wikipedia contributors, “Matrix multiplication,” Wikipedia, The Free Encyclopedia,http://en.wikipedia.org/w/index.php?title=Matrix_multiplication&oldid=600028683

13. D. G. Voelz, Computational Fourier Optics: A MATLAB Tutorial (SPIE Press, 2011).

14. J. D. Schmidt, Numerical Simulation of Optical Wave Propagation With Examples in MATLAB (SPIE Press Monograph, 2010).

15. P. W. Milonni and A. H. Paxton, “Model for the unstable-resonator carbon monoxide electric-discharge laser,” J. Appl. Phys. 49(3), 1012 (1978). [CrossRef]

16. G. R. Hadley, “Transparent boundary condition for beam propagation,” Opt. Lett. 16(9), 624–626 (1991). [CrossRef] [PubMed]

Efficient matrix approach to optical wave propagation and Linear Canonical Transforms

Abstract

1. Introduction

2. An efficient matrix approach (EMA)

3. The Mean-Kernel Method

4. Applications and examples

5. Accuracy and error bounds

6. Flexibility and computational efficiency

7. Conclusion

Acknowledgment

References and links

Cited By

Figures (4)

Tables (1)

Equations (12)

Optics Express