Matsushima et al.

Vol. 20, No. 9 / September 2003 / J. Opt. Soc. Am. A

1755

Fast calculation method for optical diffraction
on tilted planes by use of theangular spectrum of plane waves
Kyoji Matsushima
Department of Electrical Engineering, Kansai University, 3-3-35 Yamate-cho, Suita, Osaka 564-8680, Japan

Hagen SchimmelLightTrans GmbH, Jena D-07745, Germany

Frank Wyrowski
Friedrich Schiller University Jena, Jena D-07743, Germany
Received March 22, 2003; revised manuscript received May 8, 2003; accepted May 8, 2003
A novel method for simulating field propagation is presented. The method, based on the angular spectrum ofplane waves and coordinate rotation in the Fourier domain, removes geometric limitations posed by conventional propagation calculation and enables us to calculate complex amplitudes of diffracted waves on a plane
not parallel to the aperture. This method can be implemented by using the fast Fourier transformation twice
and a spectrum interpolation. It features computation time that is comparable with that of standard calculation methods for diffraction or propagation between parallel planes. To demonstrate the method, numericalresults as well as a general formulation are reported for a single-axis rotation. © 2003 Optical Society ofAmerica
OCIS codes: 050.1960, 350.5500, 090.1970, 000.3860.

1. INTRODUCTION
Developments in computer technology in the past few decades have created a revolution in the field of wave optics.
Today our desktop computers can easily simulate optical
diffraction by an aperture. The complex amplitudes insimulations are typically sampled on an equidistant grid.
These simulations of diffraction are usually based onwell-known scalar theories such as the Rayleigh–Sommerfeld formula, the Fresnel–Kirchhoff formula, andthe Fresnel and Fraunhofer approximations of the latter.
Another important method for computing optical diffraction is to use the angular spectrum of plane waves.1 Thismethod, accompanied by the fast-Fourier-transformation(FFT) algorithm, plays an important role in numericalsimulations of optical diffraction and propagation. However, these methods share a common constraint in simulating diffraction: The screen or plane on which lightwaves diffracted by an aperture are calculated must beparallel to the aperture. No tilted screens or tilted apertures are allowed in standard formulations. To get rid of
this limitation, diffraction by a tilted aperture has been
investigated in some studies: Fraunhofer diffraction in
Refs. 2–5 and Fresnel diffraction in Ref. 6. Furthermore,in the context of creating image holograms, Leseberg has
suggested using rotational transformation of complex amplitudes in the Fourier domain in order to calculate object
waves diffracted by tilted apertures.7 However, we havefound no complete formulations or results on such rotation in the literature.
1084-7529/2003/091755-08$15.00

New methods of simulation for optical diffraction andpropagation with fewer constraints and higher performance are needed to improve optical techniques such as
diffractive optics and digital holography, which allow us
to measure8 and create9 complex amplitudes of lightwaves emitted by three-dimensional objects. In this paper we treat the formulation and numerical simulation ofdiffraction by an aperture onto an inclined plane, i.e., wecalculate complex amplitudes at sampling points placedon a plane not parallel to the aperture. Although our
major interest lies in three-dimensional displays,10 thetechnique presented is useful in any wave-optical simulation concerning tilted planes.

2. BASIC THEORY AND FORMULATION
A. Basic Idea and Concept
Two coordinate systems are used in this investigation.
One is the (x, y, z) source coordinate system, in which
the complex amplitude g(x, y) of a light wave is initiallydefined in the (x, y, 0) plane. This plane is referred to as
the source plane in this paper. The other system is thereference coordinate system (x̂, ŷ, ẑ) at a distance d fromthe source coordinates, in which the complex amplitudef(x̂, ŷ) of the propagated field is finally obtained in the(x̂, ŷ, 0) plane, which is referred to as the reference
plane.
Figure 1 illustrates the specified coordinate systems
with the source and reference planes and geometries cor© 2003 Optical Society of America

1756

J. Opt. Soc. Am. A / Vol. 20, No. 9 / September 2003

Matsushima et al.

G 共 u, v 兲 ⫽ F 兵 g 共 x, y 兲 其
⫽

冕冕

⬁

g 共 x, y 兲 exp关 ⫺i2 ␲ 共 ux ⫹ vy 兲兴 dxdy

(1)

⫺⬁

and the inverse transformation by
g 共 x, y 兲 ⫽ F ⫺1 兵 G 共 u, v 兲 其⫽Fig. 1. Definition of coordinate system and schemes for (a) atilted screen and (b) a tilted aperture.

冕冕⬁G 共 u, v 兲 exp关 i2 ␲ 共 ux ⫹ vy 兲兴 dudv.

(2)

⫺⬁

Equation (2) can be interpreted as a superimposition ofplane waves:u 共 x, y, z ⫽ 0; u, v 兲 ⫽ G 共 u, v 兲 exp关 i2 ␲ 共 ux ⫹ vy 兲兴 .
(3)
G(u, v) expresses the complex amplitude of the planewave. The same plane wave with a complex amplitude A
and a wave vector k can be defined in general as冋

u 共 x, y, z ⫽ 0; k兲 ⫽ A exp i

2␲
␭

册

共 k xx ⫹ k yy 兲 ,

(4)

with
k⫽
Fig. 2. Geometry for (a) a tilted screen and (b) a tilted aperture.

responding to (a) a tilted screen and (b) a tilted aperture.
In (a), the source plane is perpendicular to the z axis,
while the reference plane is not. The inverse situation isdepicted in (b). Figure 2 shows schemes (a) and (b) from
Fig. 1 depicted for a common aperture placed in thesource plane. Figure 2(a) illustrates the situation for atilted screen and Fig. 2(b) for a tilted aperture.
In both schemes the calculation of diffraction consistsof two essential steps. For the inclined reference plane,light diffracted by the aperture is propagated to an intermediate plane (x ⬘ , ŷ) that is parallel to the source plane.
This calculation can be carried out with conventional
methods based on the traditional diffraction theory.
Next the calculated field given as complex amplitudes on(x ⬘ , ŷ) is decomposed into plane waves by Fourier transformation. The expanded plane waves are reassembled
in the tilted reference plane (x̂, ŷ) and then combined
again by using an inverse Fourier transformation. This
step is referred to as coordinate rotation in the Fourier
domain. The same procedure is generally applicable to atilted source plane, as shown in Fig. 2(b); i.e., source complex amplitudes are propagated with conventional methods, followed by coordinate rotation.
In either scheme, the most essential and significant operation is the coordinate rotation in the Fourier domain,which will be described in the following subsection.

2␲

关kx

␭

kz兴.

(5)

Therefore 兩 k兩 ⫽ 2 ␲ /␭ and k x 2 ⫹ k y 2 ⫹ k z 2 ⫽ 1. By comparing Eq. (3) with Eq. (4), wave vector k can be associated with Fourier frequencies (u, v) in the source planebyk ⫽ 2␲关u⫺2v

w 共 u, v 兲兴 ,(6)

with w(u, v) ⫽ (␭ ⫺ u ⫺ v ) .
Let us now turn to the field in the reference plane,
which is denoted f(x̂, ŷ). Its spectrum is given by
2

2 1/2

F 共 û, v̂ 兲 ⫽ F 兵 f 共 û, v̂ 兲 其⫽

冕

⬁f 共 x̂, ŷ 兲 exp关 ⫺i2 ␲ 共 ûx̂ ⫹ v̂ŷ 兲兴 dx̂dŷ.

⫺⬁

(7)

Similarly to obtaining Eq. (4), we obtaink̂ ⫽ 2 ␲ 关 ûv̂ŵ 共 û, v̂ 兲兴 ,(8)with ŵ(û, v̂) ⫽ (␭ ⫺2 ⫺ û 2 ⫺ v̂ 2 ) 1/2.
The wave vectors given in Eqs. (6) and (8) are defined
respectively, in the (x, y, z) and (x̂, ŷ, ẑ) coordinate systems. Therefore if the two fields represent the samewave, the wave vector expressed in one coordinate system
can be transformed into the other by using ordinary
coordinate-transformation procedures.
Suppose that T is a matrix used to transform thesource coordinates into the reference coordinates; then
the wave vectors can be transformed into each other byk ⫽ T⫺1 k̂.

k̂ ⫽ Tk,
B. Coordinate Rotation in the Fourier Domain
In this subsection, assume that d ⫽ 0, i.e., the two coordinates share the origin, so that we can concentrate on
the formulation of the coordinate rotation. If the twodimensional complex amplitude g(x, y, 0) is given in thesource plane, the spectrum G(u, v) is defined by

ky

(9)

When the inverse transformation matrix is given as
⫺1

T

⫽

冋

册

a1

a2

a3

a4

a5

a6 ,

a7

a8

a9

(10)

Matsushima et al.

Vol. 20, No. 9 / September 2003 / J. Opt. Soc. Am. A

we can associate the Fourier frequencies in the source coordinates with those in the reference coordinates as follows:u ⫽ ␣ 共 û, v̂ 兲 ⫽ a 1 û ⫹ a 2 v̂ ⫹ a 3 ŵ 共 û, v̂ 兲 ,v ⫽ ␤ 共 û, v̂ 兲 ⫽ a 4 û ⫹ a 5 v̂ ⫹ a 6 ŵ 共 û, v̂ 兲 .

(11)

f 共 x̂, ŷ 兲 ⫽ F ⫺1 兵 F 共 û, v̂ 兲 兩 J 共 û, v̂ 兲 兩 其 .

(12)

It appears as if it is possible to calculate the complex amplitude in the reference plane by inverse transformationof the spectrum F(û, v̂). However, a simple inverse Fourier transformation of F(û, v̂) does not provide correct results, because the total energy of the field is not conserved
after rotational transformation, owing to the nonlinearityof the transformation in Eq. (11).
The energy of a plane wave in Eq. (3) in source coordinates, which have Fourier frequencies (u, v) to (u⫹ du, v ⫹ dv), is proportional to 兩 G(u, v)exp关i2␲(ux⫹ vy)兴兩2dudv, and therefore the total energy of the field isgiven by
E total ⬀冕冕⬁兩 G 共 u, v 兲 兩 dudv.
2

(13)

⫺⬁

In the transformed spectrum F(û, v̂) ⫽ G(␣ (û, v̂),␤ (û, v̂)), the energy factor of a plane wave with Fourier
frequencies (û, v̂) to (û ⫹ dû, v̂ ⫹ dv̂) is written as兩 G(␣ 共 û, v̂ 兲 , ␤ 共 û, v̂ 兲 )兩 2 兩 J 共 û, v̂ 兲 兩 dûdv̂,

where the Jacobian J(û, v̂) is defined as
J 共 û, v̂ 兲 ⫽⳵␣ ⳵␤⳵ û ⳵ v̂

⳵␣ ⳵␤

⫺

⳵ v̂ ⳵ û⫽ 共 a 2a 6 ⫺ a 3a 5 兲ûŵ 共 û, v̂ 兲

⫹ 共 a 3a 4 ⫺ a 1a 6 兲v̂
ŵ 共 û, v̂ 兲

⫹ 共 a 1a 5 ⫺ a 2a 4 兲,
(14)which is needed to compensate for the nonlinearity of therotational transformation. Thus the total energy is written as
E total ⬀冕冕⬁兩 G(␣ 共 û, v̂ 兲 , ␤ 共 û, v̂ 兲 )兩 2 兩 J 共 û, v̂ 兲 兩 dûdv̂.

⫺⬁J 共 û, v̂ 兲 ⯝ 共 a 1 a 5 ⫺ a 2 a 4 兲

f 共 x̂, ŷ 兲 ⫽冕冕

⬁F 共 û, v̂ 兲 exp关 i2 ␲ 共 ûx̂ ⫹ v̂ŷ 兲兴⫺⬁⫻ 兩 J 共 û, v̂ 兲 兩 dûdv̂.

(16)

This is rewritten with the inverse Fourier transformationand yields(18)

is a good approximation.
On the other hand, when the waves are paraxial in thesource coordinates, the Jacobian remaining after the approximation is a little complex, as described in Appendix
A, but all terms that include the Fourier frequencies canbe ignored. As a result, the integration of the rotatedspectrum for the paraxial waves is reduced to the inverse
Fourier transformationf 共 x̂, ŷ 兲 ⬀ F ⫺1 兵 F 共 û, v̂ 兲 其 .

(19)

We refer to this as the paraxial approximation in the coordinate rotation.

3. SINGLE-AXIS ROTATION OF THE
REFERENCE COORDINATESIn this section we treat the single-axis rotation of the reference coordinates, of which the scheme is depicted inFig. 1(a), to verify the basic idea of rotating the coordinates in the Fourier domain.
A. Analytical Procedure
Suppose that a source field is given by g(x, y) and its
spectrum is given again by G(u, v) in the source plane,
but in this subsection the reference plane is assumed tobe apart from the source plane; i.e., d ⬎ 0. In addition,
suppose that G d (u, v) denotes the source spectrum in anintermediate plane, which is parallel to the source plane,
and that its coordinates share the origin with the reference coordinates. One can choose an analytical or a numerical method from some conventional algorithms offield propagation to calculate fields in the intermediate
plane, because this propagation is translational and is independent of the step in the coordinate rotation. In thispaper we choose a standard formula of the angular spectrum of plane waves given as follows:
G d 共 u, v 兲 ⫽ G 共 u, v 兲 exp关 i2 ␲ d 共 ␭ ⫺2 ⫺ u 2 ⫺ v 2 兲 1/2兴 .
(20)
The transformation matrix used to rotate coordinates
on the y axis with the angle of ␸ is given as(15)
To obtain complex amplitudes of the field in the reference coordinates that posses the same energy as thatin the source coordinates, the plane waves
F(û, v̂)exp关i2␲(ûx̂ ⫹ v̂ŷ)兴 must be integrated as follows:(17)

When waves are paraxial in the reference coordinates[i.e., fields propagate in the ẑ direction as shown in Fig.
1(b)], û and v̂ are much smaller than ŵ(û, v̂). As a result, the first and second terms of the Jacobian in Eq. (14)can be ignored. In that case,Thus the spectrum in the reference coordinates is givenby the spectrum in the source coordinates according to
F 共 û, v̂ 兲 ⫽ G(␣ 共 û, v̂ 兲 , ␤ 共 û, v̂ 兲 ).

1757

⫺1

T

⫽

冋

cos ␸

0

sin ␸

0

1

0

⫺sin ␸

0

cos ␸

册

.

(21)

By substituting elements of matrix (21) into Eqs. (11) and(14), we can write the angular spectrum in the reference
coordinates and the Jacobian, respectively, as follows:
F 共 û, v̂ 兲 ⫽ G d 共 û cos ␸ ⫹ ŵ 共 û, v̂ 兲 sin ␸ , v̂ 兲 ,J 共 û, v̂ 兲 ⫽ cos ␸ ⫺ûŵ 共 û, v̂ 兲

sin ␸ .

(22)
(23)

1758

J. Opt. Soc. Am. A / Vol. 20, No. 9 / September 2003

Matsushima et al.

nates. Assume that the center frequency is u 0 (⫽0) in
the source spectrum, as shown in Fig. 4(a), and thatû 0 (⫽0) is the carrier frequency after rotation, as shownin Fig. 4(b). For an efficient conversion by FFT, anothercoordinate u ⬘ ⫽ û ⫺ û 0 should be introduced, as shownin Fig. 4(c). This results in complex amplitudes in thereference coordinates, which are given as
f 共 x̂, ŷ 兲 ⫽ F ⫺1 兵 F 共 u ⬘ , v̂ 兲 兩 J 共 u ⬘ , v̂ 兲 兩 其 exp关 i2 ␲ û 0 x 兴 .

Fig. 3. Curves for transforming Fourier frequency u in source
coordinates into û in reference coordinates; v ⫽ 0.

Therefore complex amplitudes in the reference plane are
given by the inverse Fourier transformation

再f 共 x̂, ŷ 兲 ⫽ F ⫺1 G d 共 û cos ␸ ⫹ ŵ 共 û, v̂ 兲 sin ␸ , v̂ 兲冏⫻ cos ␸ ⫺ûŵ 共 û, v̂ 兲

sin ␸

冏冎

.

(24)

In the case of paraxial waves or small rotation angles,
the second term on the right-hand side of Eq. (23) can beignored. As a result, complex amplitudes in the paraxial
approximation (19) are written as(26)

The carrier frequency û 0 is possibly too high to besampled in a given sampling scheme in the reference
plane. Even in this case, the intensity distribution
兩 f(x̂, ŷ) 兩 2 can be easily obtained, because the carrier frequency is included only as a phase factor of exp (i2␲û 0x).
Moreover, when one wants to consider the factor of thecarrier frequency, it can be multiplied after inverse FFT.
C. Numerical Examples for Single-Axis Rotation
We prepared three numerical examples for the source amplitude g(x, y) to verify the theory and evaluate numerical accuracy. We calculated complex amplitudes in a reference plane at z ⫽ d, which was rotated at an angle ␸ on
the y axis, as shown in Fig. 5.
Parameters of the three samples are summarized inTable 1. Sample A is a paraxial Gaussian wave with aconstant phase (i.e., the wave’s waist) in the source plane.
This is an example of a simple paraxial wave with smooth
amplitude distribution. Sample B is a plane wave diffracted by a rectangular aperture in the source plane.

f 共 x̂, ŷ 兲 ⬀ F ⫺1 兵 G d 共 û cos ␸ ⫹ ŵ 共 û, v̂ 兲 sin ␸ , v̂ 兲 其 . (25)
B. Numerical Procedure
Equation (24) is rigorous in an analytic sense. However,
the transformation in Eq. (22) causes numerical problems. One problem is the distortion of the equidistantsampling grid. For a fast calculation of diffraction forpractical purposes, the FFT algorithm is commonly used
to calculate G(u, v) from g(x, y) and f(x̂, ŷ) fromF(û, v̂). The FFT algorithm, however, works only for anequidistant sampling grid. Since the inverse functionû ⫽ ␣ ⫺1 (u, 0), as well as the function u ⫽ ␣ (û, 0) itself,
is a nonlinear function as shown in Fig. 3, even if G(u, v)
and G d (u, v) are evenly sampled, the spacing of the sampling points for the transformed spectrum is not equal.
As a result, some method of interpolation must be introduced into the sampled F(û, v̂) to calculate f(x̂, ŷ) by
FFT.
Furthermore, the transformation from Eq. (22) also
causes a shift in the center frequency in the spectrum
F(û, v̂). When u ⫽ 0 and v ⫽ 0 in the source coordinates, the corresponding reference frequency û⫽ ␣ ⫺1 (0, 0) does not equal zero. This can be interpreted
to be a carrier frequency observed in the tilted reference
coordinates. If one wants to calculate the diffracted
waves in the reference plane by FFT, this nonzero carrier
frequency should be eliminated, because the FFT algorithm does not work effectively for spectrums sampled far
from zero frequency. Figure 4 shows the procedure used
to eliminate the carrier frequency in the reference coordi-

Fig. 4. Schematic diagram of spectrum shift: (a) spectrum insource coordinates, (b) reference coordinates, and (c) shift of spectral origin.

Fig. 5. Schematic setup for numerical simulation of single-axis
rotation. One-dimensional rigorous amplitudes for comparison
are calculated at sampling points (dots).

Matsushima et al.

Vol. 20, No. 9 / September 2003 / J. Opt. Soc. Am. A

Table 1. Summary of Parameters Used for
Numerical Simulations
Sampling
g(x, y)NumberPitchdSample A, Gaussian
Sample B, rectangle aperture
Sample C, Gaussian with quadratic
phase

256
256
512

20 ␮m
20 ␮m
5 ␮m

10 mm
30 mm
10 mm

The amplitude distribution of this wave is not smooth in
the reference plane because of the fringe pattern.
Sample C is also Gaussian, but it has a spherical phase in
the source plane to simulate a Gaussian beam behind alens with a focal length of 15 mm. This phase brings
asymmetry as well as nonparaxiality into the distributionof the complex amplitudes in the tilted reference plane
placed near the focal point.
D. Numerical Accuracy
As mentioned at the beginning of Subsection 3.B, approximations are not applied in the derivation of Eq. (24); i.e.,
the equation itself is rigorous. However, the interpolation of the sampled spectrum, which is necessary for thenumerical calculation, causes numerical errors. Therefore the accuracy as well as the computation time of thenumerical calculation is dependent on the method of interpolation. According to the sampling theorem, the sincinterpolation should provide the most accurate results,
but the sinc interpolation is generally a very timeconsuming method. Therefore we examine other interpolation methods in this section.
The numerical accuracy of the coordinate rotation was
evaluated as signal-to-noise ratio (SNR) of rotated complex amplitudes defined as

SNR ⫽

冕冕

⬁

冕冕

⬁

dicular to the z axis and placed before and behind the position d; then the resultant amplitudes at sampling points
were stitched along side the tilted reference plane. Although the resulting amplitudes may also create numerical errors, the accuracy of the method is high enough toprovide a reference for comparison with complex amplitudes calculated by the proposed coordinate rotation. In
addition, calculating exact amplitudes with this method
was possible only for a one-dimensional distribution because of the long computation time.
Figure 6 shows SNRs measured numerically for a rotation angle of 30° in various interpolation algorithms of theFourier spectrum. ‘‘nth’’ and ‘‘Cubic 8’’ indicate thenth-order iterated linear interpolation and the cubicinterpolation algorithms at 8 points,11 respectively. The
cubic interpolation almost always results in the best
SNRs. These SNRs exceeded 50 dB for all samples withan angle of 30°.
In Fig. 7, SNRs measured in the cubic interpolation are
shown in relation to the rotation angles. In samples A
and B, the paraxial approximation of Eq. (25) gives thesame SNRs as Eq. (24), whereas in sample C, the SNR
calculated with the approximation decreases with an increase in rotation angle. The reason for this is clear:samples A and B are paraxial waves, whereas sample C,
which has a definite beam waist, deviates from theparaxial condition. Thus the accuracy of numerical cal-兩 f 共 x̂, ŷ 兲 兩 2 dx̂dŷ

,

(27)

兩 f 共 x̂, ŷ 兲 ⫺ ␣ f rig共 x̂, ŷ 兲 兩 dx̂dŷ
2where f rig(x̂, ŷ) is the distribution of complex amplitudes
calculated by another rigorous method described below.
In most cases, magnitude of the field is of no concern, andaccuracy in the field form is of interest to us. Therefore ascale factor ␣ maximizing SNR is defined as follows:

␣⫽

1759

冕冕
冕冕
⬁

Fig. 6. SNRs in one-dimensional amplitudes for several interpolation algorithms. Here nth and Cubic 8 stand for nth iterative linear interpolation and cubic interpolation at 8 points, respectively. The rotation angle is 30°.

* 共 x̂, ŷ 兲 dx̂dŷ
f 共 x̂, ŷ 兲 f rig

⬁

.

(28)

兩 f rig共 x̂, ŷ 兲 兩 2 dx̂dŷ

A rigorous calculation of complex amplitudes in thetilted reference plane is a serious problem because of its
mathematical complexity. In almost all cases it seems to
be impossible to obtain analytically rigorous solutions of
the field in a tilted plane, and therefore we adopted a numerical approach. Figure 5 shows how we obtained theexact amplitudes in the tilted plane in this paper. Complex amplitudes were calculated on many planes perpen-

Fig. 7. SNR versus rotation angle.
points is used in calculation.

Cubic interpolation at 8

1760

J. Opt. Soc. Am. A / Vol. 20, No. 9 / September 2003

Matsushima et al.

Fig. 8. Two-dimensional distribution of amplitudes and phase calculated on tilted reference planes. The phase factor due to the carrier
frequency appearing on the tilted plane is eliminated in depicted phase images.

tions of the FFT consume most of the processing time andthe rest of the time is spent on interpolation and propagation, with computational costs that are proportional ton.

4. DISCUSSION

Fig. 9. Computation time for rotation and translational propagation. Measurements were made by using a Pentium III processor with 1 GHz frequency.

culation obtained by Eq. (25) drops with increasing degreeof the rotation angle.
E. Two-Dimensional Results and Computation Time
Figure 8 shows the two-dimensional distributions of complex amplitudes for numerical samples A–C. The left
and right columns for each sample depict amplitude andphase distribution, respectively. Note that the carrier
frequency û 0 mentioned in Subsection 3.B has been eliminated in the phase distributions.
The computation time for rotation accompanied bytranslational propagation was measured on an Intel Pentium III processor at 1 GHz. Figure 9 shows the time atdifferent numbers of sampling points. In the measurement, the single-axis rotation was again calculated withcubic interpolation. The results of the computation time
approximately agree with those of FFT, i.e., n log n for nsamplings. This is attributed to the fact that two itera-

As shown in Fig. 3, the curves for transformation of frequency u into û open upward with positive values of ␸and open downward with negative ␸. Therefore thecurves have a minimum value or a maximum value except for the case ␸ ⫽ 0. In some cases of large-angle rotation, this causes a ‘‘fold in the spectrum’’ in the reference coordinates. For example, when ␸ ⫽ 90° in Fig. 3,
two source frequencies ⫹u ⬘ and ⫺u ⬘ are assigned to an
identical frequency û ⬘ of the reference spectrum. Thismeans that the source spectrum is folded at its origin; i.e.,
both negative and positive frequency ranges in the source
coordinates are projected onto the same range of the reference frequency.
The spectral fold originates in a group of plane waves
that propagates from the region ẑ ⬎ 0 to the regionẑ ⬍ 0 in reference coordinates, as shown in the following.
In the single-axis rotation, the frequency transformation
is written asû ⫽ u cos ␸ ⫺ w 共 u, v 兲 sin ␸ ,(29)ŵ ⫽ u sin ␸ ⫹ w 共 u, v 兲 cos ␸ .

(30)

The maximum and the minimum values are given at u
satisfying dû/du ⫽ 0, and thus
u sin ␸ ⫹ w 共 u, v 兲 cos ␸ ⫽ 0,(31)

Matsushima et al.

Vol. 20, No. 9 / September 2003 / J. Opt. Soc. Am. A

1761

measured on a Pentium III processor at 1 GHz. The computation time was similar to that of FFT, for example,
35.6 s in the case of 20482 sampling points.
Finally, a physical interpretation for the spectral fold
was discussed in relation to the propagation direction ofplane waves in the reference coordinate system.

APPENDIX A: PARAXIAL APPROXIMATION
IN THE SOURCE COORDINATES
Fig. 10. Schematic depiction of the spectral fold. Each sidewith respect to the maximum or minimum value in frequency
transformation curves corresponds to the sign of the k̂ z elementof wave vectors in the reference coordinates.

and the z element of the wave vector in the reference coordinates is given by substituting Eq. (30) into Eq. (8) as
follows:k̂ z 共 u, v 兲 ⫽ 2 ␲ 关 u sin ␸ ⫹ w 共 u, v 兲 cos ␸ 兴 .

(32)

If u ⫽ u fold satisfies the condition defined in Eq. (31),u fold leads to k̂ z (u fold , v) ⫽ 0. Therefore the sign of k̂ zchanges at u ⫽ u fold when u increases or decreases. The
plane wave with a frequency 兩 u 兩 ⬍ 兩 u fold兩 in the source
coordinates consequently propagates in the directionẑ ⬎ 0, and vice versa. This relation is depicted in Fig.
10. Accordingly, the spectral fold means that two planewaves having an opposite sign of k̂ z to each other, i.e.,having a mirrorlike symmetry with respect to the reference plane, contribute to the same plane.
This spectral fold does not occur in our numerical examples, because the spectrum is folded at 兩 u fold兩⫽ w(0, v)cos ␸ in the single-axis rotation, and thus
兩 u fold兩 ⫽ cos ␸/␭ in v ⫽ 0. In the examples presented, thebroadest source spectrum is provided in u ⫽ ⫾0.063/␭
from Sample C, whereas 兩 u fold兩 ⫽ 0.174/␭ even in the
maximum rotation angle of 80°.

When the inverse rotation matrix T⫺1 is defined by Eq.
(10), assume that the rotation matrix is given byT⫽兩 T⫺1 兩

冋

册

A1

A2

A3

A4

A5

A6 .

A7

A8

A9

(A1)

Here A i is a cofactor of the matrix T⫺1 , and A i ’s in thethird column are associated with the matrix T⫺1 by

冏

a2

a3

a5

a6

冏

a1

a3

a4

a6

冏

a1

a2

a4

a5

A 3 ⫽ 共 ⫺1 兲 共 1⫹3 兲

A 6 ⫽ 共 ⫺1 兲 共 3⫹2 兲

A 9 ⫽ 共 ⫺1 兲 共 3⫹3 兲冏⫽ a 2a 6 ⫺ a 3a 5 ,

(A2)

冏

⫽ a 3a 4 ⫺ a 1a 6 .

(A3)

冏

⫽ a 1a 5 ⫺ a 2a 4 .

(A4)Thus the Jacobian (14) is rewritten by using the cofactors:
J 共 û, v̂ 兲 ⫽ A 3

û
ŵ

⫹ A6

v̂
ŵ

⫹ A9 .

(A5)

Moreover, the first term of the Jacobian is rewritten bytransforming frequencies (û, v̂, ŵ) into (u, v, w)
through the matrix T:û
A3

5. CONCLUSION
A method for calculating complex amplitudes of a diffracted wave on an inclined screen as well as on a wave
diffracted by an inclined aperture was suggested. Thismethod features a computation time that is comparablewith that of standard methods in which a reference plane
is parallel to a source plane. The procedure of thismethod can be summarized as follows: First calculate
the spectrum of the source amplitude by FFT, then rotate
the coordinate system in which the spectrum is given, andfinally reintegrate the spectrum in rotated coordinates byusing inverse FFT with a Jacobian because of the nonlinearity of the rotational transformation of the spectrum.
In the examples of single-axis rotation, SNRs of complex amplitudes on the tilted reference plane are measured by comparing them with rigorous results. In the
one-dimensional distribution of complex amplitudes,
these SNRs exceeded 50 dB in almost all examples androtation angles. Two-dimensional distributions of complex amplitudes in the rotated reference plane were calculated for these examples, and the computation time was

1

ŵ

⫽ A3

A 1u ⫹ A 2v ⫹ A 3w
A 7u ⫹ A 8v ⫹ A 9w

.

(A6)

If a field is paraxial in the source coordinate, i.e., if itpropagates in almost the z direction, the frequency w ismuch larger than u and v. Therefore
û
A3

ŵ

⫽ A3

A 1 u/w ⫹ A 2 v/w ⫹ A 3
A 7 u/w ⫹ A 8 v/w ⫹ A 9

⯝

A 32
A9

,

(A7)

is a good approximation. With the same procedure, thesecond term of the Jacobian (A5) is approximated to
v̂
A6

ŵ

⯝

A 62
A9

.

(A8)

As a result, when fields are paraxial in the source coordinates, the paraxial approximation of coordinate rotation
is given by
J 共 û, v̂ 兲 ⯝A 32 ⫹ A 62 ⫹ A 92
A9

.

(A9)

1762

J. Opt. Soc. Am. A / Vol. 20, No. 9 / September 2003

Matsushima et al.

ACKNOWLEDGMENT
This work is partly supported by Kansai University HighTechnology Research Center.

7.
8.

REFERENCES
1.
2.
3.
4.
5.
6.

J. W. Goodman, Introduction to Fourier Optics, 2nd ed.
(McGraw-Hill, New York, 1996), Chap. 3.10.
G. Harburn, J. K. Ranniko, and R. P. Williams, ‘‘An aspectof phase in Fraunhofer diffraction patterns,’’ Optik (Stuttgart) 48, 321–328 (1977).
S. Ganci, ‘‘Fourier diffraction through a tilted slit,’’ Eur. J.
Phys. 2, 158–160 (1981).
K. Patorski, ‘‘Fraunhofer diffraction patterns of tilted planar objects,’’ Opt. Acta 30, 673–679 (1983).
H. J. Rabal, N. Bolognini, and E. E. Sicre, ‘‘Diffraction by atilted aperture: coherent and partially coherent cases,’’
Opt. Acta 32, 1309–1311 (1985).
D. Leseberg and C. Frère, ‘‘Computer-generated holograms

9.
10.

11.

of 3-D objects composed of tilted planar segments,’’ Appl.
Opt. 27, 3020–3024 (1988).
D. Leseberg, ‘‘Computer-generated three-dimensional image holograms,’’ Appl. Opt. 31, 223–229 (1992).
Y. Takaki and H. Ohzu, ‘‘Hybrid holographic microscopy:
visualization of three-dimensional object informationby use of viewing angles,’’ Appl. Opt. 39, 5302–5308 (2000).
K. Matsushima and M. Takai, ‘‘Recurrence formulas for fastcreation of synthetic three-dimensional holograms,’’ Appl.
Opt. 39, 6587–6594 (2000).
K. Matsushima, H. Schimmel, and F. Wyrowski, ‘‘New creation algorithm for digitally synthesized holograms in surface model by diffraction from tilted planes,’’ in Practical
Holography XVI and Holographic Materials VIII, S. A. Benton, S. H. Stevenson, and T. J. Trout, eds., Proc. SPIE 4659,
53–60 (2002).
T. M. Lehmann, C. Gönner, and K. Spitzer, ‘‘Large-sized local interpolators,’’ in Proceedings of The International Association of Science and Technology for Development(IASTED) Conference on Computer Graphics and Imaging(ACTA Press, Calgary, Alberta, Canada, 1999), pp. 156–
161.