ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

ISA Transactions journal homepage: www.elsevier.com/locate/isatrans

Nonstationary multiscale turbulence simulation based on local PCA Alessandro Beghi a,1, Angelo Cenedese a,2, Andrea Masiero b,3 a b

Department of Information Engineering, University of Padova, via Gradenigo 6/B, 35131 Padova, Italy Interdepartmental Research Center of Geomatics, University of Padova, via dell0 Università 16, 35020 Legnaro, PD, Italy

art ic l e i nf o

a b s t r a c t

Article history: Received 31 July 2013 Received in revised form 18 September 2013 Accepted 14 December 2013

Turbulence simulation methods are of fundamental importance for evaluating the performance of control strategies for Adaptive Optics (AO) systems. In order to obtain a reliable evaluation of the performance a statistically accurate turbulence simulation method has to be used. This work generalizes a previously proposed method for turbulence simulation based on the use of a multiscale stochastic model. The main contributions of this work are: first, a multiresolution local PCA representation is considered. In typical operating conditions, the computational load for turbulence simulation is reduced approximately by a factor of 4, with respect to the previously proposed method, by means of this PCA representation. Second, thanks to a different low resolution method, based on a moving average model, the wind velocity can be in any direction (not necessarily that of the spatial axes). Finally, this paper extends the simulation procedure to generate, if needed, turbulence samples by using a more general model than that of the frozen flow hypothesis. & 2014 ISA. Published by Elsevier Ltd. All rights reserved.

Keywords: Adaptive optics Multiscale stochastic systems Signal processing Fast Fourier transforms Principal component analysis

1. Introduction The technological advancements of the last few decades opened the way to the use of control systems in a large number of high-end technological applications. In particular, their recent introduction in astronomical and radio telescope applications, [1–5], led to great improvements in the performance of the developed systems. In the specific case of astronomical telescopes, Active and Adaptive Optics systems have been introduced to compensate the undesirable effect of permanent manufacturing errors, long time aberrations, and atmospheric turbulence [6–9]. The Active Optics system intervenes to preserve the mirror optimal shape against static or slowly varying environmental factors, such as manufacturing errors, gravity due to telescope inclination, wind, and thermal effects. On the other hand, Adaptive optics [10–14] operate on a much shorter timescale to compensate for factors that affect the image at faster timescales (1/100th second or even less): these are usually caused by the atmosphere and are not easily corrected with primary mirrors. The light beam wavefront incident to the telescope is corrugated by the fluctuations of the atmosphere refraction index: light beams arriving on the telescope lens are delayed because of random changes in the

E-mail addresses: [email protected] (A. Beghi), [email protected] (A. Cenedese), [email protected] (A. Masiero). 1 Tel.: þ39 0498277626. 2 Tel.: þ39 0498277677. 3 Tel.: þ39 3200269583.

atmosphere refraction index. As a consequence, the real resolution of the telescope is not determined by the telescope lens diameter, but it is limited by the degradation of the image due to the atmosphere refraction index fluctuations. Thus, Adaptive Optics (AO) systems are used in large ground telescopes to estimate and compensate the effect of the atmosphere by properly commanding a set of deformable mirrors. The goal of the control unit in the AO systems is quite challenging due to the large dimension of the system (deformable mirrors of large telescopes are actuated by thousands of piezosystems [16]), the modeling uncertainties and its nonlinearity, and the high working frequencies (so, in practice, the deformable mirrors only work in a feedback control loop). These considerations motivate an accurate testing procedure to compare different control strategies and evaluate their performance. According to this aim, this paper deals with the problem of statistically accurate simulation of the effect of atmospheric turbulence on ground telescope measurements. Models commonly used to simulate the influence of atmospheric turbulence on ground telescope observations are usually based on the frozen flow hypothesis [1]: at short time scales the temporal evolution of the turbulence can be approximated as a rigid shift over the telescope aperture. Based on the frozen flow hypothesis, several approaches have been proposed in the literature for simulating turbulence in such a way to accurately reproduce its spatial statistical characteristics (Kolmogorov or Von Karman statistics [1,15]). Among the wide range of proposed methods it is possible to distinguish two particular categories: methods based on the fast Fourier transform (FFT) [17], and recent methods based on linear dynamic systems [18–20]. The first are

0019-0578/$ - see front matter & 2014 ISA. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.isatra.2013.12.026

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

attractive because of their high accuracy, however their memory and computational requirements for long simulations are prohibitive. Then, much of the interest in the latter is due to their ability in producing long, possibly infinite, sequences of phase screens. Unfortunately, the computational complexity of methods [18–20] increases (approximately) quadratically with the turbulence resolution along one spatial direction. Thus, they cannot be used when dealing with both long exposure and high resolution simulations. A multiscale stochastic model has been proposed in [21] to allow the synthesis of long and accurate high resolution turbulent phase sequences with acceptable computational and memory requirements, i.e. O(n). In the approach proposed in [21] the turbulence is simulated at low resolution exploiting dynamic models as in [18–20], whereas the high resolution representation is obtained by means of a multiscale stochastic model [22–26]. In the phase screen simulation literature, similar works based on different multiscale models were also considered in [27,10]. This work considers a generalization of the method for turbulence simulation proposed in [21,28]. First, a PCA local representation has been used in combination with the multiscale stochastic model proposed in [21]. In typical working conditions, this representation allows us to reduce the computational load required by the simulation algorithm of approximately a factor of 2. Furthermore, a global time reduction of approximately a factor of 4 can be obtained by combining the use of the PCA representation with the reduction of the order of the multiscale model at the highest resolution scale. Second, a different low resolution model has been considered with respect to [21]. The Moving Average (MA) model proposed here allows us to evolve the turbulence with a wind velocity vector not aligned with one of the two spatial axes,4 without any significant computational complexity issue. Finally, different from [21] and [28], the simulation of nonstationary turbulence is considered as well. The frozen flow model of the atmosphere is valid at time scales of the order of tens/ hundreds of milliseconds [29]. Thus, for long turbulence simulation, a more general model should be considered. In order to improve the statistical accuracy in long turbulence simulations, this paper proposes a more general non-frozen flow simulation model. Furthermore, it also considers the case of efficiently deal with changes of the value of certain spatial characteristics of the turbulence (i.e. changes of the value of the Fried parameter). The paper is organized as follows: Section 2 introduces the turbulence spatio-temporal statistical characteristics. Section 3 presents the PCA local representation of the turbulent phase. Then, Section 4 introduces the multiscale stochastic model used for the simulation of stationary atmospheric turbulence. Section 5 considers the case of nonstationary turbulence simulation. Finally, in Section 6 the proposed procedure is tested on some simulations.

2. Turbulent phase characterization 2.1. Spatial characteristics From a spatial point of view, the turbulent phase is usually assumed to be zero-mean stationary and spatially homogeneous. It is supposed to be normally distributed [30], hence the second order statistics are sufficient to completely describe its statistical properties. Let u and v be two unit vectors indicating two orthogonal spatial directions, as in Fig. 1, and let ϕðu; v; tÞ be the value of the turbulent phase on the point (u,v) at time t on the telescope aperture plane, where u and v are the coordinates of the 4

In [21] the wind direction was restricted to be along one of the spatial axes.

Fig. 1. An example of phase screen. Coordinate axes on the telescope image domain are plotted in red. The two points, (u,v) and ðu0 ; v0 Þ, are separated by a distance r on the telescope aperture plane. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

point along u and v. Then, the covariance between two values of the turbulence, ϕðu; v; tÞ and ϕðu0 ; v0 ; tÞ, depends only on the distance, r, between the two points: C ϕ ðrÞ ¼ E½ϕðu; v; tÞϕðu0 ; v0 ; tÞ, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 ðu; v; u0 ; v0 Þ, such that r ¼ ðu  u0 Þ2 þ ðv  v0 Þ2 . According to the Von Karman theory, the shape of the spatial covariance function, C ϕ ðÞ, is completely characterized by the values of two physical parameters, r0, the Fried parameter, and L0, the outer scale [1,15]: C ϕ ðrÞ ¼

 5=3     L0 η 2π r 5=6 2π r K 5=6 ; L0 r0 2 L0

ð1Þ

where K  ðÞ is the MacDonald function (modified Bessel function of the third type), and η is a constant [1,15,21]. The outer scale of spatial coherence, L0, is a parameter that characterizes the greatest size of turbulence eddies, whereas the Fried parameter, r0, corresponds to the diameter of a circular area over which the variance of the phase distortion (without any AO correction) is 1 rad2 [1]. Notice that this, 1 rad2, is commonly accepted as the maximum tolerable distortion to have usable astronomic observations. Roughly speaking, the Fried parameter is related to the strength of the turbulence. In particular, when the telescope resolution is limited by the atmospheric turbulence (e.g. for large telescopes without AO correction), then the resolution corresponds to that of a telescope of diameter equal to r0. Notice that, mostly because of the presence of wind and of fluctuations of the atmosphere refraction index, in practical applications the spatial statistical characteristics of the turbulence, and in particular the value of r0, are not constant in long time intervals. This case will be considered in Section 5.

2.2. Temporal characteristics From a temporal point of view, the turbulence is typically modeled as the superposition of a finite number ℓ of independent layers, moving at different altitudes, with different energies and velocities. Let ψ i ðu; v; tÞ be the value of the ith layer at point (u,v) and time t: the ith layer models the atmosphere from hi  1 to hi meter high, where hl Z⋯ Z hi Z hi  1 Z ⋯ Zh0 ¼ 0. Then, ϕðu; v; tÞ

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

3

is given by ℓ

ϕðu; v; tÞ ¼ ∑ γ i ψ i ðu; v; tÞ;

ð2Þ

i¼1

where fγ i g are suitable normalized coefficients, i.e. ∑li ¼ 1 γ 2i ¼ 1. The coefficients fγ 2i g have the physical meaning of strengths (or normalized energies) of the layers. The layers are assumed to be characterized by similar spatial characteristics, i.e. all the layers are spatially described by the same covariance function C ψ ðrÞ: ð3Þ C ψ ðrÞ ¼ E½ψ i ðu; v; tÞψ i ðu0 ; v0 ; tÞ ¼ C ϕ ðrÞ; i ¼ 1; …; ℓ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ ðu  u0 Þ2 þ ðv  v0 Þ2 . A commonly agreed assumption (frozen flow hypothesis [1,29]) considers that each layer translates in front of the telescope pupil with constant velocity vi ¼ vi;u u þ vi;v v:

ψ i ðu; v; t þk=f s Þ ¼ ψ i ðu  vi;u k=f s ; v vi;v k=f s ; tÞ;

ð4Þ

where fs is the sampling rate, and k=f s is a delay multiple of 1=f s . In the simulation of stationary turbulent phase, in Section 4, the frozen flow hypothesis (4) is assumed to hold. In addition, r0 and L0 are fixed to constant values. However, the behavior of the real turbulence actually deviates from that established by the frozen flow approximation (4): in order to take into account of such deviations, the following, more general, model is considered in the simulation of nonstationary turbulence (Section 5):

ψ i ðu; v; t þ1=f s Þ ¼ αψ i ðu  vi;u =f s ; v  vi;v =f s ; tÞ þ βei ðu; v; tÞ;

ð5Þ

where fei gi ¼ 1;…;ℓ have the same spatial statistical characteristics of ϕ and fψ i g. ei is independent of ei0 , 8 i0 a i, and of ψ i ðu; v; t 0 Þ, 8 ðu; vÞ, t 0 rt. α is a scalar factor, 0 o α r 1,pthat determines the temporal ffiffiffiffiffiffiffiffiffiffiffiffiffi decorrelation of the layer, and β ¼ 1  α2 . In the nonstationary case considered in Section 5, r0, is allowed to be time varying, as well.

2.3. Remark Since the layers are independent, each of them can be simulated separately and then combined by means of (2). Hence, without loss of generality, hereafter the problem of simulating a single turbulent layer is considered instead of that of the overall turbulence simulation. Then, the velocity of such layer is vðtÞ ¼ vu ðtÞu þ vv ðtÞv, and its strength γ 21 ¼ 1. The turbulent phase screen (Fig. 1) is assumed to be represented as an r  c matrix containing the turbulent phase values evaluated on a discrete grid, where the physical dimension of each pixel in the matrix is ps  ps . Then, the temporal turbulence evolution is equivalent to generating new rows and columns of the phase screen matrix and properly shifting the window corresponding to the telescope aperture of ðvv ðtÞ=f s ; vu ðtÞ=f s Þ pixels at each sampling period5 (see Fig. 2). Since the problem of simulating the temporal turbulence evolution can be reduced to that of generating a phase screen of larger size, then, to simplify the notation, hereafter we will omit the time coordinate t from equations, and we focus on the goal of generating an r  c phase screen. Then the temporal evolution can be straightly derived from the spatial extension of the current phase screen and (4) (in Section 4) or (5) (in Section 5). 5 For high accuracy simulations, vv =f s and vu =f s are restricted to be integer multiples of ps.

Fig. 2. Temporal evolution of the turbulence is obtained by generating new rows and columns of the phase screen matrix and shifting of v=f s pixels the window corresponding to the telescope aperture.

3. Multiscale spatial representation based on local PCA A multiscale representation of the turbulent phase is obtained by means of a local PCA decomposition. PCA is a statistical tool used to convert a set of possibly correlated random variables into a set of uncorrelated variables by means of a linear transformation. The transformation is defined to successively maximize the variance of the ith component, for i starting from 1 to the number of considered variables [31]. From an algorithmic point of view, the PCA transformation can be computed by applying the Singular Value Decomposition to the covariance matrix Σ of the vector formed by the considered random variables [31,32]. Specifically, Σ ¼ US2 U > , where U contains the spatial bases associated to the PCA decomposition, and S is a diagonal matrix of singular values in descending order. Let the order l of the representation level be larger for scales associated to higher resolution. Let M be the scale associated to the highest resolution, then define xM ðu; vÞ ¼ ϕðups ; vps Þ; 8 ðu; vÞ A Z2 , where, with a slight abuse of notation, u and v are two spatial integer indexes. Conversely, scale 0 is the lowest resolution scale. At each scale the 2D domain is bisected in both directions with respect to the representation at the previous resolution level. The coefficients characterizing the turbulence at the spatial position (u,v) at scale l  1 (e.g. xl  1 ðu; vÞ) are obtained from those at scale l (e.g. xl ð2u; 2vÞ) by means of a local PCA decomposition. Consider, as in Fig. 3(a), the values of xl in a proper 2  2 pixel spatial neighborhood of ð2u; 2vÞ. The PCA representation of such pixels is given by the values of ½xl  1 ðu; vÞ xl  1;1 ðu; vÞ xl  1;2 ðu; vÞ xl  1;3 ðu; vÞ > (see Fig. 3(b)): 2 3 2 3 xl  1 ðu; vÞ xl ð2u; 2vÞ 6 7 6x ðu; vÞ 7 6 xl ð2u; 2v þ 1Þ 7 6 7 6 7 ¼ U l Sl V l ðu; vÞ ¼ U l 6 l  1;1 7 ð6Þ 6 xl ð2u þ1; 2vÞ 7 6 xl  1;2 ðu; vÞ 7; 4 5 4 5 xl  1;3 ðu; vÞ xl ð2u þ 1; 2v þ1Þ where Ul is the 4  4 unitary matrix containing the spatial PCA bases,6 Sl is a diagonal matrix of singular values ðsl;1 ; sl;2 ; sl;3 ; sl;4 Þ in descending order, and V l ðu; vÞ is a unit vector. The four components in the score vector Sl V l ðu; vÞ provide the values of h i> xl  1 ðu; vÞ xl  1;1 ðu; vÞ xl  1;2 ðu; vÞ xl  1;3 ðu; vÞ , where xl  1 is a low pass representation of the turbulent phase, while xl  1;1 , xl  1;2 and xl  1;3 provide the spatial details to obtain xl from xl  1 . 6 Since xM ðu; vÞ ¼ ϕðu; vÞ; 8 ðu; vÞ, and for each l the coefficients representing the turbulence at scale l  1 can be linearly obtained from those at scale l, then the second order statistics of xl can be obtained by means of linear combinations of values from (1).

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

Fig. 3. (a) Grid of turbulent phase values at scale l. (b) PCA representation of the turbulent phase values at scale l corresponding to the shaded area in (a). (c) Grid of turbulent phase values at scale l  1. (a) and (c) The same turbulent phase area, but at different resolutions.

1

modified as follows: 2 3 2 3 xl  1 ðu; vÞ xl ð2u; 2vÞ 6x 6 7 ðu; vÞ 7 6 7 6 xl ð2u; 2v þ 1Þ 7 7 6 7 C C l 6 l  1;1 6 xl  1;2 ðu; vÞ 7: 6 xl ð2u þ 1; 2vÞ 7 4 5 4 5 xl  1;3 ðu; vÞ xl ð2u þ 1; 2v þ 1Þ

energy of the 4th PC

0.8

ð7Þ

0.6

where 2 0.4

6 6 U ð:; 1 : 3Þ Cl ¼ 6 6 l 4

0.2

0 1

2

3

4

5

6

7

scale Fig. 4. Energy s2l;4 of the 4th component in the multiscale PCA representation varying the scale l from 1 to M¼ 7. Turbulence parameters: r 0 ¼ 0:11 m, L0 ¼ 60 m, ps ¼ 8:2 mm.

Then, recursively applying (6), the values of xl  1 are used to obtain the values of the representation at scale l  2 (Fig. 3(c)). Notice that, similar to wavelet transforms [33], for each l, the size of the domain of xl  1 is reduced by a factor of 2 in both directions with respect to that of xl. It is worth to notice that, while fxl ; xl;1 ; xl;2 ; xl;3 g are related to the particular turbulence realization, fU l g and fsl;1 ; sl;2 ; sl;3 ; sl;4 g refer only to the turbulence spatial statistics. Specifically, for an isotropic and homogeneous turbulence fU l g do not depend on (u, b) and on l, i.e. U l ¼ U l0 ¼ U, 8 l. Since Ul is unitary, then (6) can be inverted by using U l> ¼ U l 1 : by recursively using such linear relation, fxl ; xl;1 ; xl;2 ; xl;3 ; …g can be linearly expressed in terms of xM, and vice versa. When l is sufficiently large the last component of the PCA representation (6) has typically much lower energy with respect to the others (i.e. sl;4 b sl;3 r sl;2 r sl;1 ). Applying this decomposition to case studies of turbulence with statistical characteristics of common interest, it is possible to show that typically s2l;4 , the energy of xl  1;3 , is usually much larger at the first scales than at the highest resolution ones, e.g. Fig. 4. This suggests that in such a case xl  1;3 can be neglected without affecting the statistical accuracy of the turbulence representation. Thus, the PCA representation can be exploited to reduce the number of values describing the current turbulent phase and, consequently, the simulation time. In order to take into account of the xl  1;3 components neglected in the PCA representation at certain scales, (6) is

0

3

07 7 7 07 5 0

if xl  1;3 is neglected, whereas Cl ¼Ul otherwise. ^ ðu; vÞ be the representation of ϕðu; vÞ obtained by discardLet ϕ ing the 4th component at level l, 8 l Z l. Then, by recursively ^ can be expressed as follows: applying (6) and (7), ϕ and ϕ    v ; Ml 2 2 l¼0 h¼1     u v 0 ; Ml ; þ ξ0 x0 2M  l 2 M 1

3

ϕðu; vÞ ¼ ∑ ∑ ξhl xl;h



u

Ml

ð8Þ

and    v ; Ml ; 2 2 l¼l |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

M1

ϕ^ ðu; vÞ ¼ ϕðu; vÞ  ∑ ξ3l xl;3



u

M l

ð9Þ

ɛðu;vÞ

where fξl gl;h are proper coefficients7 that can be computed from fC l g and (u,v), and the quantity ɛðu; vÞ has been introduced. ^ ðu; vÞ ¼ 0, and E½ðɛðu; vÞÞ2  r Δϕ Then, E½ɛðu; vÞ ¼ E½ϕðu; vÞ  ϕ immediately follows from the Cauchy–Schwarz inequality [34], where h

Δϕ ¼

M

∑

1

M l l ¼ l þ1 4

M

ss

M

s2l;4 þ ∑

∑

l ¼ l þ1

l0 ¼ l þ 1 l0 a l

2

l l0 0: 2M  l  l

ð10Þ

^ ðu; vÞϕ ^ ðu0 ; v0 Þ be the realized covarLet C^ ϕ ðu  u0 ; v  v0 Þ ¼ E½ϕ iances, i.e. C^ ϕ ðu u0 ; v  v0 Þ ¼ E½ðϕðu; vÞ  ɛðu; vÞÞðϕðu0 ; v0 Þ  ɛðu0 ; v0 ÞÞ:

ð11Þ

By means of the Cauchy–Schwarz inequality, it is also possible to derive the following bound on the difference between the realized covariances C^ ϕ ðu  u0 ; v  v0 Þ and the theoretical ones 7

It is possible to show that jξl j ¼ 1=2M  l ; 8 ðl; hÞ. h

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4.1. Lowest level MA model

6

Bound on covariance function error (%)

5

Let C 0ϕ ðu; vÞ be the spatial covariance of x0 evaluated at spatial distance (u,v). Then x0 is modeled as a spatial MA process, where the process x0 at point (u,v) is obtained as a linear combination of values of the random process ε0 in a spatial neighborhood of (u,v) [35]:

5

4

x0 ðu; vÞ ¼ ∑ θ0 ðku ; kv Þε0 ðu ku ; v  kv Þ;

3

ð13Þ

ku ;kv

2

1

0

0

5

10

15

20

25

Separation [m] Fig. 5. Upper bound on the realized covariance error computed by means of (12) as a function of the spatial separation (r in Eq. (12)). Turbulence parameters: r 0 ¼ 0:11 m, L0 ¼ 60 m, ps ¼ 8:2 mm, M ¼7, l ¼ 5.

C ϕ ðrÞ, where r ¼ jðu  u0 ; v  v0 Þj: jC^ ϕ ðu  u0 ; v  v0 Þ  C ϕ ðrÞj ¼  E½ɛðu; vÞϕðu0 ; v0 Þ  E½ϕðu; vÞɛðu0 ; v0 Þ þE½ɛðu; vÞɛðu0 ; v0 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 ΔϕC ϕ ð0Þ þ Δϕ:

ð12Þ

3.1. Numerical example For instance, let M ¼7, l ¼ 5, and the turbulence parameters r 0 ¼ 0:11 m, L0 ¼ 60 m, ps ¼ 8:2 mm. Then, Fig. 5 shows the bound on the error of the reproduced covariances computed by means of (12), the bound is shown in percent with respect to C ϕ ðrÞ. Since the bound provided by (12) is constant, and C ϕ ðrÞ vanishes as the separation r becomes large, then (12) (normalized with respect to C ϕ ðrÞ) provides a useful bound only for r r r, where r is a proper threshold distance, e.g. r ¼ 25 m in Fig. 5.

4. Multiscale simulation of stationary turbulence The goal of this section is that of providing an efficient method for the simulation of stationary turbulence (e.g. r0, L0, wind velocities and energies are fixed to constant values, the frozen flow hypothesis (4) holds). Since the turbulent phase has Gaussian statistics [30], xM ðu; vÞ ¼ ϕðu; vÞ; 8 ðu; vÞ, and fxl ðu; vÞ; xl;1 ðu; vÞ; xl;2 ðu; vÞ; xl;3 ðu; vÞg can be obtained from xM by means of linear transformations, then they all have Gaussian statistics that can be computed from that of the Von Karman model (1). Then, the synthesis procedure can be summarized as follows:

 Generate x0, the turbulence at the lowest resolution, matching the second order statistics computed from (1).

 For each l, 0 r l o M  1, use a proper stochastic model that takes as input xl and provides as output a realization of xl þ 1 that matches the theoretical second order statistics computed from (1). The following subsections present the low-resolution and the multiscale stochastic models that allow us to successfully perform the above steps. For simplicity of exposition, the border effect is discarded, i.e. at each scale the domain of the process is assumed to be an infinite grid.

where ε0 is a zero-mean Gaussian white-noise process and θ0 ð; Þ are the coefficients of the MA process. Since theffi turbulence covariance (1) vanishes for large values of pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ u2 þ v2 , then C 0ϕ ðu; vÞ vanishes as well. Then, the sum in (13) can be practically thought as the sum of a finite number of terms: the range of values of indexes in the sum (13) is assumed to be ku ¼  δ0 ; …; δ0 , and kv ¼  δ0 ; …; δ0 . The spatial correlations C 0ϕ ðu; vÞ, u ¼  δ0 ; …; δ0 , v ¼  δ0 ; …; δ0 can be expressed as the 2D convolution product of θ0 ð; Þ with itself in reverse order. The computation of proper values of the parameters fθ0 ðÞg to make the covariances of the MA process match the theoretical ones can be solved in the (spatial) frequency domain (by using the fast Fourier Transform) as shown in [28]. Then, the simulation of the turbulence at the lowest resolution level is obtained by generating random numbers fε0 ðu; vÞg and filtering them as in (13). Since the generation of new portions of low resolution turbulence is made by using the local model (13), which requires only that the underlying process ε0 has been generated in the region of interest (and at its boundaries), then it allows us to evolve the turbulence in spatial directions different from those of the axes u and v. 4.2. MAR–MA model for level l According to (6), the model at level l, l ¼ f1; …; Mg, takes xl  1 as input and aims at generating xl  1;1 , xl  1;2 and xl  1;3 such that the reconstructed xl matches the theoretical statistics derived from (1) by means of linear operations. In this paper a MAR–MA stochastic model similar to that of [21] is used. The MAR–MA model for xl  1;1 is as follows: xl  1;1 ðu; vÞ ¼ ∑ al;ku ;kv xl ðu  ku ; v  kv Þ ku ;kv

0

0

0

0

þ ∑ θl ðku ; kv Þεl;1 ðu  ku ; v  kv Þ; k0u ;k0v

8 ðu; vÞ;

ð14Þ

where the coefficients fal;ku ;kv g allow a spatial prediction of the values of xl  1;1 ðu; vÞ, given those of xl  1 in a spatial neighborhood of (u,v) (MAR process); the values of εl;1 are randomly sampled from a zero-mean white-noise Gaussian distribution, and then are 0 0 filtered using the coefficients fθl ðku ; kv Þg (MA process, similar to that of Section 4.1). The range of values of ku and kv is  dl ; …; dl , 0 0 while the range of ku , kv is  δl ; …; δl [36]. The MAR–MA stochastic models for xl  1;2 and xl  1;3 are similar to (14). The reader is referred to [21] for a detailed description of the models. Then, the synthesis procedure is summarized in Table 1. In the procedure of Table 1 Steps 5.A, 5.B, 5.C and 6 are performed only if the xl;3 component is not neglected in the PCA representation. 5. Nonstationary turbulence In practical applications, the frozen flow hypothesis (4) is an approximation of the real behavior of the turbulence phase, that,

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

Table 1 Synthesis procedure at scale l. 1. MAR prediction of xl;1 by using xl as in (14) 2. Generate εl;1 and filter it by means of the MA coefficients fθl;1 g as in (14) 3. MAR prediction of xl;2 by using xl and xl;1 : 3.A. Filter xl with the corresponding coefficients in the MAR model for xl;2 3.B. Filter xl;1 with the corresponding coefficients in the MAR model for xl;2 4. Generate εl;2 and filter it by means of the MA coefficients fθl;2 g 5. MAR prediction of xl;3 by using xl, xl;1 and xl;2 : 5.A. Filter xl with the corresponding coefficients in the MAR model for xl;3 5.B. Filter xl;1 with the corresponding coefficients in the MAR model for xl;3 5.C. Filter xl;2 with the corresponding coefficients in the MAR model for xl;3 6. Generate εl;3 and filter it by means of the MA coefficients fθl;3 g 7. Compute xl  1 by combining xl, xl;1 , xl;2 , xl;3

however, deviates from it. Furthermore, changes in the spatial characteristics of the turbulence (e.g. r0, L0) can occur. Then, this section considers the case of nonstationary turbulence simulation. Specifically, Section 5.1 considers the case of changes in the value of the Fried parameter r0, that is quite common in real applications. In the presentation of Section 5.1 the frozen flow model (4) is used, while Section 5.2 extends the simulation model to the non-frozen flow case (5). 5.1. Variation of the value of the Fried parameter r0 In this subsection the case of changes in the value of r0 is considered. Let the value change from r0 to r 00 . Even after such change the turbulent phase is simulated as described in Sections 4.1 and 4.2, however the coefficients of the MAR–MA models and of the low resolution MA model have to be changed correspondingly. The new values of the model coefficients can be obtained by using the general procedure as described in [21,28]. However, if L0 is fixed to a constant value, then the new values of the coefficients (associated to a change of the value of r0) of the models of Sections 4.1 and 4.2 can be conveniently computed from the previous ones. Consequently to the change from r0 to r 00 , the covariance function becomes C 0ϕ ðrÞ ¼ ðr 0 =r 00 Þ5=3 C ϕ ðrÞ, 8 r Z 0, i.e. it is a scaled version of the previous one. It is simply to prove that the fU l g matrices are invariant with respect to r0. Similarly, the fC l g matrices do not change if the 4th PCA component is neglected for the same values of the scale l. Thus, the covariances of the new process values fxl ; xl;1 ; xl;2 ; xl;3 gl ¼ 0;…;M  1 are scaled version of the previous ones, where, in analogy with the covariance function C 0ϕ ðÞ, the scale factor is ðr 0 =r 00 Þ5=3 . Thus, from the equation of the minimum variance predictor [35], it immediately follows that the coefficients of the MAR predictors are invariant with respect to changes of the Fried parameter. Instead, the coefficients of the MA models are not invariant to changes of the Fried parameter. However, their values can easily be updated from the previous ones. For instance, consider the new 0 values of the MA coefficients fθ0 ðÞg of the low resolution model. Then, it is easy to prove that ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s  ffi r 0 5=3 θ00 ðku ; kv Þ ¼ θ0 ðku ; kv Þ; 8 ðku ; kv Þ: ð15Þ r 00 Finally, when switching from one set of values for the model parameters to another set, it is necessary to ensure the spatiotemporal continuity of the generated turbulence. This can be considered at each scale and for each component separately. Without loss of generality, let us consider the xl  1;1 component. Then, ensuring the spatio-temporal continuity corresponds to satisfy the following equation in the spatial neighborhood region

N where the change occurs x0l  1;1 ðu; vÞ ¼ xl  1;1 ðu; vÞ;

8 ðu; vÞ A N;

ð16Þ

where xl  1;1 ðu; vÞ can be expressed as in (14), while x0l  1;1 ðu; vÞ is as follows: x0l  1;1 ðu; vÞ ¼ ∑ al;ku ;kv xl ðu  ku ; v  kv Þ ku ;kv

0

0

0

0

0

þ ∑ θl ðku ; kv Þε0l;1 ðu ku ; v  kv Þ: 0

ð17Þ

0

ku ;kv

Since the MAR part of the above expression corresponds to that in (14), then the continuity is ensured by imposing 0

∑ θl ðku ; kv Þεl;1 ðu  ku ; v  kv Þ ¼ ∑ θl ðku ; kv Þε0l;1 ðu  ku ; v  kv Þ:

ku ;kv

ð18Þ

ku ;kv

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5=3ffi 0 r 0 =r 00 θ0 ðku ; kv Þ, 8 ðku ; kv Þ, then the soluSince θ0 ðku ; kv Þ ¼ tion of the above equation is simply sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  0 5=3 r0 0 εl;1 ðu; vÞ ¼ εl;1 ðu; vÞ: ð19Þ r0 Similar considerations can be repeated for ε0l;2 and ε0l;3 . 5.2. Non-frozen flow temporal evolution It has already been introduced, in (5), a non-frozen flow model that is here discussed more in detail. In (5), the factor α r 1 determines the temporal decorrelation of the turbulence. Specifically, for α o 1, even a static layer (with zero velocity) has temporal covariances that decrease exponentially, i.e. proportional to αk=f s , where k=f s is a multiple of the sampling period. At each time instant t, turbulent phase simulation by means of (5) is achieved with the following procedure:

 new values of ψi are generated according to the procedure 

described in Section 4 to take into account of the turbulence shift due to the wind velocity, ei ðu; v; tÞ, for each (u,v) in the telescope aperture, is randomly generated, independent of ψi, by using the procedure described in Section 4 (ei have the same statistical characteristics of ψi).

The above procedure allows us to obtain a more accurate simulation of the turbulent phase, however it is unattractive from the computational point of view. Indeed, at each time instant ei has to be generated over all the telescope aperture, while ψi is generated only over a portion of the domain proportional to the wind velocity. The minimum increase of the computational complexity is for v=f s Z d: even in such case the computational load of the simulation increases by a factor of 2 with respect to the frozen flow simulation obtained by means of (4). Notice that even more general models than (5) can be considered for the temporal evolution of the turbulence: for instance, an ARMA (autoregressive moving average) model can be used to obtain a more general temporal spectrum of the turbulent phase, however, the use of such model can take to computational complexity and the memory requirement issues (both of them increase linearly with respect to the order of the considered ARMA model).

6. Numerical simulations 6.1. Computational burden reduction As stated in Section 3, the PCA representation allows us to reduce the overall number of coefficients describing the current value of the turbulent phase (by neglecting the last component of the PCA representation at the highest resolution scales).

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

7

8

structure function [rad2]

5000 4000 3000 2000 1000 0

0

5

10

15

20

25

30

35

40

difference of reproduced normalized structure functions

normalized structure function error

0.015 6000

0.01 0.005 0

−0.005 −0.01 −0.015

0

5

10

separation [m]

15

20

25

30

35

40

separation [m]

x 10−6

6 4 2 0 −2 −4 −6 −8

0

5

10

15

20

25

30

35

40

45

separation [m]

Fig. 6. (a) Comparison of the theoretical structure function (red dotted line) with the sample structure function obtained by means of the complete MAR–MA model (blue solid line), the reduced complexity MAR–MA (black dashed line), the MAR–MA model obtained discarding last component 8 l (green dash–dotted line), estimated from a 40 m  20,000 m phase screen. (b) Error with respect to the theoretical structure function. (c) Difference between the structure functions obtained by the complete and the reduced complexity MAR–MA model. Turbulence parameters: r0 ¼0.11 m, L0 ¼ 60 m. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Consequently, this leads to a reduction of the computational time of the simulation procedure. Since the simulation of both stationary and nonstationary turbulence basically relies on the use of the simulation procedure described in Section 4, then in both cases the PCA representation allows us to obtain a similar computational time reduction: without loss of generality such computational complexity reduction will be estimated for the stationary case in the following. With respect to the procedure summarized in Table 1, a realistic assumption is that Steps 1, 2, 3.A, 3.B, 4, 5.A, 5.B, 5.C and 6 have approximately the same computational cost, while the computational load of Step 7 can usually be neglected. Since the lowest order model of Section 4.1 does not significantly influence the overall computational complexity of the turbulence simulation algorithm, then it will be discarded in the following analysis. Then, discarding the last component of the PCA representation at all the scales allows us to reduce the computational time of 4=9  44%, approximately, since Steps 5.A, 5.B, 5.C and 6 are not necessary when the last component is discarded. As suggested in Section 3, in order to preserve the statistical accuracy of the synthesis method while reducing the simulation computational load, the last component of the PCA representation can be discarded only for scales l Zl, that are the most expensive ones from a computational point of view. For instance, let M¼7, and l ¼ 5, then the computational time reduction is approximately8:

γ

1 ∑M l ¼ 5 4 Ml

4

 19 ∑M l¼0

4

γ

 42%;

ð20Þ

Ml

where γ is a proper constant such that γ =4M  l is the computational time for processing each of the filtering steps at scale l in the turbulent phase generation procedure. Thus, discarding xl;3 for l ¼ fM  2; M  1g allows us to obtain a computational time reduction similar to that obtained discarding xl;3 for all l. On the other hand, the error on the realized covariance function results to be much smaller than that obtained by discarding xl;3 , 8 l. Fig. 5 shows the upper bound (5) on the realized covariance function obtained by setting the turbulence parameters as in the case of Fig. 4. 8

Discarding operations with computational cost typically much lower than filtering the (multiscale) turbulent phase with filters of size dl  dl .

Furthermore, if δ1 ¼ δ2 ¼ ⋯ ¼ δM , then approximately 75% of the total computational time is spent for filtering at the highest resolution. However, if the neighborhood size at scale M is 0 0 reduced, i.e. δM ¼ δM =2, while δl ¼ δl for l o M, then the cost of filtering at scale M is approximately reduced by a factor of 4 (while the cost of operations at the other scales is unchanged). Thus reducing the neighborhood size by a factor of 2 just at scale M allows a computational time reduction on the overall algorithm of, approximately, 34  75%  56%. If this computational complexity reduction is combined with that obtained by means of the PCA representation in the example of Figs. 4 and 5, then an overall  73% time reduction is obtained.

6.2. Accuracy The statistical accuracy of the version of the algorithm proposed above,9 namely the reduced complexity MAR–MA model hereafter, will be compared in the following with the complete model version in two case studies. Astronomers commonly describe the turbulent phase spatial statistical characteristics by means of the structure function, Dϕ ðÞ, that can be defined as follows: Dϕ ðrÞ ¼ 2ðC ϕ ð0Þ  C ϕ ðrÞÞ:

ð21Þ

Figs. 6 and 7 compare the theoretical turbulence statistics with the sample ones estimated from 40 m  20,000 m phase screens generated by means of: the complete MAR–MA model (the last component of the PCA representation is never discarded), the reduced complexity MAR–MA model previously suggested (the last component of the PCA representation is discarded at the two highest resolution scales, and at the highest resolution a reduced spatial neighborhood is considered in the MAR–MA), and the reduced complexity MAR–MA model where the last component of the PCA representation is discarded at all the scales. Specifically, Figs. 6(a) and 7(a) compare the theoretical and the sample structure functions, whereas Figs. 6(b) and 7(b) compare the error on the sample structure functions with respect to the theoretical ones. Finally, Figs. 6(c) and 7(c) show the difference 9 The last component of the PCA representation, xl;3 , is discarded for scales l Z M  2, and neighborhood size at the highest resolution scale is reduced, i.e. δ0M ¼ δM =2, while δ0l ¼ δl for l o M.

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

1.5

2000

1500

1000

500

0

difference of reproduced normalized structure functions

normalized structure function error

structure function [rad2]

0.015 0.01 0.005 0

−0.005 −0.01

5

10

15

20

25

30

1 0.5 0 −0.5 −1 −1.5

−0.015

0

x 10−5

0

5

10

separation [m]

15

20

25

separation [m]

30

0

5

10

15

20

25

30

separation [m]

Fig. 7. (a) Comparison of the theoretical structure function (red dotted line) with the sample structure function obtained by means of the complete MAR–MA model (blue solid line), the reduced complexity MAR–MA (black dashed line), the MAR–MA model obtained discarding last component 8 l (green dashed-dot line), estimated from a 40 m  20,000 m phase screen. (b) Error with respect to the theoretical structure function. (c) Difference between the structure functions obtained by the complete and the reduced complexity MAR–MA model. Turbulence parameters: r0 ¼ 0.15 m, L0 ¼ 45 m. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

between the structure functions obtained by the complete and the reduced complexity MAR–MA model. The simulation parameters are set to typical values of interest for extremely large telescopes: telescope diameter is set to d ¼40 m, L0 ¼ 60 m and r 0 ¼ 0:11 m in Fig. 6, while L0 ¼ 45 m and r 0 ¼ 0:15 m in Fig. 7. The pixel size is set to 0.01 m  0.01 mm, and the total number of scales is set to M¼6.

6.3. Discussion As expected, the reduced complexity (with the last component of the PCA representation discarded at the two highest resolution scales) and the full model provide practically equivalent results (Figs. 6(c) and 7(c)): the sample structure function estimated by synthesized data is very close to the theoretical one, while the computational time of the reduced complexity model is approximately 1/4 of that of the complete model: thus, the proposed model allows a particularly relevant time reduction without compromising the statistical accuracy. On the other hand, discarding the last PCA component at all the scales leads to a nonsignificant time reduction with respect to (20), and, as shown in Figs. 6(b) and 7(b), it reduces the statistical accuracy of the model. Hence, a good tradeoff between computational time reduction and statistical accuracy in turbulence realization can be obtained by describing the turbulent phase with all its principal components at the coarse resolution scales, while discarding the last PCA component only at the scales with highest resolution, that are the most expensive ones from the computational point of view (i.e. xl;3 is discarded for l Zl, with l r M  1). Furthermore, the synthesis procedure can take advantage of the low resolution model of Section 4.1 to dynamically evolve the turbulence in whatever direction: while the low level model in [21] accurately reproduces the turbulence statistics only if the turbulent phase moves along either the u or the v direction, the local model proposed here removes such requirement (actually, to ensure a good statistical accuracy, it only requires that vv =f s and vu =f s are integer multiples of ps). To conclude, two observations are now in order for what concerns the case of nonstationary turbulent phase simulation: first, when the outer scale L0 is fixed to a constant value, the proposed simulation procedure deals very efficiently with changes in the value of the Fried parameter r0. Second, the non-frozen flow model (5) has been proposed to generalize the frozen flow turbulence model (4). This allows a more realistic simulation of

the turbulence, however, at the price of an increase in the computational complexity of the simulation.

7. Conclusions This paper considers the problem of simulating high resolution phase screens. The contribution of this work is the generalization of the method proposed in [21,28] in several directions. First, thanks to the use of a multiscale local PCA representation, it significantly reduces the computational load of the method presented in [21]. Furthermore, the low resolution MA model proposed in this paper allows to dynamically evolve the turbulent phase indifferently in all the spatial directions, with computational cost and memory requirements similar to those in [21]. Finally, this paper considers the problem of non-stationary turbulence simulation: the proposed method deals efficiently with changes of the Fried parameter, and considers the non-frozen flow turbulence simulation as well.

References [1] Roddier F. Adaptive optics in astronomy. Cambridge, UK: Cambridge University Press; 1999. [2] Wiberg D, Max C, Gavel D. Geometric view of adaptive optics control. J Opt Soc Am A 2005;22(5):870–80. [3] Gendron E, Léna P. Astronomical adaptive optics. 1: modal control optimization. Astron Astrophys 1994;291:337–47. [4] Maneri E, Gawronski W. LQG controller design using GUI: application to antennas and radio-telescopes. ISA Trans 2000;39(2):243–64. [5] Gawronski W, Ahlstrom H, Bernardo A. Analysis and performance of the control systems of the NASA 70-meter antennas. ISA Trans 2004;43(4):597–610. [6] Hardy J. Active optics: a new technology for the control of light. Proc IEEE 1978;66(6):651–97. [7] Schipani P, D0 Orsi S, Fierro D, Marty L. Active optics control of VST telescope secondary mirror. Appl Optics 2010;49(16):3199–207. [8] Poyneer L, Véran J. Predictive wavefront control for adaptive optics with arbitrary control loop delays. J Opt Soc Am A 2008;25(7):1486–96. [9] Correia C, Raynaud H, Kulcsár C, Conan J. On the optimal reconstruction and control of adaptive optical systems with mirror dynamics. J Opt Soc Am A 2010;27(2):333–49. [10] Thiébaut É, Tallon M. Fast minimum variance wavefront reconstruction for extremely large telescope. J Opt Soc Am A 2010;27(5):1046–59. [11] Béchet C, Tallon M, Tallon-Bosc I, Thiébaut É, Le Louarn M, Clare R. Optimal reconstruction for closed-loop ground-layer adaptive optics with elongated spots. J Opt Soc Am A 2010;27(11):A1–8. [12] Rosensteiner M. Wavefront reconstruction for extremely large telescopes via CuRe with domain decomposition. J Opt Soc Am A 2012;29(11):2328–36.

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i

A. Beghi et al. / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎ [13] Shatokhina I, Obereder A, Rosensteiner M, Ramlau R. Preprocessed cumulative reconstructor with domain decomposition: a fast wavefront reconstruction method for pyramid wavefront sensor. Appl Optics 2013;52(12):2640–52. [14] de Visser C, Verhaegen M. Wavefront reconstruction in adaptive optics systems using nonlinear multivariate splines. J Opt Soc Am A 2013;30(1):82–95. [15] Conan R. Modelisation des effets de l0 echelle externe de coherence spatiale du front d0 onde pour l0 observation a haute resolution angulaire en astronomie [Ph.D. thesis]. Univ. Nice, France: Nice Sophia Antipolis; 2000. [16] Luck R, Agba E. On the design of piezoelectric sensors and actuators. ISA Trans 1998;37(1):65–72. [17] McGlamery B. Computer simulation studies of compensation of turbulence degraded images. In: Proceedings of SPIE, vol. 74, 1976, p. 225–33. [18] Assemat F, Wilson R, Gendron E. Method for simulating infinitely long and non stationary phase screens with optimized memory storage. Optic Express 2006;14(3):988–99. [19] Beghi A, Cenedese A, Masiero A. Stochastic realization approach to the efficient simulation of phase screens. J Opt Soc Am A 2008;25(2):515–25. [20] Fried D, Clark T. Extruding Kolmogorov-type phase screen ribbons. J Opt Soc Am A 2008;25(2):463–8. [21] Beghi A, Cenedese A, Masiero A. Multiscale stochastic approach for phase screens synthesis. Appl Optics 2011;50(21):4124–33. [22] Luettgen M, Karl W, Willsky A, Tenney R. Multiscale representations of Markov random fields. IEEE Trans Signal Process 1993;41(12):3377–96, http: //dx.doi.org/10.1109/78.258081. [23] Benveniste A, Nikoukhah R, Willsky A. Multiscale system theory. IEEE Trans Circuits Syst I: Fundam Theory Appl 1994;41(1):2–15, http://dx.doi.org/ 10.1109/81.260214.

9

[24] Daoudi K, Frakt A, Willsky A. Multiscale autoregressive models and wavelets. IEEE Trans Inf Theory 1999;45(3):828–45, http://dx.doi.org/10.1109/18.761321. [25] Irving W, Willsky A. A canonical correlations approach to multiscale stochastic realization. IEEE Trans Autom Control 2001;46(10):1514–28, http://dx.doi.org/ 10.1109/9.956048. [26] Frakt A, Willsky A. Computationally efficient stochastic realization for internal multiscale autoregressive models. Multidimens Syst Signal Process 2001;12 (2):109–42, http://dx.doi.org/10.1023/A:1011184728562. [27] Lane R, Glindemann A, Dainty J. Simulation of a Kolmogorov phase screen. Waves Random Media 1992;2(3):209–24. [28] Beghi A, Cenedese A, Masiero A. Multiscale phase screens synthesis based on local PCA. In: Proceedings of the 10th IEEE International Conference on Control and Automation. Hangzhou, China: IEEE ICCA, 2013; 2013. p. 1862–7. [29] Schöck M, Spillar E. Method for a quantitative investigation of the frozen flow hypothesis. Appl Optics 2000;17(9):1650–8. [30] Roddier F. The effects of atmospheric turbulence in optical astronomy. Progr Optics 1981;19:281–376. [31] Jackson J. A user0 s guide to principal components. New York, USA: John Wiley & Sons Inc.; 1991. [32] Golub G, Loan CV. Matrix computations. Baltimore, MD: Johns Hopkins University Press; 1989. [33] Mallat S. A wavelet tour of signal processing. San Diego, CA, USA: Academic Press; 1999. [34] Shiryaev A. Probability. Springer-Verlag New York, NJ, USA: Springer; 1995. [35] Söderström T. Discrete-time stochastic systems. Springer-Verlag London, UK: Springer; 2002. [36] Sadabadi M, Shafiee M, Karrari M. Two-dimensional ARMA model order determination. ISA Trans 2009;48(3):247–53.

Please cite this article as: Beghi A, et al. Nonstationary multiscale turbulence simulation based on local PCA. ISA Transactions (2014), http://dx.doi.org/10.1016/j.isatra.2013.12.026i