20 January 2014 Fast and accurate simulation of the turbulent phase screen using fast Fourier transform
Author Affiliations +
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.
Xiang: Fast and accurate simulation of the turbulent phase screen using fast Fourier transform

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

(1)

ϕFFT(mΔ,nΔ)=m=M/2M/21n=M/2M/21[Ra(m,n)+iRb(m,n)]Φtheory(mΔ,nΔ)exp[i2π(mm+nn)/M],

(2)

Φtheory(mΔ,nΔ)=0.49r053[(mΔ)2+(nΔ)2+(2πL0)2]116(Δ)2,
where m, n, m, n=M/2,,M/21 are the integer indices, Ra(m,n) and Rb(m,n) are the independent Gaussian random numbers with zero-mean and unit-variance, i=1 denotes the imaginary unit, Φtheory(mΔ,nΔ) is the discrete power spectrum of turbulence, Δ=2π/D and Δ=D/M are the sampling intervals in the spatial frequencies and spatial domains, respectively, L0 is the turbulence outer scale, and r0 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 methods4,8,10 or weighted subharmonic methods11,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 times compared to the standard FFT-based method and the compensation processes are complicated. Charnotskii13 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.

Xiang14 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

(3)

Dtheory(r)=2[Btheory(0)Btheory(r)],
where Btheory(r) is the theoretical phase autocorrelation and r is the separation distance between two points.

The theoretical phase autocorrelation function can be expressed as7

(4)

Btheory(r)=(L0/r0)5/325/6π8/3Γ(11/6)·[(245)Γ(65)]5/6(2πrL0)5/6K5/6(2πrL0),forr>0,
where K5/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

(5)

Btheory(0)=(L0/r0)5/325/6π8/3Γ(11/6)[(245)Γ(65)]5/621/6Γ(5/6).

For the FFT-based phase screen, the expected phase autocorrelation matrix is an inverse discrete Fourier transform (IDFT) of the power spectrum matrix8

(6)

Bexp(mΔ,nΔ)=m=M/2M/21n=M/2M/21Φtheory(mΔ,mΔ)exp[i2π(mm+nn)/M].

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

(7)

ΦFFT(mΔ,nΔ)=1M2m=M/2M/21n=M/2M/21Btheory[(m2+n2)Δ]exp[i2π(mm+nn)/M].

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 Bexp will be consistent with the theoretical phase autocorrelation Btheory, 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 L0/D<1/2, the simulation error is very small. For the case of large L0/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 105 phase screens. For the case of L0=1m, the simulated curve agrees well with the theoretical curve, it is hard to distinguish them.

OE_53_1_016110_f001.png

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.

Fig. 2

Comparison of power spectrum density given by Eqs. (7) and (2). Simulation parameters are D=2m, M=256, r0=0.2m, and L0=1m.

OE_53_1_016110_f002.png

Theoretically, the power spectrum should be real and non-negative. 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.

Fig. 3

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

OE_53_1_016110_f003.png

3.

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 L0/D is very similar to that with a small L0/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

(8)

φtilt(x,y)=(θxx+θyy)σtilt,
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 x- or y-directions.

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

(9)

Dtilt(r)=[ϕtilt(x1,y1)ϕtilt(x2,y2)]2=r2σtilt2,

(10)

Btilt(r)=Btilt(0)r2σtilt2/2,
where Btilt(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

(11)

BFFT(r)=Btheory(r)Btilt(r)Btheory(D/2)+Btilt(D/2).

For the case of infinite turbulence outer scale (L0), the phase autocorrelation function Btheory cannot be obtained by Eq. (4). In this case, Btheory can be obtained by

(12)

Btheroy(r)=Btheroy(0)6.88(r/r0)5/3/2,
where Btheroy(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 BFFT(r) is zero at r=D/2. This condition can be expressed as

(13)

[BFFT(D/2)BFFT(D/2ε)]/ε=0,
where ε is a very short distance <Δ/100.

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

(14)

Btheory(D/2)Btheory(D/2ε)=Btilt(D/2)Btilt(D/2ε).

By inserting Eq. (10) into Eq. (14), then σtilt2 can be obtained as

(15)

σtilt2=[Btheory(D/2ε)Btheory(D/2)]ε(Dε)/2.

When r>D/2, set the phase autocorrelation BFFT(r) to zero

(16)

BFFT(r)=0,forr>D/2.

From Eq. (16), it is obvious that the simulated phase screen is only valid in the range of rD/2. For a rectangular screen with Dx×Dy sizes, the valid range is rmin(Dx/2,Dy/2), so this simulation method is meaningless for the rectangular screens.

We can obtain the discrete power spectrum ΦFFT using the above phase autocorrelation BFFT to replace Btheory 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

(17)

ΦFFT(mΔ,mΔ)=0,forΦFFT(mΔ,mΔ)<0.

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

(18)

ϕ(mΔ,nΔ)=ϕFFT(mΔ,nΔ)+ϕtilt(mΔ,nΔ).

4.

Evaluation of the Simulation Error

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

(19)

Dϕexp(mΔ,nΔ)=DFFTexp(mΔ,nΔ)+Dtilt[(m2+n2)Δ],
where DFFTexp is the expected phase structure function of the FFT-based phase screen ϕFFT, which can be expressed as

(20)

DFFTexp(mΔ,nΔ)=2[BFFTexp(0,0)BFFTexp(mΔ,nΔ)],
where BFFTexp is given by8

(21)

BFFTexp(mΔ,nΔ)=m=M/2M/21n=M/2M/21ΦFFT(mΔ,mΔ)exp[i2π(mm+nn)/M].

The expected relative simulation error can be defined as

(22)

re_errD(mΔ,nΔ)=Dϕexp(mΔ,nΔ)/Dtheory[(m2+n2)Δ]1.

In the range of rD/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 L0/D. When L0/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 r0.

Fig. 4

The expected errors as a function of the separation distance mΔ, D=2m.

OE_53_1_016110_f004.png

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

(23)

errB(mΔ,nΔ)=BFFTexp(mΔ,nΔ)Btheory[(m2+n2)Δ].

After predistortion, the phase autocorrelation can be expressed as

(24)

Bhighcom(mΔ,nΔ)=Btheory[(m2+n2)Δ]C(mΔ,nΔ)errB(mΔ,nΔ),
where C(mΔ,nΔ) is the predistortion coefficient. The discrete power spectrum ΦFFT after high frequency error compensation can be obtained by using Bhighcom to replace Btheory 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.

Fig. 5

The expected errors distribution diagram of the phase autocorrelation. Simulation parameters are M=2048, D=2m, L0=100m, and r0=0.2m.

OE_53_1_016110_f005.png

A Gaussian shaped predistortion coefficients matrix is suggested as

(25)

C(mΔ,nΔ)=Aexp{[(mΔ)2+(nΔ)2]/W2}.

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Δ after high frequency compensation, D=2m.

OE_53_1_016110_f006.png

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(|re_errD|) represents the maximum absolute value in the relative error matrix re_errD. For M2048, the maximum residual error is <0.13%.

Fig. 7

The maximum residual errors after high frequency error compensation.

OE_53_1_016110_f007.png

6.

Simulation Results

Figure 8 shows a sample of the simulated phase screen with the parameters of D=2m, M=256, L0=20m, r0=0.2m, A=1.5, and W=D/4.

Fig. 8

Sample of the simulated phase screen.

OE_53_1_016110_f008.png

Figure 9 shows the simulated phase structure function (a) and the error (b) averaged over 2×106 simulated phase screens. When rD/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.

OE_53_1_016110_f009.png

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 l0 cannot be easily simulated. The MATLAB code is be available via e-mail.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 61071117 and 61275077) and Scientific and Technological Research Program of Chongqing Municipal Education Commission (Grant No. KJ130515).

References

1. B. L. McGlamery, “Computer simulation studies of compensation of turbulence degraded images,” Proc. SPIE 74, 225–233 (1976).PSISDG0277-786X http://dx.doi.org/10.1117/12.954724 Google Scholar

2. N. Roddier, “Atmospheric wavefront simulation using Zernike polynomials,” Opt. Eng. 29(10), 1174–1180 (1990).OPEGAR0091-3286 http://dx.doi.org/10.1117/12.55712 Google Scholar

3. K. A. Winick, “Atmospheric turbulence-induced signal fades on optical heterodyne communication links,” Appl. Opt. 25(11), 1817–1825 (1986).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.25.001817 Google Scholar

4. R. G. LaneA. GlindemannJ. C. Dainty, “Simulation of a Kolmogorov phase screen,” Waves Random Media 2(3), 209–224 (1992).WRMEEV0959-7174 http://dx.doi.org/10.1088/0959-7174/2/3/003 Google Scholar

5. C. M. HardingR. A. JohnstonR. G. Lane, “Fast simulation of a Kolmogorov phase screen,” Appl. Opt. 38(11), 2161–2170 (1999).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.38.002161 Google Scholar

6. H.-L. Wuet al., “Statistical interpolation method of turbulent phase screen,” Opt. Express 17(17), 14649–14664 (2009).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.17.014649 Google Scholar

7. F. AssématR. W. WilsonE. Gendron, “Method for simulating infinitely long and non stationary phase screens with optimized memory storage,” Opt. Express 14(3), 988–999 (2006).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.14.000988 Google Scholar

8. E. M. JohanssonD. T. Gavel, “Simulation of stellar speckle imaging,” Proc. SPIE 2200, 372–383 (1994).PSISDG0277-786X http://dx.doi.org/10.1117/12.177254 Google Scholar

9. R. Frehlich, “Simulation of laser propagation in a turbulent atmosphere,” Appl. Opt. 39(3), 393–397 (2000).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.39.000393 Google Scholar

10. B. J. HermanL. A. Strugala, “Method for inclusion of low-frequency contributions in numerical representation of atmospheric turbulence,” Proc. SPIE 1221, 183–192 (1990).PSISDG0277-786X http://dx.doi.org/10.1117/12.18342 Google Scholar

11. G. Sedmak, “Implementation of fast-Fourier-transform-based simulations of extra-large atmospheric phase and scintillation screens,” Appl. Opt. 43(23), 4527–4538 (2004).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.43.004527 Google Scholar

12. J. RecolonsF. Dios, “Accurate calculation of phase screens for the modelling of laser beam propagation through atmospheric turbulence,” Proc. SPIE 5891, 589107 (2005).PSISDG0277-786X http://dx.doi.org/10.1117/12.618476 Google Scholar

13. M. Charnotskii, “Sparse spectrum model for a turbulent phase,” J. Opt. Soc. Am. A 30(3), 479–488 (2013).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.30.000479 Google Scholar

14. J. Xiang, “Accurate compensation of the low-frequency components for the FFT-based turbulent phase screen,” Opt. Express 20(1), 681–687 (2012).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.20.000681 Google Scholar

Biography

Jingsong Xiang received his PhD in optics from University of Electronic Science and Technology of China, Chengdu, China, in 2007, an MS degree in optical engineering in 1999, and a BS degree in precision instrumentation in 1996, both from University of Chongqing, Chongqing, China. Currently, he is an associate professor at Chongqing University of Posts and Telecommunications, Chongqing, China. His research interests include free-space optical communication, atmospheric turbulence, and optical fiber communication.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Jingsong Xiang, "Fast and accurate simulation of the turbulent phase screen using fast Fourier transform," Optical Engineering 53(1), 016110 (20 January 2014). https://doi.org/10.1117/1.OE.53.1.016110
JOURNAL ARTICLE
7 PAGES


SHARE
Back to Top