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...
