Fast and accurate simulation of the turbulent phase screen using fast Fourier transform

Abstract. This work describes an accurate method for simulating turbulent phase screens. The phase screen is divided into a fast Fourier transform (FFT)-based screen and a tilt screen. The simulation of the FFT-based screen is different from that of the standard method. In the simulation, the discrete power spectrum of the turbulence is obtained from the discrete Fourier transform of the phase autocorrelation matrix, not from the theoretical power spectrum. This method avoids the drawbacks of the undersampling of the low frequency and high frequency components which occurs in the standard FFT-based method. The maximum error in the phase structure function can be reduced to <0.13%, and the additional execution time increases by only several percents. This method is only suitable for square screens.


Introduction
Simulation of the turbulent phase screens is very important for studying the light propagation and imaging through atmospheric turbulence. There are several methods for simulating turbulent phase screens, e.g., the discrete Fourier transform (DFT) method, 1 the Zernike polynomials method, 2 and the covariance-based method. 3 The Zernike polynomial method is very fast, but it has the drawbacks of inaccurate simulation of the high-spatial frequency turbulence components and difficulty of simulating the effects of the outer scale of turbulence. In the covariance-based method, the simulation is very accurate, but the supported points of the screen are very small. Some interpolation methods can be used to improve the resolution [4][5][6] or to extend the screen size, 7 but the interpolation decreases the accuracy and speed. The DFT method can take the advantage of the fast Fourier transform (FFT) computational performance. This FFT-based method is simple, fast, and less constrained by the computer memory size, it has the ability to generate an extra-large phase screen, and has been widely used for the simulation of the laser beam propagation through atmospheric turbulence. The main drawbacks of the FFT-based screen are the undersampling of the low spatial frequency components and the periodicity induced by the FFT algorithm.
For a square screen with D × D sizes and M × M points, according to Refs. 8 and 9, the FFT-based phase screen can be written as where m, n, m 0 , n 0 ¼ −M∕2; : : : ; M∕2 − 1 are the integer indices, R a ðm 0 ; n 0 Þ and R b ðm 0 ; n 0 Þ are the independent Gaussian random numbers with zero-mean and unit-variance, i ¼ ffiffiffiffiffiffi −1 p denotes the imaginary unit, Φ theory ðm 0 Δ 0 ; n 0 Δ 0 Þ is the discrete power spectrum of turbulence, Δ 0 ¼ 2π∕D and Δ ¼ D∕M are the sampling intervals in the spatial frequencies and spatial domains, respectively, L 0 is the turbulence outer scale, and r 0 is the Fried parameter. Equation (1) can be implemented by using an FFT. One FFT operation can generate two independent phase screens, the real part of ϕ is one and the imaginary part is another.
The FFT-based screen is deficient in the low-spatial frequency and high-spatial frequency turbulence components. The low frequency components have important influence on the low-order turbulent effects. For higher accuracy simulation, the low frequency components should be compensated. Over the past two decades, subharmonic methods 4,8,10 or weighted subharmonic methods 11,12 have been the main solutions to compensate the low frequency components. The subharmonic compensation screens have smaller frequency sampling intervals in the low frequency region, so the discretization error of the power spectrum is reduced. In the subharmonic method of Johansson and Gavel, 8 the low frequency error can be reduced to about 5%. In the weighted subharmonic method of Sedmak, 11 the discrete subharmonic power spectrum at low frequency has been weighted to compensate the low frequency error, and the power spectrum Φ theory around the folding frequency π∕Δ has also been weighted to compensate the high frequency error. After those compensations, the error from the low frequency to high frequency can be reduced to about 1%, but the additional execution time increases by more than two *Address all correspondence to: Jingsong Xiang, E-mail: xiangjs1923@163 .com times compared to the standard FFT-based method and the compensation processes are complicated. Charnotskii 13 proposed a sparse spectrum-based simulation method to reduce the low frequency error. The sampling intervals in the frequency domain are not constant, its exponential increases as frequency increases. In the low frequency region, the sampling intervals are very small and the discretization error is reduced. But this method cannot use the FFT algorithm.
All of those methods are processed in the power spectrum domain. Because the power spectrum is discretely sampled, no matter how the power spectrum is processed the discretization error always exists. The ladder-like discrete power spectrum does not agree well with the theoretical continuous power spectrum and the simulation error always exists.
Xiang 14 proposed a low frequency compensation method. The compensation screen is generated by the covariancebased method. This method avoids the discretization error of the power spectrum in the low frequency region, and the low frequency error can be reduced to <0.1%, but the high frequency error still exists and the execution time increases about 50%.
In this paper, we proposed a very simple and fast method to generate the square phase screen with negligible low frequency and high frequency errors. The screen is divided into an FFT-based screen and a tilt screen. The generating of the FFT-based screen is unlike the conventional method, the discrete power spectrum is obtained from the DFT of the phase autocorrelation matrix, not from the theoretical power spectrum given by Eq. (2). By this means, the discretization error of the power spectrum is avoided.

FFT-Based Phase Screen Simulation According
to the Phase Autocorrelation Function Generally, we use the phase structure function to evaluate the simulation accuracy. The theoretical phase structure function is given by where B theory ðrÞ is the theoretical phase autocorrelation and r is the separation distance between two points. The theoretical phase autocorrelation function can be expressed as 7 where K 5∕6 ð·Þ is the modified Bessel function of the third kind and Γð·Þ is the gamma function. When r ¼ 0, the theoretical phase autocorrelation function is given by For the FFT-based phase screen, the expected phase autocorrelation matrix is an inverse discrete Fourier transform (IDFT) of the power spectrum matrix 8 Conversely, we also can get a power spectrum matrix by using a DFT to the theoretical phase autocorrelation matrix If the theoretical power spectrum Φ theory in Eqs. (1) and (6) is replaced by Φ FFT given by Eq. (7), according to the inverse relationship of Eqs. (6) and (7), and assuming no other errors occurred (an invalid assumption), then the expected phase autocorrelation B exp will be consistent with the theoretical phase autocorrelation B theory , and the simulation error of phase structure function will be eliminated. This is the starting point of our new method. Figure 1 shows the phase structure functions of the phase screen simulated from Eqs. (1) and (7). For the case of L 0 ∕D < ∼1∕2, the simulation error is very small. For the case of large L 0 ∕D, the simulation error is very large. This condition is similar to that of the traditional FFT-based method. If there is no additional processing, this simulation method has no especial advantage over the traditional FFTbased method.
Let us see the reason of the error. The power spectrum Φ FFT given by Eq. (7) is not completely consistent with the theoretical spectrum Φ theory given by Eq. (2), as shown in Fig. 2. Power spectrum Φ FFT oscillates around the theoretical spectrum Φ theory , and some elements in Φ FFT may be negative, whereas the power spectrum in Eq. (1) is required to be non-negative. The negative elements in Φ FFT induce simulation error.
Theoretically, the power spectrum should be real and nonnegative. But here, the power spectrum Φ FFT is obtained from the DFT of the autocorrelation matrix. The autocorrelation matrix cannot contain all autocorrelation information (as shown in Fig. 3), so negative values may appear in Φ FFT , resulting in simulation errors.

Preprocessing of the Phase Autocorrelation
Function In order to eliminate the negative values in Φ FFT , the phase autocorrelation function should be preprocessed. First, we can extract some piston and tilt components from the phase autocorrelation function, and these piston and tilt components are chosen to ensure that the rest of the phase autocorrelation function and its derivative are both zero at r ¼ D∕2. Second, set the phase autocorrelation function to zero when r > D∕2. After such treatments, the phase autocorrelation function with a large L 0 ∕D is very similar to that with a small L 0 ∕D, and the number of negative values in Φ FFT can be greatly reduced. In the following, we describe these steps in detail.
The extracted tilt screen can be expressed as where θ x and θ y are the independent Gaussian random numbers with zero-mean and unit-variance, and σ tilt is the standard deviation of the random tilt angle in the xor ydirections. The phase structure function and the autocorrelation function of the tilt phase screen can be expressed as where B tilt ð0Þ can take on any value, which only influences the piston component, but does not influence the shape of the phase screen. The extracted piston component is chosen to ensure that the rest of the phase autocorrelation is zero at r ¼ D∕2. After extracting out the tilt and piston components, the remaining phase autocorrelation is given by For the case of infinite turbulence outer scale (L 0 → ∞), the phase autocorrelation function B theory cannot be obtained by Eq. (4). In this case, B theory can be obtained by B theroy ðrÞ ¼ B theroy ð0Þ − 6.88ðr∕r 0 Þ 5∕3 ∕2; (12) where B theroy ð0Þ can take on any value, which does not influence the shape of the phase screen.
The tilt screen ϕ tilt should be chosen to ensure that the derivative of B FFT ðrÞ is zero at r ¼ D∕2. This condition can be expressed as where ε is a very short distance <Δ∕100. Inserting Eq. (11) into Eq. (13) leads to By inserting Eq. (10) into Eq. (14), then σ 2 tilt can be obtained as When r > D∕2, set the phase autocorrelation B FFT ðrÞ to zero B FFT ðrÞ ¼ 0; for r > D∕2: (16) From Eq. (16), it is obvious that the simulated phase screen is only valid in the range of r ≤ D∕2. For a rectangular screen with D x × D y sizes, the valid range is r ≤ minðD x ∕2; D y ∕2Þ, so this simulation method is meaningless for the rectangular screens.
We can obtain the discrete power spectrum Φ FFT using the above phase autocorrelation B FFT to replace B theory in Eq. (7). Because some elements in Φ FFT may still be negative, we throw away the negative values to ensure that the discrete power spectrum is non-negative Finally, using the above discrete power spectrum Φ FFT to replace Φ theory in Eq. (1), then the FFT-based phase screen ϕ FFT is obtained.
The final desired phase screen is the sum of the FFT-based screen and the tilt screen, which can be expressed as ϕðmΔ; nΔÞ ¼ ϕ FFT ðmΔ; nΔÞ þ ϕ tilt ðmΔ; nΔÞ: (18)

Evaluation of the Simulation Error
The expected phase structure function of the simulated phase screen is given by where D exp FFT is the expected phase structure function of the FFT-based phase screen ϕ FFT , which can be expressed as where B exp FFT is given by 8 The expected relative simulation error can be defined as re errDðmΔ; nΔÞ¼ D exp ϕ ðmΔ; nΔÞ∕D theory h ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi In the range of r ≤ D∕2, the simulation error is mainly induced by the negative values in Φ FFT obtained by Eq. (7), and essentially no other errors appear in the simulation.
The expected simulation errors are shown in Fig. 4. The high frequency error is a few percents and the low frequency error is near zero. The high frequency errors increase with an increasing M or L 0 ∕D. When L 0 ∕D > 50, the errors are reaching saturation and no longer increase noticeably. The saturation errors are about 4.5%, 2.4%, and 1.1% for M ¼ 2048, 1024, and 512, respectively. The relative error is independent of r 0 .

Compensation of the High Frequency Error
Although this level of simulation accuracy is adequate for most applications, here we will introduce a predistortion method to compensate for the high frequency error. Before using the IDFT to obtain the discrete power spectrum, we distorted the phase autocorrelation value in the opposite direction of the error.
The error of the phase autocorrelation is given by After predistortion, the phase autocorrelation can be expressed as where CðmΔ; nΔÞ is the predistortion coefficient. The discrete power spectrum Φ FFT after high frequency error compensation can be obtained by using B highcom to replace B theory in Eq. (7). Two points should be noted here. The first is that a constant coefficient is not a good choice. The second is that the phase autocorrelation function should be as smooth as possible after predistortion. Figure 5 shows the errors distribution of errB, the error is dependent on azimuth and the points with large error are located near the center of matrix errB, so only the elements in the central region of errB require predistortion.
A Gaussian shaped predistortion coefficients matrix is suggested as CðmΔ; nΔÞ ¼ A expf−½ðmΔÞ 2 þ ðnΔÞ 2 ∕W 2 g: For different A and W, the performance of the high frequency error compensation has some differences, as shown  in Fig. 6. The residual high frequency errors can be reduced by orders of magnitude after high frequency error compensation. Figure 7 shows the maximum relative errors after high frequency error compensation with parameters of A ¼ 1.5 and W ¼ D∕4. The ordinate maxðjre errDjÞ represents the maximum absolute value in the relative error matrix re errD. For M ≤ 2048, the maximum residual error is <0.13%. Figure 8 shows a sample of the simulated phase screen with the parameters of D ¼ 2 m, M ¼ 256, L 0 ¼ 20 m, r 0 ¼ 0.2 m, A ¼ 1.5, and W ¼ D∕4. Figure 9 shows the simulated phase structure function (a) and the error (b) averaged over 2 × 10 6 simulated phase screens. When r ≤ D∕2 the error is very small, only about 0.1%, but when r > D∕2, the error increases rapidly.

Simulation Results
The speed of this method is very fast. On a Lenovo computer running MATLAB R2010a with 2.0 GHz Pentium dual core processor and 2 GB memory, the preparation time and the generating time for one phase screen with M ¼ 512 are about 0.74 and 0.077 s, respectively, whereas those of the standard FFT-based screen are about 0.09 and 0.073 s, respectively. The preparation time of this method is about eight times that of the standard FFT-based method, which is mainly due to the low computational efficiency of the phase autocorrelation matrix using Eq. (4), but this preparation work only needs to be done once. Without accounting for the preparation time, the execution time of generating one screen increases about 6% compared to the standard FFT-based screen.

Conclusion
Compared to other turbulent phase screen simulation methods, this new simulation method has obvious advantages in accuracy and speed. For M < 2048, the overall simulation error from the low spatial frequency to the high spatial frequency can be easily reduced to <0.13%, which is much less than that of any other FFT-based phase screens. The additional execution time only increases about several percents    compared to that of the standard FFT-based phase screen, which is much less than that of any other compensation methods based on the FFT.
The main drawback is that this method is only suitable for square screens, more specifically, the circular screens, and is not suitable for rectangular screens. Another drawback is that the effects of the turbulence inner scale l 0 cannot be easily simulated. The MATLAB code is be available via e-mail.