1.

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,¹ the Zernike polynomials method,² and the covariance-based method.³ 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–6 or to extend the screen size,⁷ 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\times D$ sizes and $M\times M$ points, according to Refs. 8 and 9, the FFT-based phase screen can be written as

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,⁸ the low frequency error can be reduced to about 5%. In the weighted subharmonic method of Sedmak,¹¹ the discrete subharmonic power spectrum at low frequency has been weighted to compensate the low frequency error, and the power spectrum ${\mathrm{\Phi }}_{\text{theory}}$ around the folding frequency $\pi /\mathrm{\Delta }$ 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 times compared to the standard FFT-based method and the compensation processes are complicated. Charnotskii¹³ 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¹⁴ proposed a low frequency compensation method. The compensation screen is generated by the covariance-based 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.

2.

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

The theoretical phase autocorrelation function can be expressed as⁷

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⁸

Conversely, we also can get a power spectrum matrix by using a DFT to the theoretical phase autocorrelation matrix

If the theoretical power spectrum ${\mathrm{\Phi }}_{\text{theory}}$ in Eqs. (1) and (6) is replaced by ${\mathrm{\Phi }}_{\mathrm{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_{\text{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<\sim 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 FFT-based method.

Fig. 1

Phase structure functions of the simulated screen as a function of the separation distance. The curves are averaged over ${10}^5$ phase screens. For the case of $L_0=1\mathrm{m}$, the simulated curve agrees well with the theoretical curve, it is hard to distinguish them.

Let us see the reason of the error. The power spectrum ${\mathrm{\Phi }}_{\mathrm{FFT}}$ given by Eq. (7) is not completely consistent with the theoretical spectrum ${\mathrm{\Phi }}_{\text{theory}}$ given by Eq. (2), as shown in Fig. 2. Power spectrum ${\mathrm{\Phi }}_{\mathrm{FFT}}$ oscillates around the theoretical spectrum ${\mathrm{\Phi }}_{\text{theory}}$, and some elements in ${\mathrm{\Phi }}_{\mathrm{FFT}}$ may be negative, whereas the power spectrum in Eq. (1) is required to be non-negative. The negative elements in ${\mathrm{\Phi }}_{\mathrm{FFT}}$ induce simulation error.

Fig. 2

Comparison of power spectrum density given by Eqs. (7) and (2). Simulation parameters are $D=2\mathrm{m}$, $M=256$, $r_0=0.2\mathrm{m}$, and $L_0=1\mathrm{m}$.

Theoretically, the power spectrum should be real and non-negative. But here, the power spectrum ${\mathrm{\Phi }}_{\mathrm{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 ${\mathrm{\Phi }}_{\mathrm{FFT}}$, resulting in simulation errors.

Fig. 3

The diagrams of the phase autocorrelation. Simulation parameters are $D=2\mathrm{m}$, $M=256$, and $r_0=0.2\mathrm{m}$. The left is for the case of $L_0/D=0.5$, the autocorrelation matrix contains almost all autocorrelation information (when $r>D/2$, $B(r)\approx 0$). The right is for the case of $L_0/D=2.5$, the autocorrelation matrix cannot contains all autocorrelation information, some autocorrelation information for $r>D/2$ is lost.

3.

Preprocessing of the Phase Autocorrelation Function

In order to eliminate the negative values in ${\mathrm{\Phi }}_{\mathrm{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 ${\mathrm{\Phi }}_{\mathrm{FFT}}$ can be greatly reduced. In the following, we describe these steps in detail.

The extracted tilt screen can be expressed as

The phase structure function and the autocorrelation function of the tilt phase screen can be expressed as

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\to \infty$), the phase autocorrelation function $B_{\text{theory}}$ cannot be obtained by Eq. (4). In this case, $B_{\text{theory}}$ can be obtained by

The tilt screen ${\varphi }_{\text{tilt}}$ should be chosen to ensure that the derivative of $B_{\mathrm{FFT}}(r)$ is zero at $r=D/2$. This condition can be expressed as

Inserting Eq. (11) into Eq. (13) leads to

By inserting Eq. (10) into Eq. (14), then ${\sigma }_{\text{tilt}}^2$ can be obtained as

When $r>D/2$, set the phase autocorrelation $B_{\mathrm{FFT}}(r)$ to zero

From Eq. (16), it is obvious that the simulated phase screen is only valid in the range of $r\le D/2$. For a rectangular screen with $D_x\times D_y$ sizes, the valid range is $r\le \min (D_x/2,D_y/2)$, so this simulation method is meaningless for the rectangular screens.

We can obtain the discrete power spectrum ${\mathrm{\Phi }}_{\mathrm{FFT}}$ using the above phase autocorrelation $B_{\mathrm{FFT}}$ to replace $B_{\text{theory}}$ in Eq. (7). Because some elements in ${\mathrm{\Phi }}_{\mathrm{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 ${\mathrm{\Phi }}_{\mathrm{FFT}}$ to replace ${\mathrm{\Phi }}_{\text{theory}}$ in Eq. (1), then the FFT-based phase screen ${\varphi }_{\mathrm{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

4.

Evaluation of the Simulation Error

The expected phase structure function of the simulated phase screen is given by

The expected relative simulation error can be defined as

In the range of $r\le D/2$, the simulation error is mainly induced by the negative values in ${\mathrm{\Phi }}_{\mathrm{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$.

Fig. 4

The expected errors as a function of the separation distance $m\mathrm{\Delta }$, $D=2\mathrm{m}$.

5.

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

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 $\mathrm{err}B$, the error is dependent on azimuth and the points with large error are located near the center of matrix $\mathrm{err}B$, so only the elements in the central region of $\mathrm{err}B$ require predistortion.

Fig. 5

The expected errors distribution diagram of the phase autocorrelation. Simulation parameters are $M=2048$, $D=2\mathrm{m}$, $L_0=100\mathrm{m}$, and $r_0=0.2\mathrm{m}$.

A Gaussian shaped predistortion coefficients matrix is suggested as

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.

Fig. 6

The expected residual errors as a function of the separation distance $m\mathrm{\Delta }$ after high frequency compensation, $D=2\mathrm{m}$.

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 (|\mathrm{re}\_\mathrm{err}D|)$ represents the maximum absolute value in the relative error matrix $\mathrm{re}\_\mathrm{err}D$. For $M\le 2048$, the maximum residual error is $<0.13\%$.

Fig. 7

The maximum residual errors after high frequency error compensation.

6.

Simulation Results

Figure 8 shows a sample of the simulated phase screen with the parameters of $D=2\mathrm{m}$, $M=256$, $L_0=20\mathrm{m}$, $r_0=0.2\mathrm{m}$, $A=1.5$, and $W=D/4$.

Fig. 8

Sample of the simulated phase screen.

Figure 9 shows the simulated phase structure function (a) and the error (b) averaged over $2\times {10}^6$ simulated phase screens. When $r\le D/2$ the error is very small, only about 0.1%, but when $r>D/2$, the error increases rapidly.

Fig. 9

The phase structure function (a) and the error (b) in diagonal direction. Other parameters are $M=512$, $A=1.5$, and $W=D/4$.

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.

7.

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.