Efficient computation of quadratic-phase integrals in optics

Abstract

We present a fast N log N time algorithm for computing quadratic-phase integrals. This three-parameter class of integrals models propagation in free space in the Fresnel approximation, passage through thin lenses, and propagation in quadratic graded-index media as well as any combination of any number of these and is therefore of importance in optics. By carefully managing the sampling rate, one need not choose N much larger than the space-bandwidth product of the signals, despite the highly oscillatory integral kernel. The only deviation from exactness arises from the approximation of a continuous Fourier transform with the discrete Fourier transform. Thus the algorithm computes quadratic-phase integrals with a performance similar to that of the fast-Fourier-transform algorithm in computing the Fourier transform, in terms of both speed and accuracy. A quadratic-phase system is a unitary system whose output g͑u͒ is related to its input f͑u͒ through a quadratic-phase integral:where ␣, ␤, and ␥ are real parameters. This relationship is also known by other names, including "linear canonical transforms." [1][2][3][4] This input-output relation can model a broad class of optical systems, including thin lenses, sections of free space in the Fresnel approximation, sections of quadratic graded-index media, and arbitrary concatenations of any number of these. 2,3Computation of the Fresnel diffraction integral, which is a special case of Eq. (1) with ␣ = ␤ = ␥, has received the greatest attention because it describes the propagation of light in free space (see Refs. 5 and 6 and the references therein). As the input-output relationship represented by the Fresnel integral is space invariant and takes the form of a convolution, one can compute it by taking the Fourier transform (FT) of the input, multiplying by the transfer function of free space, and inverse transforming, with overall time ϳN log N. Here we present an algorithm with which to compute Eq. (1) in ϳN log N time, despite the fact that the relationship represented by the more general Eq. (1) is not space invariant and is not a convolution. It is important to underline that here N is chosen close to the space-bandwidth product of the set of input signals, which is usually the smallest possible value of N that can be chosen in terms of information-theoretic considerations. Therefore the algorithm is highly efficient and, although we do not present a formal proof of this, is quite possibly nearly the most efficient that is possible. Indeed, the distinguishing feature of the present approach is the care with which sampling and space-bandwidth product issues are handled. Straightforward use of conventional numerical methods can result in inefficiencies either because their complexity is larger than ϳN log N or because the highly oscillatory quadraticphase kernel in Eq. (1) forces N to be chosen much larger than the space-bandwidth product of the signals.In the method presented here, the samples of the output are computed in terms of the samples of the input, in the sa...

Fundamental structure of Fresnel diffraction: natural sampling grid and the fractional Fourier transform

“…It is known that the Fresnel integral can be decomposed into a FRT, followed by magnification, followed by chirp multiplication [1,[3][4][5][6]:…”

Optimal representation and processing of optical signals in quadratic-phase systems

“…One of the possible decompositions involves three stages. The first is a FRT operation, the second is a magnification operation, and the final stage is a chirp multiplication operation [28][29][30][31][32]:…”

Signal processing issues in diffraction and holographic 3DTV

“…The output of a linear canonical transform can be expressed as [10,27] gðxÞ ¼ e Àjap=4 e Àjpqx 2 f sc ðxÞ, ð23Þ…”

“…Since fractional Fourier transformation corresponds to pure rotation in the space-frequency plane (phase space), neither the spatial nor frequency extent of the signal is altered, and the algorithm of [25] is able to compute the FRT in $N log N time with N being comparable to the space-bandwidth product of f ðxÞ, regardless of the high oscillations of the kernels in question. Furthermore, it is shown in [27] that since the only approximation involved in these computations is that arising from approximate computation of a continuous Fourier transform using a DFT, the accuracy obtained in computing arbitrary linear canonical transforms is likewise comparable to that obtained in using the FFT to compute continuous Fourier transforms. In other words, this algorithm computes linear canonical transforms representing general optical systems, with a performance similar to that of the fast Fourier transform algorithm in computing the Fourier transform, both in terms of speed and accuracy.…”