## 3 Enhanced Interpolation Approach
### Decomposition of EPSD Matrix Based on Interpolation  
The phases \(\theta_{jk}(\omega)\) generally admitted the requirement that \(\theta_{jk}(\omega)=\theta_{jl}(\omega)-\theta_{kl}(\omega)\) for the ground motion [1]. Thus, the decomposition of the complex matrix \(\mathbf{S}(\omega,t)\) could be simplified as the decomposition of its real matrix \(|\mathbf{S}(\omega,t)|\) to obtain an immediate decomposition [12], i.e.,  
\[|\mathbf{S}(\omega,t)|=|\mathbf{H}(\omega,t)||\mathbf{H}(\omega,t)|^{T} \tag{14}\]  
\[\theta_{jk}(\omega)=\theta_{jk}(\omega),j>k \tag{15}\]  
where  
\[|\mathbf{\Gamma}(\omega)||\mathbf{S}(\omega,t)|=\begin{bmatrix}S_{11}(\omega,t)&|S_{12}(\omega,t)|&\cdots&|S_{1n}(\omega,t)|\\ |S_{21}(\omega,t)|&S_{22}(\omega,t)&\cdots&|S_{2n}(\omega,t)|\\ \vdots&\vdots&\ddots&\vdots\\ |S_{n1}(\omega,t)|&|S_{n2}(\omega,t)|&\cdots&S_{nn}(\omega,t)\end{bmatrix} \tag{16}\]  
\[|\mathbf{H}(\omega,t)|=\begin{bmatrix}H_{11}(\omega,t)&0&\cdots&0\\ |H_{21}(\omega,t)|&H_{22}(\omega,t)&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ |H_{n1}(\omega,t)|&|H_{n2}(\omega,t)|&\cdots&H_{nn}(\omega,t)\end{bmatrix} \tag{17}\]  
Based on Equations (2) and (3), \(|\mathbf{S}(\omega,t)|\) can be further represented as  
\[|\mathbf{S}(\omega,t)|=\mathbf{D}(\omega,t)|\mathbf{\Gamma}(\omega)|\mathbf{D }^{\mathrm{T}}(\omega,t) \tag{18}\]  
where  
\[\mathbf{D}(\omega,t){=}\text{diag}[\sqrt{S_{11}(\omega,t)},\sqrt{S_{22}( \omega,t)},\cdots,\sqrt{S_{nn}(\omega,t)}] \tag{19}\]  
and the lagged coherence matrix is expressed by  
\[|\mathbf{\Gamma}(\omega)|=\begin{bmatrix}1&|\gamma_{12}(\omega)|&\cdots&| \gamma_{1n}(\omega)|\\ |\gamma_{21}(\omega)|&1&\cdots&|\gamma_{2n}(\omega)|\\ \vdots&\vdots&\ddots&\vdots\\ |\gamma_{n1}(\omega)|&|\gamma_{n2}(\omega)|&\cdots&1\end{bmatrix} \tag{20}\]  
Because \(|\mathbf{\Gamma}(\omega)|\) is also a Hermitian matrix, it can be decomposed as  
\[|\mathbf{\Gamma}(\omega)|=\mathbf{B}(\omega)\mathbf{B}^{\mathrm{T}}(\omega) \tag{21}\]  
\[\mathbf{B}(\omega)=\begin{bmatrix}\beta_{11}(\omega)&0&\cdots&0\\ \beta_{21}(\omega)&\beta_{22}(\omega)&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ \beta_{n1}(\omega)&\beta_{n2}(\omega)&\cdots&\beta_{nn}(\omega)\end{bmatrix} \tag{22}\]  
On the basis of Equations (14), (18) and (21), we had  
\[|\mathbf{H}(\omega,t)|=\mathbf{D}(\omega,t)\mathbf{B}(\omega) \tag{23}\]To sum up, the decomposed EPSD \(H_{jk}(\omega,t)\) can be obtained by  
\[H_{jk}(\omega,t)=\sqrt{S_{jj}(\omega,t)}\beta_{jk}(\omega)e^{i\theta_{jk}(\omega)} \tag{24}\]  
Obviously, the Cholesky decomposition along time and frequency has been simplified as that of along frequency.  
As shown in Equation (21), the Cholesky decomposition of the lagged coherence matrix had to be performed at each frequency. When a large number of ground motion simulation points were involved, the calculation became prohibitive. As a result, the interpolation approximation strategy is commonly utilized to decrease the Cholesky decomposition times. Considering that the elements of \(|\mathbf{\Gamma}(\omega)|\) are uniformly varying and continuous functions of \(\omega\) for the seismic motion, the uniformly distributed frequency interpolation points \(\hat{\omega}\), are selected and calculated with  
\[\hat{\omega}_{r}=(\omega_{u}-\omega_{1})\bigg{(}\frac{r-1}{N_{\hat{\omega}}-1 }\bigg{)}+\omega_{1}\ \ r=1,2,\cdots,N_{\hat{\omega}} \tag{25}\]  
where \(N_{\hat{\omega}}\) is the frequency interpolation point number; \(N_{\hat{\omega}}\ll N\). As long as the decomposition results at the interpolation points \(\mathbf{B}(\hat{\omega})\) are acquired, the interpolation technique, such as the spline interpolation, Lagrange interpolation and Hermite interpolation, may be used to obtain the results at additional frequency points. The Hermite interpolation was utilized in this study because it ensured continuity and differentiability at the interpolation locations.