Open Access
17 June 2021 Generalized approach to compensate for low- and high-frequency errors in fast Fourier transform-based phase screen simulations
Sorabh Chhabra, Jyotirmay Paul, Anamparambu N. Ramaprakash, Avinash Surendran
Author Affiliations +

Fast Fourier transform-based phase screen simulations give accurate results only when the screen size (G) is much larger than the outer scale parameter (L0). Otherwise, they fall short in correctly predicting both the low and high frequency behaviors of turbulence-induced phase distortions. Subharmonic compensation is a commonly used technique that aids in low-frequency correction but does not solve the problem for all values of screen size to outer scale parameter ratios (G  /  L0). A subharmonics-based approach will lead to unequal sampling or weights calculation for subharmonics addition at the low-frequency range and patch normalization factor. We have modified the subharmonics-based approach by introducing a Gaussian phase autocorrelation matrix that compensates for these shortfalls. We show that the maximum relative error in structure function with respect to theoretical value is as small as 0.5% to 3% for (G  /  L0) ratio of 1/1000 even for screen sizes up to 100 m diameter.



Accurately simulating the atmospheric turbulence behavior is well recognized as very challenging. For a variety of purposes such as the design and development of adaptive optics systems, speckle imaging techniques, and atmospheric propagation studies, it is essential to simulate good atmospheric phase screen models. Methods based on Zernike polynomial expansions,1 fast Fourier transform (FFT)-based methods,28 and low-frequency optimization method9 have been in use for this purpose. The Zernike polynomial method, which is widely in use, has a limitation due to the maximum number of coefficients needed for accurate compensation. The optimization method, which compensates accurately for low-frequency part of the spectrum using unequal sampling and unequal weight in low-frequency region, does not cover high-frequency deficiencies. Among these, FFT-based methods are computer memory size friendly and widely accepted. But, FFT operators assume uniform sampling for the non-uniformly distributed phase power spectrum, which can lead to underestimation in the low- and high-frequency out of band regions, as shown in Fig. 1. Thus, it has limitations in recreating the true phase power spectrum. To compensate for low-frequency components, Johansson and Gavel3 suggested employing the modified subharmonics equation (an adaptation from Lane et al.10), which works well up to an infinite outer scale length. Sedmak6 later compared the performance of this method with that of Lane et al.10 by actually calculating the phase structure function from the simulated screen. He improved upon Lane et al.10 by employing different fine-tuned subharmonic weights for different G/L0 ratios. Results from his analysis show that these FFT-based simulations are accurate for large screen size (G) to outer scale parameter (L0) ratios. For a screen size of G=200  m and outer scale of L0=25  m, the maximum relative error in the simulation approaches 1%. Our simulations demonstrate that the errors from low-frequency components start shooting up once we move to smaller G/L0 ratios, even after compensating with modified subharmonics.

Fig. 1

Comparison between simulation band and full band.


In Fig. 1, we show11 a typical situation where the simulation band (1G1Δ) is actually smaller than full band (1L01l0), where Δ is the sampling size defined as the ratio of screen size G to sampling number N and l0 is the inner scale parameter. In practice, the simulations are often curtailed at the low-frequency end, to a few times the optical beam size (say as determined by the telescope or laser beam diameter), whereas at the high frequency end, they often extend to only a few times that determined by the Fried parameter r0. Clearly, the larger the simulation band to full band ratio, the more accurate the simulated results will be.

On the one hand, the apertures of upcoming and future astronomical telescopes are often of the same order or even larger than the typical median outer scale sizes of about 20 to 25 m.12 On the other hand, wavefront sensing and compensation technologies are fast progressing that Nyquist sampling at r0 scales even for large aperture telescopes are becoming quite possible. Thus, atmospheric turbulence simulations have to deal with a wide range in a multi-dimensional parameter space.

For working with very small apertures relative to the outer scale, it may appear that we need to simulate only a relatively small screen size. But cutting out small apertures from a larger screen introduces deviation from phase structure function due to misrepresentation of low-frequency components present in the small screen power spectrum.

In this paper, we present an approach and a corresponding algorithm to deal with phase screen simulations for a wide range of G/L0 ratios, using the FFT-based method.13 Our technique builds upon the modified subharmonic approach of Johansson and Gavel3 and is inspired by Jingsong Xiang’s work.7 It works well for space- and time-invariant, zero intermittency atmospheric turbulence. Section 2 explains how to obtain phase autocorrealtion matrix using phase power spectrum, Sec. 3 presents the algorithm part to compensate for the remaining error in phase structure function calculation, Sec. 4 steps through the implementation of the algorithm with the help of a flow chart, Sec. 5 covers the validation of the technique using results from simulated phase screens, and Sec. 6 provides the concluding remarks.


Obtaining Phase Autocorrelation Matrix Using Phase Power Spectrum

The 2D phase structure function and phase autocorrelation matrix are related as follows:14

Eq. (1)

where Bϕ(m,n) is the phase autocorrelation matrix and (m,n) is the coordinates along x and y axis. The 2D phase autocorrelation matrices for the FFT-based phase screen and the modified subharmonic method by Johansson and Gavel3 are represented as follows.

Eq. (2)


Eq. (3)

where fFFT2(m,n) and fSUB2(m,n) are the von-Kármán spectrum and subharmonic power spectrum as explained by Johansson and Gavel. (Nx,Ny) are sample points, p is the p’th subharmonic, and Np is the total number of subharmonics. Set fFFT=0, for (m,n)=(0,0) and fSUB=0, for (m,n)=(1,0) and (0,1) as originally proposed by Lane et al.10 There will be an overlap between subharmonic energy sample and secondary lobes from first sample of high-frequency spectrum or harmonic sample during subharmonic addition. Earlier this leakage of energy has been dealt using patch normalization factor, where first patch of high-frequency spectrum is weighted by 0.707 for (m,n)=(±1,0) and (m,n)=(0,±1) and 0.866 for (m,n)=(±1,±1) in the original method of Johansson and Gavel.3 Similarly, the original method of Lane et al.,10 Sedmak6 proposed the corresponding weights to be 0.935 and 0.998, respectively. Our simulations show that these weights do not fit perfectly for different G/L0 ratios and hence need to be tuned on a case by case basis. We have made our approach independent from these weights assignments. The weight factor has been set equal to 1 in our approach. Section 3 explains this approach in detail.

The 2D phase autocorrelation matrix after compensating with subharmonics is represented as

Eq. (4)



Algorithm to Compensate for Residual Error in Phase Structure Function

To calculate the remaining error in the final Bϕ(m,n), Eq. (4) is converted to phase structure matrix Dϕ(m,n) with the help of Eq. (1) with the assumption that BϕFFT(0,0) and BϕSUB(0,0) are zero because we are not concerned about the piston component. This gives the following equation

Eq. (5)

where Dtheory(m,n) is the well-known theoretical von-Kármán phase structure matrix,3 given as follows:

Eq. (6)

where r2=(mΔ)2+(nΔ)2, Δ=G/N.

We need to compensate Dϕ so that Derror is minimized. However, simply adding error correction terms in the Dϕ matrix directly would only introduce further error into the system while taking the Fourier transform. This is because any matrix or curve in general will have higher order moments. Thus, if we take the Fourier transform of the adjusted equation, the resultant curve will have completely different moments and hence power spectrum. This is because the transition between two steps in the error matrix will not be smooth, which introduces additional errors due to Gibb’s phenomena such as overshoots. Just curve fitting with any function does not satisfy the additional requirement of leaving the power spectrum unaffected by the process. What we really need is to introduce a smoothening operator such as a Gaussian function in the phase autocorrelation matrix, which exactly compensates for Derror.

For that we have developed an iterative algorithm (see the flow chart shown in Fig. 4) and implemented it in Matlab. The algorithm looks for the perfect Gaussian curve that minimizes the Derror matrix. We use Matlab cftool to initially determine the correct 1D Gaussian matrix and later convert it into a 2D matrix by exploiting the fact that Bϕ(r), Btheory(r), and BSUB(r) all are dependent on r only and hence are center symmetric functions. We call the fitted Gaussian phase structure matrix Dgauss and the corresponding Gaussian phase autocorrelaiton matrix Bgauss [using Eq. (1)].

The final equation for Btot can then be written as

Eq. (7)


Here, we have used Bgauss=Dgauss/2 from Eq. (1). A look at the power spectrum of Btot(m,n) in Fig. 2 shows that it contains negative terms7 for the case of G/L0<1. Directly putting those frequency terms equal to zero leads to a loss in the energy spectrum. Hence, Btot(m,n) matrix needs to be preprocessed to eliminate most of these negative values in the power spectrum. Over small frequencies, piston and tip/tilt components account for most of these high magnitude negative elements. Therefore, we first extract the piston and tip/tilt components from the phase autocorrelation matrix Btot. The tip/tilt component from phase autocorrelation matrix is given as7

Eq. (8)

where σtilt2 is the variance of the random tilt angle in the x or y directions and given as follows:7

Eq. (9)


Fig. 2

Negative power spectrum values for small G/L0 ratios.


After setting, Btilt(0)=0, the remaining phase autocorrelation matrix is given as follows:

Eq. (10)


The power spectra fhigh2 and ftilt2 of the phase autocorrelation matrices Bhigh(r) and Btilt(r) are obtained by standard Fourier transformation. Figure 3, shows the remaining negative power elements present in the power spectrum of Bhigh matrix. In comparison to Fig. 2, the largest negative power contributions fall by factor of three order of magnitude. Now, we set the negative values in fhigh2 equal to zero by hand. The new error matrix is given as

Eq. (11)

where Bhigh(r) is the phase autocorrelation matrix obtained after setting the negative elements in fhigh2 to zero. The residual error that is present in the high-frequency region can then be reduced with the help of a Gaussian smoothing operator, using Matlab fmincon tool. The high-frequency compensated matrix is given as

Eq. (12)

where Hhighcomp(r) is the smoothening operator, multiplied with error matrix to reduce the high frequency errors. fmincon gives the optimized parameter for smoothening operator by calculating the final error in the Dϕ(r) matrix with respect to Dtheory(r).

Fig. 3

Residual negative power spectrum values after removing tip/tilt from Btot for small G/L0 ratios.



Implementation of the Compensation Algorithm

In this section, we explain the error compensation algorithm with the help of the flow chart shown in Fig. 4. Brief explanations of each of the steps from L1 to L12 are given below.

  • L1: Input screen size G, outer scale size L0, Fried parameter r0, and number of samples N.

  • L2: Initialize the algorithm with Np, the total number of subharmonics, and extrapolation factor EF, both ranging from 1 to 10. The EF factor is relevant while performing curve fitting, e.g., EF=2 means curve fitting will work from 3 to N/2 points and later that curve will be extrapolated from 1 to N/2 points (example shown in Sec. 4.1 and Fig. 6)

  • L3: Obtain BϕFFT and BϕSUB based upon L1 and L2 parameters.

  • L4: Add the matrices that were calculated in L3 layer, call that Bϕ.

  • L5: Obtain Dϕ from Bϕ using Eq. (1) and also produce Dtheory matrix based upon L1 and L2. Then obtain error matrix Derror using Eq. (5).

  • L6: Obtain 1D array from Derror matrix from the center and give as input to the curve fitting tool cftool, which works on 1D data. The output of the tool will be a best fitted curve in terms of sum of Gaussian’s (SOG), called Dgauss.

  • L7: The final expression for Btot is shown Eq. (7).

  • L8: Extract tilt component from Btot matrix using Eqs. (8) and (9).

Fig. 4

Flow chart for error compensation.


fmincon: High Frequency error Optimization

  • L9L11: Obtain error matrix Bhigherr after setting negative elements in the power spectrum to zero. Compensate for high frequency error by multiplying error matrix with smoothening operator-SOG. Calculate the maximum remaining error in structure function matrix relative to the Dtheory matrix. Thus fmincon will give parameters for the smoothing operator, which gives the lowest possible residual error.

  • L12: Update the entry for Np and EF to next value, evaluation from L2L12 would go in loop, and maximum relative error (MRE) value stored in vector form. At the end of the iterations, min entry will get extracted out from stored vector and accepted for final analysis.

Table 1 shows the result of curve fitting using cftool for different cases of G/L0 and N, which demonstrates that the Gaussian error matrix can compensate for a wide range of G/L0 ratios and under different sampling constraints.

Table 1

Result of curve fitting against Gaussian function for different cases of G/L0 and N in terms of MRE for fixed r0=0.2  m.




To illustrate the robustness of the above algorithm, we have taken an example with G=80  m, say for a large future telescope, N=256, and median value of L0=20  m.

The output from the above algorithm corresponding to minimum error entry as in (step L12), has been plotted against EF=5 and Np=8. Figure 5 shows a 3D rendering of Derror matrix with a maximum separation of up to 40 m, corresponding to Eq. (4).

Fig. 5

3D Derror matrix for case G=80  m, L0=20  m, r0=0.2  m, and Nx=Ny=256.


Figure 6 shows 1D Derror matrix (radial section from 3D Derror matrix) along with 1D fitted curve Dgauss including the extrapolated part, for a maximum separation of up to 40 m.

Fig. 6

1D Derror matrix fitted against Dgauss matrix, along with an extrapolated part of the curve. Here G=80  m, L0=20  m, r0=0.2  m, and Nx=Ny=256.


Lastly, MRE values are stored against 500 entries corresponding to Np ranging from 1 to 10, EF ranging from 1 to 10 and SoG ranging from 2 to 6 after performing cftool fitting. This has been arranged in descending order and shown in Fig. 7, which illustrates a large set of iterations where errors are <1% and entry with minimum MRE has been picked up. Typical time required to perform each iteration for this case is 4.9  s on 2.3 GHz quad-core Intel Core i5 Macbook Pro 2018 model.

Fig. 7

Maximum relative error MRE with the maximum number of iterations for G=80  m, L0=20  m, r0=0.2  m, and Nx=Ny=256.



Validation via Phase Structure Function Calculated from Simulated Phase Screen

To obtain the phase screen ϕ(m,n) from the power spectrum, the following relation is used:7

Eq. (13)

where Ra(m,n) and Rb(m,n) are zero-mean and unity-variance Gaussian random number generator. We get ϕhigh and ϕtilt, by replacing f with fhigh and ftilt, which are square roots of the power spectrum corresponding to autocorrelation matrix Bhighcomp and Btilt, respectively.

For validation, we consider scenarios of apertures up to 40 m, i.e., G=80  m, at a median L0=20  m for two different sampling levels N=256 and 512. The phase structure function, defined as an ensemble average of differences of phases at various separation,14 has been averaged over 100 K independent frames. The relative error in phase structure function is calculated as follows:

Eq. (14)


Here, Dϕsim(r) is the phase structure function from the simulated phase screen. The magnitude of the peak relative error max(|err(r)|) is <1.6% for N=256 and <0.5% for N=512 as shown in Fig. 8.

Fig. 8

Left: Compares simulated structure function with respect to theoretical structure function for maximum separation of G/2 for two different cases N=256 and 512, for fixed G=80  m, r0=0.2  m, and L0=20  m. Right: Calculates the magnitude of relative error in simulated structure function for maximum separation of G/2, for both cases.


We also illustrate the performance (shown in Fig. 9) with parameters G=1  m, L0=100  m, and 1000 m, N=128, r0=0.2  m, which cover the extreme cases (very low G/L0 ratios) which leads to the maximum error in the simulation. The magnitude of the peak relative error max(|err(r)|) is <1.6% for L0=100  m and <1.8% for L0=1000  m. Figure 10 shows one realization of the corresponding phase screen plots for L0=100  m.

Fig. 9

(a) Compares simulated structure function with respect to theoretical structure function for maximum separation of G/2 for two different cases L0=100 and 1000 m, for fixed G=1  m, r0=0.2  m, and N=128. (b) Calculates the magnitude of relative error in simulated structure function for maximum separation of G/2, for both cases.


Fig. 10

Phase screen for case G=1  m, L0=100  m, r0=0.2  m, and Nx=Ny=128.


Figure 11 contains results of magnitude of the peak relative error in Dϕsim(r) for the case of different sampling points N=128/256/512/1024, for L0 ranges up to 1024 m, r0=0.2  m and G=2  m. Similarly, Fig. 12 contains results of magnitude of the peak relative error in Dϕsim(r) for the case of different sampling points N=128/256/512/1024, for G ranges up to 100 m, r0=0.2  m, and L0=25  m.

Fig. 11

The magnitude of the peak relative error for N=128, 256, 512, and 1024 for L0 ranges up to 1024 m. Here G=2  m and r0=0.2  m.


Fig. 12

The magnitude of the peak relative error for N=128, 256, 512, and 1024 for screen size up to 100 m. Here, L0=25  m and r0=0.2  m.


There are some outliers that have a high residual error as shown in Figs. 11 and 12, because we have not set the phase autocorrelation matrix to zero for r>G/2. The reason for this stems from the non-zero value of Bhigh(r) at r>G/2, where Bhigh(r) is formed from the removal of piston and tilt from Btot(r). This can be resolved using a better smoothening operator, which we can multiply with Bhigh so that it falls to zero progressively and not sharply. This can provide further improvement in the compensation.



In this paper, we put forward a new method to compensate for the residual error in both the low and/or high-frequency region of FFT simulated phase screens that remain even after compensating with the modified subharmonic method. This method provides accurate phase screen structure for even G/L0 ratios as small as 1/1000 plus screen sizes as large as 100 m. No patch normalization factor is needed, no need to calculate subharmonic weight coefficient10 and weights to compensate for high-frequency components, as done by Sedmak.6 While adequately large G/L0 ratios may be the natural choice for modern large telescopes, simulations that deal with applications such as laser beam propagation through turbulent atmospheres would tend to have very small G/L0 ratios. The method we propose is independent of the G/L0 ratio choice. However, we emphasize that properly sampling r0 and the high-frequency phase spectrum forces N to be at least larger than (2G/r0) and preferably up to the inner scale limit (2G/l0). Currently, we have demonstrated this technique for only circular screens. We have used a GPU processor with total number of 128 cores, such that each iteration runs independently on each core. We have fixed the number of iterations to 500, although increasing this will lead to improvement of errors in some cases. Each core takes about 0.06, 0.1, 0.36, and 1.1  min for sampling sizes of 128, 256, 512, and 1024. On the above GPU system, this translates to total computing times for error minimization of about 0.2, 0.5, 1.25, and 4  min for sampling sizes of 128, 256, 512, and 1024, respectively. Once the coefficients are determined, generating multiple phase screen realizations from the corresponding power spectrum takes a few milliseconds at most on 2.3 GHz quad-core Intel Core i5 Macbook Pro 2018 model. Then it takes <1  min to 10  min for averaging over 100k phase screens, for sampling sizes ranging from 128 to 1024.

The uniqueness of our approach is its ability to deal with any G/L0 ratio within a very broad range, in an automated iterative process with little human intervention needed for tuning of parameters. Any standard FFT-based approach (say Sedmak’s6 compensated approach) for a given computer platform is computationally fast, only if we already have determined proper measures of the various compensating components such as the patch normalization factor, subharmonic weights, and high-frequency weights. Typically, determining these compensations is where the difficulty is due to mathematical complexity, algorithmic limitations, and/or computational power requirements. Our algorithm accomplishes the determination of the required compensation in very little time, with fairly reasonable computational power while at the same time keeping the residual errors competitively low using an appropriate compensator. Other existing FFT based approaches have limitations in their operable G/L0 range. For example, Xiang et al.7 offer a very computationally fast approach but does not apply subharmonic compensation. Zhang et al.9 does not consider compensation for the high-frequency error, thus leaving a residual error of more than 100% in the high-frequency region. Sedmak’s6 approach needs the determination of accurate subharmonic weights for different G/L0 ratios. The accuracy of our method from low-frequency to the high-frequency range is between 0.5% and 3% for G/L0 as low as 1/1000 and screen size up to 100 m in diameter.


We would like to thank Sedmak for providing insights into the nature of atmospheric phase power spectrum through private communication. We also thank Xiang for sharing his MATLAB code that calculates the phase structure function quickly for a large number of phase screens. We acknowledge usage of IUCAA’s Pegasus cluster computer for running multiple independent iterations in parallel.



N. A. Roddier, “Atmospheric wavefront simulation and Zernike polynomials,” Proc. SPIE, 1237 668 –679 (1990). PSISDG 0277-786X Google Scholar


B. J. Herman and L. A. Strugala, “Method for inclusion of low-frequency contributions in numerical representation of atmospheric turbulence,” Proc. SPIE, 1221 183 –192 (1990). PSISDG 0277-786X Google Scholar


E. M. Johansson and D. T. Gavel, “Simulation of stellar speckle imaging,” Proc. SPIE, 2200 372 –383 (1994). PSISDG 0277-786X Google Scholar


G. Sedmak, “Performance analysis of and compensation for aspect-ratio effects of fast-Fourier-transform-based simulations of large atmospheric wave fronts,” Appl. Opt., 37 (21), 4605 –4613 (1998). APOPAI 0003-6935 Google Scholar


B. L. McGlamery, “Computer simulation studies of compensation of turbulence degraded images,” Proc. SPIE, 0074 225 –233 (1976). PSISDG 0277-786X Google Scholar


G. Sedmak, “Implementation of fast-Fourier-transform-based simulations of extra-large atmospheric phase and scintillation screens,” Appl. Opt., 43 (23), 4527 –4538 (2004). APOPAI 0003-6935 Google Scholar


J. Xiang, “Fast and accurate simulation of the turbulent phase screen using fast Fourier transform,” Opt. Eng., 53 (1), 016110 (2014). Google Scholar


J. Xiang, “Accurate compensation of the low-frequency components for the FFT-based turbulent phase screen,” Opt. Express, 20 (1), 681 –687 (2012). OPEXFF 1094-4087 Google Scholar


D. Zhang et al., “Accurate simulation of turbulent phase screen using optimization method,” Optik, 178 1023 –1028 (2019). OTIKAJ 0030-4026 Google Scholar


R. Lane et al., “Simulation of a Kolmogorov phase screen,” Waves Random Media, 2 (3), 209 –224 (1992). WRMEEV 0959-7174 Google Scholar


G. Sedmak, (2019). Google Scholar


A. Ziad, “Review of the outer scale of the atmospheric turbulence,” Proc. SPIE, 9909 99091K (2016). PSISDG 0277-786X Google Scholar


S. Chhabra et al., “Gaussian phase autocorrelation as an accurate compensator for FFT-based atmospheric phase screen simulations,” Proc. SPIE, 11448 114487U (2020). PSISDG 0277-786X Google Scholar


F. Roddier, Adaptive Optics in Astronomy, Cambridge University Press(1999). Google Scholar


Sorabh Chhabra is a PhD student at Inter University Center for Astronomy and Astrophysics (IUCAA), Pune. He received his B. Tech degree in electronics and communication from Delhi Technological University (formally Delhi College of Engineering) in 2016 and joined for his PhD in the same year in the Instrumentation Department at IUCAA.

Biographies of the other authors are not available.

© 2021 Society of Photo-Optical Instrumentation Engineers (SPIE)
Sorabh Chhabra, Jyotirmay Paul, Anamparambu N. Ramaprakash, and Avinash Surendran "Generalized approach to compensate for low- and high-frequency errors in fast Fourier transform-based phase screen simulations," Journal of Astronomical Telescopes, Instruments, and Systems 7(2), 025007 (17 June 2021).
Received: 25 January 2021; Accepted: 1 June 2021; Published: 17 June 2021 Logo
Cited by 2 scholarly publications.
Get copyright permission  Get copyright permission on Copyright Marketplace
Fourier transforms




Atmospheric propagation

Atmospheric turbulence

Device simulation


Back to Top