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,2–8 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 ratios. Results from his analysis show that these FFT-based simulations are accurate for large screen size () to outer scale parameter () ratios. For a screen size of and outer scale of , 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 ratios, even after compensating with modified subharmonics.
In Fig. 1, we show11 a typical situation where the simulation band is actually smaller than full band , where is the sampling size defined as the ratio of screen size to sampling number and 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 . 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 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 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:143 are represented as follows. 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 and and 0.866 for 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 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
Algorithm to Compensate for Residual Error in Phase Structure Function
To calculate the remaining error in the final , Eq. (4) is converted to phase structure matrix with the help of Eq. (1) with the assumption that and are zero because we are not concerned about the piston component. This gives the following equation3 given as follows:
We need to compensate so that is minimized. However, simply adding error correction terms in the 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 .
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 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 , , and all are dependent on only and hence are center symmetric functions. We call the fitted Gaussian phase structure matrix and the corresponding Gaussian phase autocorrelaiton matrix [using Eq. (1)].
The final equation for can then be written as
Here, we have used from Eq. (1). A look at the power spectrum of in Fig. 2 shows that it contains negative terms7 for the case of . Directly putting those frequency terms equal to zero leads to a loss in the energy spectrum. Hence, 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 . The tip/tilt component from phase autocorrelation matrix is given as77
After setting, , the remaining phase autocorrelation matrix is given as follows:
The power spectra and of the phase autocorrelation matrices and are obtained by standard Fourier transformation. Figure 3, shows the remaining negative power elements present in the power spectrum of 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 equal to zero by hand. The new error matrix is given as
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 to are given below.
fmincon: High Frequency error Optimization
Table 1 shows the result of curve fitting using cftool for different cases of and N, which demonstrates that the Gaussian error matrix can compensate for a wide range of ratios and under different sampling constraints.
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 , say for a large future telescope, , and median value of .
The output from the above algorithm corresponding to minimum error entry as in (step ), has been plotted against and . Figure 5 shows a 3D rendering of matrix with a maximum separation of up to 40 m, corresponding to Eq. (4).
Figure 6 shows 1D matrix (radial section from 3D matrix) along with 1D fitted curve including the extrapolated part, for a maximum separation of up to 40 m.
Lastly, MRE values are stored against 500 entries corresponding to ranging from 1 to 10, ranging from 1 to 10 and 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 and entry with minimum MRE has been picked up. Typical time required to perform each iteration for this case is on 2.3 GHz quad-core Intel Core i5 Macbook Pro 2018 model.
Validation via Phase Structure Function Calculated from Simulated Phase Screen
To obtain the phase screen from the power spectrum, the following relation is used:7
For validation, we consider scenarios of apertures up to 40 m, i.e., , at a median for two different sampling levels 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:
Here, is the phase structure function from the simulated phase screen. The magnitude of the peak relative error is for and for as shown in Fig. 8.
We also illustrate the performance (shown in Fig. 9) with parameters , , and 1000 m, , , which cover the extreme cases (very low ratios) which leads to the maximum error in the simulation. The magnitude of the peak relative error is for and for . Figure 10 shows one realization of the corresponding phase screen plots for .
Figure 11 contains results of magnitude of the peak relative error in for the case of different sampling points , for ranges up to 1024 m, and . Similarly, Fig. 12 contains results of magnitude of the peak relative error in for the case of different sampling points , for ranges up to 100 m, , and .
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 . The reason for this stems from the non-zero value of at , where is formed from the removal of piston and tilt from . This can be resolved using a better smoothening operator, which we can multiply with 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 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 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 ratios. The method we propose is independent of the ratio choice. However, we emphasize that properly sampling and the high-frequency phase spectrum forces to be at least larger than and preferably up to the inner scale limit . 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 , , , and 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 , , , and 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 to 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 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 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 ratios. The accuracy of our method from low-frequency to the high-frequency range is between 0.5% and 3% for 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.
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.