## 1 Introduction  
Flexible long-span bridges are susceptible to extreme wind. As one of the wind-induced structural responses, bridge buffeting has been the concern of engineers for a long time [1]. The buffeting response, which is induced by turbulent winds, can be analyzed in both time [2] and frequency domain [3]. It is noted that the linear hypothesis is the prerequisite of the buffeting analysis in the frequency domain. Since the aerodynamic and geometrical effects may cause flexibilities and non-linear behaviors in long-span bridges, frequency-domain buffeting analysis is limited to most of long-span bridges [4; 5; 6]. As for time domain analysis, numerical integration is selected to solve the buffeting responses, hence, the non-linear and aerodynamic behaviors of long-span bridges can be considered [7; 8]. To conduct time domain analysis, it is significant to accurately simulate the turbulent wind fields of wind-sensitive structures, since wind-induced responses are sensitive to the simulated wind loads [9].  
The spectral representation method (SRM), which is an accurate and simplified method, is widely employed to generate the turbulent wind fields in structural wind engineering [5; 10]. According to the stationary hypothesis, the turbulence is typically regarded as a zero-mean stationary stochastic process over a given time, which is determined by the power spectral density (PSD) [11; 12]. However, field measurements indicate that extreme wind shows strong non-stationary features [13; 14]. Therefore,it is significant to conduct non-stationary wind field simulation for analyzing non-stationary buffeting in the time domain. [PERSON] [15] proposed the general SRM-based framework to generate the non-stationary and multivariate stochastic processes, in which the evolutionary power spectral density (EPSD) was utilized to characterize the non-stationary processes. It is revealed in this work that the fast Fourier transform (FFT) cannot be implemented due to the time-dependent EPSD matrix. Thus, the simulation efficiency is restricted when using the general non-stationary SRM.  
To take advantage of the FFT technique, some efforts have been made. [PERSON] and [PERSON] [16] proposed the FFT-aided SRM, in which the decomposed EPSD matrix is fitted by various polynomial functions. [PERSON] [17] introduced the proper orthogonal decomposition to decompose the time-frequency EPSD, which leads to a higher simulation efficiency. [PERSON] et al. [18] used the Taylor series expansion to enhance the simulation efficiency of the non-stationary SRM. [PERSON] et al. [19] also developed a FFT-aided SRM to simulate the non-stationary turbulent wind field of the bridge deck based on the non-negative matrix factorization (NMF). Compared with other FFT-aided approaches, the NMF-based method ensures the physical meaning of the factorized low rank matrices, and achieves straightforwardness in the simulation. However, the existing NMF-based method can only simulate horizontal non-stationary turbulence, which is a flawed approach when the vertical turbulent wind field must be considered.  
In this study, the wavelet-based method is utilized to estimate both the along-wind and vertical EPSDs using the data measured at the Runyang Suspension Bridge (RSB) site during Typhoon Wipha. The existing NMF-based FFT-aided SRM, which cannot simulate the vertical turbulent wind field, is updated to consider both the horizontal and vertical ones. The non-stationary turbulent wind fields for the girder and tower of the RSB are generated using the updated NMF-based SRM. All the program was executed in Matlab due to its convenience in matrix computation. The reliability of the updated method is verified by comparing the estimated correlation functions and the target ones.