An algorithm to simulate nonstationary and non-Gaussian stochastic processes

We proposed a new iterative power and amplitude correction (IPAC) algorithm to simulate nonstationary and non-Gaussian processes. The proposed algorithm is rooted in the concept of defining the stochastic processes in the transform domain, which is elaborated and extend. The algorithm extends the iterative amplitude adjusted Fourier transform algorithm for generating surrogate and the spectral correction algorithm for simulating stationary non-Gaussian process. The IPAC algorithm can be used with different popular transforms, such as the Fourier transform, S-transform, and continuous wavelet transforms. The targets for the simulation are the marginal probability distribution function of the process and the power spectral density function of the process that is defined based on the variables in the transform domain for the adopted transform. The algorithm is versatile and efficient. Its application is illustrated using several numerical examples.


Introduction
Observed time histories of the seismic ground motions [31], wind velocity [48], wave height [32], etc. fluctuate randomly in time and space. The time histories are used as the input to carry out the structural analysis. Since the available recorded time histories of the random phenomena are limited, their synthetics are generated and used in practice. The simulation is based on the theory of stochastic processes [7,22,25]. The simulated ground motions are used to assess seismic risk of infrastructures such as dams and portfolio of buildings (e.g., [3,24]). The simulated winds are employed to assess the performance structures and infrastructure system [55,60].
For stationary Gaussian processes, and evolutionary processes [38], the simulation can be carried out using the spectral representation method [23,47], developed based on the ordinary Fourier transform (FT). A stationary process is defined by its power spectral density (PSD) function, and an evolutionary process is defined by the evolutionary PSD that is a function of an amplitude modulation function. The evolutionary process with timedependent amplitude modulation is widely used in generating seismic ground motions [1,10] and fluctuating wind velocity for high-intensity wind events [5,17,20,21].
Masters and Gurley [27] proposed an iterative spectral correction algorithm to simulate the stationary non-Gaussian processes, where the spectral representation method is used in each iteration to generate the time history. They showed that their algorithm outperforms the SRM-based simulation techniques presented in Yamazaki and Shinozuka [58], Gurley and Kareem [15], Grigoriu [14], and Deodatis and Micaletti [11]. It is noted that an algorithm similar to the spectral correction algorithm, namely the iterative amplitude adjusted Fourier transform (IAAFT) algorithm, was proposed by Schreiber and Schmitz [43,44] in the context of generating surrogate for statistical hypothesis testing. The use of the translation process for the stationary non-Gaussian process proposed in Grigoriu [14] was extended for the nonstationary processes by others, including Ferrante et al. [12] and Shields et al. [46].
The evolutionary PSD is often assessed using (timedependent) windowed Fourier transform, such as the short-time Fourier transform [6]. The resolution of such a transform is controlled by the width of the window. As the width of the window increases, a better resolution is obtained at the low frequencies, and the resolution in time deteriorates. A good resolution in both time and frequency (i.e., scale) can be obtained by applying the continuous wavelet transforms (WT) [9,35]. A procedure to estimate the evolutionary PSD by applying the continuous WT was proposed by Spanos and Failla [50]. However, an algorithm that directly applies the continuous WT to simulate the nonstationary stochastic processes with a prescribed wavelet spectrum or time-scale PSD was unavailable. Recently, an iterative algorithm was presented by Chavez and Cazelles [4] to generate surrogate for statistical hypothesis testing. We will point out, in the following sections, a potential weakness of the algorithm, as well as the link between this algorithm and an interesting way of defining nonstationary processes in the wavelet domain introduced by Maraun et al. [26]. The lack of continuous WT-based algorithm to simulate time histories is partly due to that the use of continuous WT does not lead to the decomposed signal to be represented by a set of orthogonal basis functions. Rather than using the continuous WT, the application of the discrete WT and wavelet packet transform that have orthogonal basis functions is presented in Gurley and Kareem [16], and Yamamoto and Baker [57]. The resolution obtained by using these discrete transforms is less refined than that obtained by using the continuous WT.
The phase information in WT is local, while the phase information in the Fourier transform refers to the harmonics at zero time [52]. Stockwell et al. [53] (see also [37]) developed the S-transform (ST) that provides the time-frequency representation of the analyzed signal. It is a hybrid of continuous WT and windowed FT. The Stransform provides frequency-dependent resolution. Similar to the continuous WT, ST does not lead to a decomposed signal to be represented by a set of orthogonal basis functions. Stockwell [52] proposed a discrete orthonormal S-transform. The simulation of the seismic ground motions by using the discrete orthonormal Stransform or combination with ST was presented in Cui and Hong [8] and Hong and Cui [19]. However, an algorithm by using ST alone to simulate nonstationary stochastic processes has not been reported in the literature.
There are other techniques used to simulate the nonstationary processes. These include the application of autoregressive moving-average (ARMA) [42], Karhunen-Loéve expansion [36,49,51], and polynomial chaos [41], and Hilbert-Huang transform [56]. The ARMA uses the recursive relations that connect the random quantity to be simulated at successive time increments. The application of the Karhunen-Loéve expansion uses the eigenfunctions that are obtained by solving Fredholm integral equation of the second kind. In the Hilbert-Huang transform, a random process is decomposed into intrinsic mode functions by the method of empirical mode decomposition. A review of these simulation procedures is beyond the present study since these techniques involve varieties of mathematical concepts and algorithms.
In the present study, we exam and extend the definition of the nonstationary processes in the transform domain. We proposed an iterative power and amplitude correction (IPAC) algorithm to simulate nonstationary and non-Gaussian processes. The algorithm could be viewed as an extension of IAAFT [43] and the spectral correction algorithm [27] and is rooted in the concept of defining the stochastic processes in the transform domain. In particular, we provide details of using the proposed algorithm with FT, ST and WT, where the energy distribution in the transform domain that satisfies energy preservation is prescribed, and the marginal probability distribution function of the process is given. We provide numerical examples to show the proposed algorithm and compare the simulated time histories obtained by using the ST-based and (continuous) WT-based approach.

Fourier transform, S-transform, and wavelet transforms
This section summarizes some basic properties of FT [6,30], ST [37,53], and continuous WT [9]. Only the continuous WT, including its discretized form (which differs from the discrete wavelet transform), is used in the present study. Unless otherwise indicated, WT refers to the continuous WT and its discretized form in the following. The summary provides the basis for the proposed iterative simulation algorithm to be described in the next sections.
Let x(t) denote a realization of a stochastic process such as the ground motion record, X(t). FT of x(t), and its inverse (IFT) can be expressed aŝ and, where FT(•) and IFT(•) denote the FT and IFT operations, the subscript associated with these operators indicates the domain or the index where the operation is carried out;xð f Þ denotes FT of x(t); f is the frequency in Hz,xð f Þ ¼x Ã ð−f Þ, and * denotes the complex conjugate. A symbol or function with a circumflex is used to represent its FT throughout the present study. If x(t) is given in the discrete form x(jΔ t ), j = 0, …, N − 1, with a sampling time interval Δ t and the duration T, T = NΔ t , the (discretized) FT pair is given by, and, where Δ f = 1/T, and the operators FT(•) and IFT(•) that are used for continuous FT are used for discrete FT as well. It is considered implicitly in the following that the numerical calculations ofxðpΔ f Þ and x(jΔ t ) are to be carried by using the fast Fourier transform (FFT) [30] for computational efficiency. Moreover, the notation {•} N is used for the collection of its argument of length N. For example, {x(jΔ t )} N represents all x(jΔ t ) for j = 0, …, N − 1.
ST of x(t) is defined as [37,53], where x S ð f ; τÞ is the ST coefficient, ST(·) denotes the Stransform of its argument, and τ is the center of the window function w(f, τ − t) defined as, The parameter κ in Eq. (6) controls the effective width of the window in ST. It can be shown [53] that, and, where IST(·) is the inverse S-transform (IST). Using Eqs. (7) and (8), the discretized version of x(t) and x S ð f ; τÞ, represented by x(jΔ t ) and x S ðqΔ f ; pΔ t Þ pair, can be written as, and, indicating that the evaluation of the ST coefficients at (pΔ f , qΔ t ) and its inverse at jΔ t is based on FT.
WT is defined as [9,35], where x W ðs; τÞ is the wavelet coefficient, the operator WT(•) denotes WT, ψ(⋅) is the mother wavelet and, s is the scaling or dilation factor, and τ is the translation or position parameter. Eq. (11) can be expressed as [9,35], to facilitate its computation by using FFT for signals given in the discretized form. If the admissibility condi- can be obtained using the following inverse WT [9], where IWT(•) is the inverse of WT(•). If ψ(t) = ψ * (−t), Eq. (13) becomes, Moreover, if the analytical waveletcomplex-valued wavelet function that its FT is null for negative frequencyis used, Eq. (13) can be expressed in Morlet formulation [45], Eqs. (14), (15) and (12) can be written in the following discretized form, and, where c 0 and s 0 are parameters for the numerical computation; K is the total number of scales considered for the numerical integration; T L = − L l Δ t and T U = L U Δ t define the lower and upper limits for the integral over time In the following, we restrict ourselves to the real-valued signal and the analytical wavelets or wavelets with ψ(t) = ψ * (−t).
Gaussian process, power spectral density, and defining process in the transform domain According to the spectral representation method [47] with the use of FFT [59], a sample of a Gaussian stationary process, x(t), can be simulated by transforming Gaussian white noise w(t) to the Fourier domain, multiplying it with an intensity function jŷð f Þj, and transforming it back to the time domain. That is, where e iθ F ðwðtÞÞ ¼ ηðFT ðwðtÞÞÞ , in which the function Since, by definition, the double-sided PSD function of the process x(t) with duration T, S F x ð f Þ, is given by, it indicates that given the target PSD function S F x ð f Þ , one could define a stationary Gaussian process in the Fourier domain by assigning The samples of the process so defined can be obtained using, and the expected PSD of the sampled signals equals the prescribed S F x ð f Þ. The use of the definition given in Eq. (20) preserves the energy of x(t) according to Parseval's theory. We note that by assigning jŷð f Þj equal to MðtÞjxð f Þj, Eq. (19) becomes, which simulates a uniformly amplitude modulated evolutionary process [38]. Such a process has an evolutionary PSD function equals jMðtÞj 2 ðxð f Þx Ã ð f Þ=T Þ, and M(t) is the amplitude modulation function, which will be considered to be positive. However, the use of |y(f)| equal to Mðt; f Þjxð f Þj in Eq. (19) does not lend itself to be interpreted as a proper inverse Fourier transform because the modulation function depends on the frequency. This reduces the computational efficiency that otherwise can be gained by using FFT; it also makes the distinction between the modulation function and intensity function more blurred. We will concentrate only on the case where the modulation function is defined outside of the transform domain. However, the consideration of modulation that depends on variables in the transform domain could be a valid assumption. Maraun et al. [26] emphasized the usefulness of using Eq. (19) to obtain samples of stationary Gaussian process, and extended it to define a class of nonstationary Gaussian processes in the wavelet domain by the wavelet multipliers jy W ðs; τÞj, indicating that an individual process is defined by its multipliers and a synthesizing wavelet pair. Samples of x(t) based on such a definition are then given as, where e iθ W ðwðtÞÞ ¼ ηðWT ðwðtÞÞÞ . We use the intensity function jy W ðs; τÞj and e iθ W ðwðtÞÞ in Eq. (23) instead of using y W ðs; τÞ and WT(w(t)) as suggested in Maraun et al. [26]. The use of e iθ W ðwðtÞÞ instead of WT(w(t)) is aimed at not biasing the energy arising from the intensity function since [WT(w(t))][WT(w(t))] * is not a constant in the wavelet domain by using WT defined in Eq. (12). The use of jy W ðs; τÞj (as well as jŷð f Þj in Eqs. (21) and (22)) is more restrictive than y W ðs; τÞ but is adequate for the proposed algorithm in the following section since we are focused on real-valued signals. However, a negatively valued intensity and complexvalued intensity may be considered for other applications.
Similar to the use of M(t) in defining the uniformly modulated evolutionary process mentioned earlier, we include M(t) in Eq. (23), to define a modulated and intensity function adjusted (MODIF) process. The intensity function gives timescale characteristics of the process, and the modulation function provides additional time-varying characteristics of the process. We further extend the concept of defining the MODIF process in the time-frequency domain according to ST, denoted as the S-domain, where samples of x(t) are given as, where y S ð f ; τÞ is an intensity function in the S-domain, and e iθ S ðwðtÞÞ ¼ ηðST ðwðtÞÞÞ. It is noted that besides the above-mentioned transforms, there are other transforms used for signal analysis and modeling; for example, the generalized Fourier family transforms [2]. Therefore, it is relevant and straightforward to conceptually generalize the approach in defining the MODIF processes in the transform domain if other transform pair is considered. The definitions lend themselves to an easily understandable and almost trivial algorithm to simulate stochastic processes: A) Sample Gaussian white noise, w(t), and calculate the normalized coefficients of w(t) in the transform domain (e.g., e iθ F ðwðtÞÞ , or e iθ W ðwðtÞÞ , or e iθ S ðwðtÞÞ if FT, or WT, or ST is used, respectively). B) Apply the inverse transform to the product of the intensity function and the normalized coefficients obtained in Step A). C) Apply the modulation function to the simulated signal from Step B).
Step C) is separated from Steps A) and B) and is not affected by the selected transformation. A critical issue of applying the MODIF process with prescribed target energy distribution is that the energy distribution of the sampled signals for given intensity function may not be readily established, except for the case where FT is used (i.e., transforms with non-redundant representation). This is because unlike the FT, both WT and ST provide redundant representation. The redundant representation results in that, in general, jy W ðs; τÞje iθ W ðwðtÞÞ and jy S ðs; τÞje iθ S ðwðtÞÞ do not represent the proper coefficients of WT and ST, respectively. In other words, jy W ðs; τÞje iθ W ðwðtÞÞ and jy S ðs; τÞje iθ S ðwðtÞÞ are not equal to x W ðs; τÞ ¼ WT ðIWT ðjy W ðs; τÞje iθ W ðwðtÞÞ ÞÞ and x S ð f ; τÞ ¼ ST ðIST ðjy S ðs; τÞje iθ S ðwðtÞÞ ÞÞ, respectively.
To see the impact of this inequality on the simulated MODIF process by using Eq. (24), we note that we can define the double-sided time-scale PSD (TSPSD) function of the simulated process x(t), S Wx ðs; τÞ, as, The use of this definition leads to energy preservation since the integral of S Wx ðs; τÞ in the wavelet domain equals the integral of |x(t)| 2 in the time domain (see Proposition 2.4.1 in Daubechies [9]). Consequently, even we assign jy W ðs; τÞj equals ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi S Wx ðs; τÞC ψ p jsj and M(t) = 1 for the simulation, the average energy of the sampled signals according to Eq. (24) will likely deviate from the specified target S Wx ðs; τÞ.
Consider that we simulate the MODIF process using Eq. (25). We can define the double-sided time-frequency PSD (TFPSD) function of the simulated process, S Sx ð f ; τ Þ, as, since the use of this definition leads to energy preservation [18] Â expð−ð2πκðζ−1ÞÞ 2 Þdζ . However, the average energy of the sampled signals by using Eq. (25) with jy S ð f ; τÞj equal to jx S ð f ; τÞj ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi S Sx ð f ; τÞD κ j f j p and M(t) = 1 will likely deviate from the specified target.
In addition to the discussed energy distortion, the application of the MODIF process is likely to lead to the samples obtained from Eqs. (22), (24) and (25) to follow a marginal cumulative distribution function (CDF) that deviates from the prescribed marginal CDF of the zeromean process F X, t (x(t)). An iterative process is proposed in the following sections to simulate the nonstationary and non-Gaussian with prescribed target PSD and CDF. The PSD functions that satisfy the energy preservation by considering the selected transform are used as the basis to describe the proposed algorithm to maintain consistency. Although this could become clumsy in some instances, it is useful in checking that a consistent transform pair is employed.

IAAFT algorithm
To develop the proposed iterative algorithm, we note that, given the observed {x(jΔ t )} N , the IAAFT algorithm was proposed by Schreiber and Schmitz [43,44] in the context of generating surrogates for statistical hypothesis testing. The algorithm repeatedly uses FT and IFT, and ranked data. This algorithm is explained using the ranking of x(jΔ t ) in the following. The PSD function S F x ð f Þ of {x(jΔ t )} N is calculated using Eq. (19) with possible smoothing. The objective of IAFFT is to generate surrogates that match the calculated S F x ð f Þ and shuffled {x(jΔ t )} N . A similar algorithm -the spectral correction algorithm -was independently designed by Masters and Gurley [27] to simulate non-Gaussian processes for the given target S F x ð f Þ and target marginal CDF F X (x(t)). A subtle difference between these two algorithms is how the prescribed target PSD function and CDF are obtained or assigned. For example, {x(jΔ t )} N is obtained through distribution mapping in the spectral correction algorithm. In IAAFT, {x(jΔ t )} N is given and shuffled. This shuffling, in the spectral correction method, can be viewed as matching the prescribed probability distribution. Once {x(jΔ t )} N is prescribed and S F x ð f Þ is calculated, by letting {ξ(j)} N equal to the ascendingly sorted {x(jΔ t )} N , the steps of the IAAFT algorithm are: 1. Sample a sequence of Gaussian white noise, w(t), of length N, calculate e iϕ p ¼ ηðFT ðwðtÞÞÞ; je iϕ p Þ and find the ranking of x PC (jΔ t ), r j , for j = 0, …, N − 1, based on the ascending order; 3. Set x AC (jΔ t ) = ξ(r j ), for j = 0, …, N − 1; and calculate e iϕ p ¼ ηðFT j ðx AC ð jΔ t ÞÞÞ, and. 4. Repeat Steps 2) to 3) until the convergence criterion is achieved.
Steps 1) and 2) are the same as Steps A) and B) described earlier that simulates a Gaussian process, except an additional ranking of x PC (jΔ t ) is carried out, which is equivalent to define the CDF as a preparation for the iteration. In general, Step 2) leads to x PC (jΔ t ) with the PSD correction but may deviate from the target CDF assigned by {ξ(j)} N , and Step 3) leads to the sampled x AC (jΔ t ) with the amplitude correction (i.e., matching CDF assigned based on {ξ(j)} N ) but may deviate from the target PSD. The iteration adjusts the PSD and CDF of the sampled time series to their corresponding targets. The tolerable differences between x PC (jΔ t ) and x AC (jΔ t ) can be used as the convergence criterion. Once convergence is achieved x PC (jΔ t ) or x AC (jΔ t ) can be used as the sampled time series. Note that since the time increment equals Δ t , the corresponding Nyquist frequency for the considered signal equals 1/(2Δ t ).
The IAAFT algorithm is designed for stationary processes. For the shuffling of {x(jΔ t )} N to simulation stationary process, it is implicitly considered that the marginal CDF of x(t) at any given time remains to be the same. Also, the PSD function for the stationary process is time-independent. The IAAFT algorithm or the spectral correction method is not applicable to simulate nonstationary processes as they have time-varying PSD and CDF.

Proposed iterative power and amplitude correction algorithm
In this section, we describe the proposed iterative power and amplitude correction (IPAC) algorithm to simulate the time history {x(jΔ t )} N of a zero-mean nonstationary non-Gaussian process. Since Δ t is assigned, N can be determined based on the length of the signal to be simulated and vice versa, and the Nyquist frequency for the sampled signal equals 1/(2Δ t ). The proposed algorithm could be viewed as an extension to the IAAFT algorithm. For the simulation, it is considered that, for M(t) = 1, the PSD function of the process that is characterized based on FT, or ST, or WT is given, and the distribution type for the marginal CDF of x(t), F X, t (x(t)), is known. Moreover, it is considered that F X, t (x(t)) can be completely defined by the zero-mean, the time-varying standard deviation, σ(t), and other prescribed distribution parameters if they are required (since, in some cases, a CDF with more than two parameters may be considered).
If FT is considered for a stationary process, the standard deviation σ(t) equals ffiffiffiffiffiffiffiffi λ F x p which is timeindependent, where λ F x , equals the integral of S F x ð f Þ over the frequency domain. Since S Sx ð f ; τÞ provides the energy distribution over the time-frequency domain, the integral of S Sx ð f ; τÞ over the frequency domain provides the energy distribution in the time domain, λ Sx ðτÞ, and the integral of S Sx ð f ; τÞ over the time domain provides the energy distribution in the frequency domain, S Sx ð f Þ . Analogously to the statistics for the stationary process, λ Sx ðτÞ represents the variance of x(τ), and σ(t) equals ffiffiffiffiffiffiffiffiffiffiffiffi ffi λ Sx ðtÞ p . Similarly, for the given S Wx ðs; τÞ , the time-varying variance λ Wx ðτÞ is given by, and σ(t) equals ffiffiffiffiffiffiffiffiffiffiffiffiffiffi λ Wx ðtÞ p . The integral of S Wx ðs; τÞ over the time domain provides the energy distribution in the scale domain, S Wx ðsÞ.
Let u(t) be a uniformly distributed random variable between 0 and 1 with its marginal CDF denoted as U(u(t)). The relation between u(t) and x(t) can be established based on the probability transformation, U(u(t)) = F X, t (x(t)). The steps in the IPAC algorithm in a pseudocode form are shown in the flowchart depicted in Fig. 1 and are described as follows:

I) Prescribe the targets and initiate the simulation process:
Sample {u(jΔ t )} N based on a random number generation algorithm for a uniformly distributed random variable between 0 and 1. Assign {p(j)} N equal to the ascendingly sorted {u(jΔ t )} N , and the intensity function j y T ð s jΔ t ÞÞ , and find the rank of p PC (jΔ t ), denoted as r j , for j = 0, …, N − 1; and calculate e iϕ T ð s  Since x AC (jΔ t ) is used in Step II.5), the distribution match (i.e., matching samples of F X (X(t)), {x(jΔ t )} N ) is ensured by design. One could replace Step II.5) with x(jΔ t ) = M(jΔ t ) × x PC (jΔ t ) without altering the results if a stringent convergence criterion is employed. As can be observed from the steps of the IPAC algorithm, the analysis, as well as the simulation, is carried out within the same transform pair. It avoids the need to map the obtained results from one type of transform into a different kind of transform (e.g., obtaining the spectrum using continuous WT and then transform it into evolutionary PSD). The algorithm emphasizes the ranking and matching of the probability values. Note that it may be attempting to replace the uniform distribution with the normal distribution for u(t). However, by doing so, it requires the use of the inverse distribution transformation in Steps II.2 and II.3) and increases computing demand.
The algorithm can be simplified if F X, t (x(t)) remains unchanged and only depends on σ(t), that is, the marginal probability distribution of z(t) = x(t)/σ(t), F Z (z(t)), is timeindependent and z(t) has zero mean and unit variance. In such a case, we calculate fζð and " x AC (jΔ t ) = ζ(r j )σ(jΔ t ) ", respectively. This avoids the use of probability distribution function during the iteration to gain extra computational efficiency. This simplified version can also be used to generate surrogate for observed {x(jΔ t )} N , which has the effect of the modulation function already removed. This is done by calculating {z(jΔ t )} N = {x(jΔ t )/σ(jΔ t )} N , and letting {ζ(j)} N equal to the ascendingly sorted {z(jΔ t )} N in Step I.1) (instead of fζð jÞg N ¼ fF −1 Z ðpð jÞÞg N ), where σ(jΔ t ) is to be calculated based on the PSD function estimated from {x(jΔ t )} N by using a preferred transform.
The usefulness of surrogate in the context of wind engineering was presented in McCullough and Kareem [28]. The proposed algorithm, when used with WT to generate surrogate, differs from that given in Chavez and Cazelles [4] for testing time-localized coherence, in that the time-varying σ(jΔ t ) is neglected in their algorithm (i.e., the amplitude adjustment is based on {x(jΔ t )} N rather than its normalized version in the IPAC algorithm). This is convenient and may likely speed up the convergence of the algorithm. However, the basis for the shuffling of {x(jΔ t )} N is unclear if the marginal probability distribution of x(jΔ t ) for a nonstationary process is assumed to be time-varying.

Numerical examples
In this section, we illustrate the proposed algorithm by generating surrogate for a given ground motion record and for a given fluctuating component of wind velocity time history of a high-intensity wind event. We apply the algorithm to sample nonstatinary ground motions for prescribed target PSD, where the target is defined based on a set of ground motion records, and the CDF is assumed to be Gaussian and non-Gaussian. The test of the proposed algorithm for esoteric mathematical models is beyond the consideration of the present study.

Generating surrogate for an earthquake record
Consider the record shown in Fig. 2. By applying ST with the window parameter κ = 1 (see Eq. (6)), the obtained TFPSD function is shown in Fig. 3a, and the calculated time-varying σ(t) is presented in Fig. 3b, showing that the TFPSD varies in time and frequency.
By applying the IPAC algorithm, a surrogate is simulated and shown in Fig. 3b. The TFPSD of this surrogate is depicted in Fig. 3c. The figure shows that the surrogate resembles the given record, and its TFPSD function resembles well that shown in Fig. 3a. As x AC (t) is used for generating surrogate (see Step II.5 in IPAC algorithm), the amplitude (or probability distribution) matching is certain, so no plot is provided. Additional test runs indicate that the convergence is usually achieved within five iterations, depending on the adopted convergence criterion. It was noted that the average TFPSD function from multiple generated surrogates tends to be smoother as compared to the TFPSD of the observed record, which is expected since the observed as well as a single sampled record are realizations of stochastic processes. The analysis based on ST is repeated but by applying WT using the generalized morse wavelets (GMWs) [33], where U(ω) is the Heaviside function, a β, γ = 2(eγ/β) β/γ , and β and γ are model parameters. For GMW, C ψ = 2a 2β, γ Γ(2β, γ)/(2γ) and C 1ψ = a β, γ Γ(β, γ)/(γ). The GMW is an analytical wavelet, and it was used to evaluate the coherence of the seismic ground motions [39]. For the numerical analysis, β = 3, γ = 20, c 0 = 0.528, s 0 = 2 0.1 and K = 91 (see Eq. (17)) are employed since these values are suggested as the default values for the algorithm implemented in MATLAB (Version 2019a). The obtained TSPSD and σ(t) of the record are shown in Fig. 4a and b, respectively. A generated surrogate is also shown in Fig. 4b with its corresponding TSPSD function depicted in Fig. 4c. An inspection of the surrogate depicted in Fig. 4b and the original record presented in Fig. 2 indicates that they exhibit similar temporal trends. The TSPSD of the surrogate resembles that of the given record. Again, the convergence is achieved within a few iterations. A comparison of σ(t) shown in Figs. 3 and 4b indicates that they are almost identical. The minor differences between them are due to that ST and WT have different time-frequency (or time-scale) decomposition.

Generating surrogate for a wind record
Now, consider a wind record presented in Fig. 5a. For simplicity, the box window with a width of 32 samples is used to find the mean wind velocity of the time-varying wind record. By removing the mean, the fluctuating component of the wind is presented in Fig. 5b. By applying ST and WT, and following the same analyses that are carried out for the ground motion record shown in the previous section, the obtained results are presented in Figs. 6 and 7. In general, the observations that can be drawn from this example are similar to those presented in the previous section for the ground motion record.

Simulating ground motions
Consider a set of 12 ground motion records oriented in the E-W direction for a seismic event that occurred on January 16, 1986, with a local magnitude of 6.1, focal depth of 10.2 km, and an epicentral distance of 25.2 km. These records are recorded by the LSST array in Lotung,   Taiwan, where the separation between any two recording sites is less than 100 m, as shown in Fig. 8a. Three records from the 15 recording sites seem to be corrupted and are not considered. The record obtained from FA-1 site is illustrated in Fig. 8b. To minimizing the wave passage effect, first, each of the remaining 11 records is time-shifted with respect to the record presented in Fig. 8b such that the sum of the product of a considered record and that shown in Fig. 8b is maximized after the shift. It is assumed that the average PSD of the considered 12 records could provide a good representation of the ground motions, at least for such type of seismic event. The calculated average TFPSD based on ST and the calculated average TSPSD based on WT are shown in Fig.  8c and d, respectively. The calculated σ(t) by using the average TFPSD and the average TSPSD presented in Fig.  8c and d are included in Fig. 8b. The obtained PSD and the standard deviation indicate the nonstationarity of the ground motions. σ(t) values obtained by using ST and WT are in very good agreement.
An assessment of the empirical probability distribution of the standardized variable z(t) = x(t)/σ(t) is carried out by considering all 12 records. The empirical distribution of z(t) by considering all 12 records is presented in Fig. 9, indicating that the empirical distribution can be fitted by a normal distribution only for the initial segment of the records. Moreover, the distribution shape is time-varying and non-Gaussian. The non-Gaussian behaviour of the ground motions is supported by the findings given in Radu and Grigoriu [40], indicating that the Gaussian assumption for the seismic ground motions records included in the PEER NGA-West dataset may not always appropriate. However, for this particular set of records, the tail of the distribution is shorter than that of the normal distribution, which differs from the longer tail behaviour suggested by Radu and Grigoriu [40].
For illustration purposes, it is assumed that the marginal probability density function of x(t) can be represented by the generalized Gaussian distribution (GGD) [29,54], where μ denotes the mean, β 0 and β 1 are positive model parameters, and Γ(·) denotes the gamma function. If β 0 equals 2, it represents the normal distribution. For β 0 > 2 and < 2, the distribution tail is lighter and heavier than that of normal distribution. The variance equals β 2 1 Γð3= β 0 Þ=Γð1=β 0 Þ , and the kurtosis coefficient equals Γ(5/ β 0 )Γ(1/β 0 )/Γ 2 (3/β 0 ).   By considering β 0 = 2 and β 1 ¼ ffiffi ffi 2 p (i.e., standard Gaussian), we use ST and the average target TFPSD function shown in Fig. 8c to sample the records using the IPAC algorithm. Since a comparison of two sampled records is irrelevant for a stochastic process, only a sampled record is illustrated in Fig. 10a. The average TFPSD function obtained from the 1000 sampled records is presented in Fig. 10b, and the calculated spectral acceleration (SA) for a damping ratio of 5% is shown in Fig. 10c for the 1000 sampled records. Similarly, we use WT and the average target TFPSD function shown in Fig. 8d to carry out the simulation. The obtained results are presented in Fig. 10d-f. The PSD functions shown in Fig. 10b and e are almost identical to their corresponding targets presented in Fig. 8c and d. The mean of SA values shown in Fig. 10c and f are in good agreement. The standard deviation of SA obtained by using ST smaller than that obtained by using WT.
We have tested the IPAC algorithm to simulate ground motions for additionally selected target PSD functions. It was observed that in some cases, when the  initial or final segment of records has less than 0.5% of total energy, the algorithm may converge very slowly or may not converge. In such a case, it is suggested that such segments with negligible energy are to be truncated. This is in agreement with common practice in earthquake engineering.
To simulate the non-Gaussian process, we consider f X, ground motions are presented in Fig. 11. A comparison of the results shown in Figs. 10 and 11 indicates that the results follow the same trends. To assess the quantitative differences between the obtained SA based on Gaussian and non-Gaussian assumptions, we calculate the ratio of the mean of SA shown in Fig. 11 (i.e., non-Gaussian case) to its corresponding value shown in Fig. 10 (i.e., Gaussian case). We do the same for the standard deviation of SA. The obtained values are presented in Fig. 12, indicating that the mean and standard deviation of SA for the non-Gaussian case with a lighter tail are smaller than those for the Gaussian case for the vibration period less than about 0.5 s. The decrease in SA by considering non-Gaussian excitation is most noticeable for a shorter vibration period. This is because stiffer structures are more sensitive to peak acceleration values. In general, the ratio based on ST is smoother than that based on WT. Note that we refrained from discussing the ratio between the SA values obtained based on ST and WT since such a comparison could be misleading because the TSPSD and TSPSD used are based only on 12 records and from the same seismic event.
We note that the simulation of a record (Gaussian or non-Gaussian) described in Figs. 10 and 11 is carried out by using a laptop with Intel(R) Core (TM) i7-7700 CPU @ 3.60GHz (6 core and 12 treads). The wall clock time for simulating a record (including I/O) is less than about 0.5 s if ST is used and is less than about 0.15 s if CWT is used. The difference in the computing time by using ST and CWT can be explained by noting that the octave scale is used in CWT.

Summary and conclusions
We elaborated on the concept of defining a modulated and intensity function adjusted (MODIF) stochastic process in the transform domain. The definitions of the stochastic processes in the transform domain lend themselves to an easily understandable and almost trivial algorithm to simulate stochastic processes. As such a simulated signal may not lead to the prescribed target power spectral density function and marginal cumulative distribution function of the process, we proposed a new iterative algorithm, called iterative power and amplitude correction (IPAC) algorithm, so the sampled signal after the iteration have the prescribed properties. Besides simulating nonstationary and non-Gaussian processes, the proposed iterative algorithm can be used to generate surrogate. The algorithm can be used with Fourier transform, S-transform, and continuous wavelet transforms. Practical illustrative numerical examples showed the applicability of the proposed algorithm by sampling surrogates for the ground motions and the fluctuating component of winds. The use of the IPAC algorithm to simulate nonstationary Gaussian and non-Gaussian ground motions based on S-transform (ST) and continuous wavelet transform (WT) is presented. The spectral accelerations are calculated using the simulated records. In general, the mean and standard deviation of SA of the simulated records based on ST and based on WT agree well despite the differences between ST and continuous WT and between the frequency-dependent window used in ST and the generalized Morse wavelet used in the continuous WT.