Sparse polynomial chaos expansions of frequency response functions using stochastic frequency transformation

Let ${\mathcal{M}}$ be a computational model with M-dimensional random inputs $\boldsymbol{X}$=$\{X_{1},X_{2},...,X_{M}\}^{\text{T}}$ and a scalar output $Y$. Further, let us denote the joint probability distribution function (PDF) of the random inputs by $f_{\boldsymbol{X}}(\boldsymbol{x})$ defined in the probability space ($\Omega$,$\mathscr{F}$, $\mathbb{P}$).

Assume that the system response $Y={\mathcal{M}}(\boldsymbol{X})$ is a second-order random variable, i.e. ${\mathbb{E}}\left[Y^{2}\right]<+\infty$ and therefore it belongs to the Hilbert space $\mathscr{H}=\mathscr{L}_{f_{\boldsymbol{X}}}^{2}(\mathbb{R}^{M},\mathbb{R})$ of $f_{\boldsymbol{X}}$-square integrable functions of $\boldsymbol{X}$ with respect to the inner product:

where ${\mathcal{D}}_{\boldsymbol{X}}$ is the support of $\boldsymbol{X}$. Further assume that the input variables are independent, i.e. $f_{\boldsymbol{X}}(\boldsymbol{x})=\prod_{i=1}^{\emph{M}}\emph{f}_{X_{i}}(x_{i})$. Then the generalized polynomial chaos representation of $Y$ reads (Xiu and Karniadakis, 2002):

in which $\tilde{u}_{\boldsymbol{\alpha}}$ is a set of unknown deterministic coefficients, $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},...,\alpha_{M})$ is a multi-index set which indicates the polynomial degree of $\psi_{\boldsymbol{\alpha}}(\boldsymbol{X})$ in each of the M input variables. $\psi_{\boldsymbol{\alpha}}$s are multivariate orthonormal polynomials with respect to the joint PDF $f_{\boldsymbol{X}}(\boldsymbol{x})$, i.e. :

where $\delta_{\boldsymbol{\alpha}\boldsymbol{\beta}}$ is the Kronecker delta. Since the input variables are assumed to be independent, these multivariate polynomials can be constructed by a tensorization of univariate orthonormal polynomials with respect to the marginal PDFs, i.e. ${\psi}_{\boldsymbol{\alpha}}(\boldsymbol{X})=\prod_{i=1}^{M}\psi_{\alpha_{i}}^{(i)}({X}_{i})$. For instance, if the inputs are standard normal or uniform variables, the corresponding univariate polynomials are Hermite or Legendre polynomials, respectively.

In practice, the infinite series in Eq. (6) has to be truncated. Given a maximum polynomial degree $p$, the standard truncation scheme includes all polynomials corresponding to the set ${\mathcal{A}}^{M,p}=\{\boldsymbol{\alpha}\in{\mathbb{N}}^{M}\;\colon\;|\boldsymbol{\alpha}|\leq p\},$ where $|\boldsymbol{\alpha}|=\sum_{i=1}^{M}\alpha_{i}$ is the total degree of polynomial $\psi_{\boldsymbol{\alpha}}$. The cardinality of the set ${\mathcal{A}}^{M,p}=\binom{M+p}{p}=P$ increases rapidly by increasing the number of parameters M and the order of polynomials $p$. However, it can be controlled with suitable truncation strategies such as $q$-norm hyperbolic truncation (Blatman and Sudret, 2010), that drastically reduce the number of unknowns when $M$ is large.

The estimation of the vector of coefficient $\tilde{u}_{\boldsymbol{\alpha}}$ can be done non-intrusively by projection (Ghanem and Ghiocel, 1998; Ghiocel and Ghanem, 2002) or least square regression methods (Blatman and Sudret, 2010; Berveiller et al., 2006). The latter is based on minimizing the truncation error $\epsilon$ via least square as follows:

This can be formulated as

Let $\boldsymbol{{\mathcal{X}}}=\{\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)},...,\boldsymbol{x}^{(N_{ED})}\}$ and $\boldsymbol{{\mathcal{Y}}}=\{{y}^{(1)}={\mathcal{M}}(\boldsymbol{x}^{(1)}),{y}^{(2)}={\mathcal{M}}(\boldsymbol{x}^{(2)}),...,{y}^{(N_{ED})}={\mathcal{M}}(\boldsymbol{x}^{(N_{ED})})\}$ be an experimental design with $N_{ED}$ space-filling samples of $\boldsymbol{X}$ and the corresponding system responses, respectively. Then, the minimization problem (9) admits a closed form solution

in which $\boldsymbol{\Psi}$ is the matrix containing the evaluations of the Hilbertian bases, that is $\boldsymbol{\Psi}_{ij}={\psi}_{\boldsymbol{\alpha}_{j}}(\boldsymbol{x}^{(i)}),i=1,2,...,N_{ED},j=1,2,...,P$.

The accuracy of PCE will be improved by reducing the effect of over-fitting in least square regression. This can be done by using sparse adaptive regression algorithms proposed in Hastie et al. (2007); Efron et al. (2004). In particular, the Least Angle Regression (LAR) algorithm has been demonstrated to be effective in the context of PCE by Blatman and Sudret (2011a).

In the case of vector-valued response, i.e. $\boldsymbol{Y}\in\mathbb{R}^{N},N>1$, the presented approach may be applied componentwise. This can make the algorithm computationally cumbersome for models with large number of random outputs. To decrease the computational cost, one can extract the main statistical features of the vector random response by principal component analysis (PCA). The concept has been adapted to the context of PCE by Blatman and Sudret (2013).

To perform sample-based PCA, let us rewrite the $\boldsymbol{{\mathcal{Y}}}$ as the combination of its mean $\bar{\boldsymbol{{\mathcal{Y}}}}$ and covariance matrix as follows:

where the $\boldsymbol{v}_{i}$’s are the eigenvectors of the covariance matrix:

and the $\boldsymbol{u}_{i}$’s are vectors such that

One can approximate $\boldsymbol{{\mathcal{Y}}}$ by the $\hat{N}$-term truncation:

Since $\bar{\boldsymbol{{\mathcal{Y}}}}$ and $\boldsymbol{v}^{\text{T}}_{i}$ are the mean and the eigenvectors of the system responses, respectively, they are independent of the realization. Therefore, PC expansion can be applied directly on the $\hat{N}\ll N$ auxiliary variables $\boldsymbol{u}_{i}$. Besides, acknowledging the fact that the PCA is an invertible transform, the original output can be retrieved directly from Eq. (14) for every new prediction of $\boldsymbol{u}$.

In this section, the method of stochastic frequency transformation is developed to address the challenge of frequency shift at eigenfrequencies due to uncertainty in the parameters. The idea is basically to apply a transformation to the system responses to maximize their similarity before building PCE, as first proposed by Mai and Sudret (2015). Here, the technique is extended and adopted into the frequency domain to obtain PCEs of the FRFs.

To this end, the following algorithm is proposed. First, an experimental design $\boldsymbol{\mathcal{X}}$ and the corresponding model responses $\boldsymbol{{\mathcal{Y}}}$ are evaluated. Each system response will be called a trajectory in the remainder of this paper. Let the frequency range of interest be discretized to $n_{\omega}$ equidistant frequencies $\Omega_{d}=\left[\omega_{1},\omega_{2},...,\omega_{n_{\omega}}\right]$. Then, the required system responses are matrices $\boldsymbol{\mathcal{H}}(\Omega_{d})\in\mathbb{C}^{n_{u}n_{y}\times n_{\omega}}$ and $\boldsymbol{\boldsymbol{\mathcal{F}}}\in\mathbb{R}^{n_{u}n_{y}\times n_{sf}}$. The matrix $\boldsymbol{\mathcal{H}}(\Omega_{d})$ is obtained by evaluating Eq. (4) at frequencies $\Omega_{d}$. The matrix $\boldsymbol{\boldsymbol{\mathcal{F}}}$ consists of all the resonance and antiresonance frequencies of the system’s input-output relations for one system realization, as follows:

in which $n_{p}$ is the number of eigenvalues of the system. Furthermore, $\{\omega_{p_{i}}i=1,2,...,n_{p}\}$ are the resonant frequencies and $\{\omega^{l}_{m_{i}}i=1,2,...,n_{p}-1,l=1,2,...,n_{u}\times n_{y}\}$ are frequencies between each two consecutive resonant frequencies at which the minimum amplitude occurs. Throughout the paper, these important frequencies, shown by red asterisks in Figure 1 for a typical frequency response, will be referred to as selected frequencies. Their number $n_{sf}$ is assumed to be constant across different realizations of the system inputs. $\{\mathcal{F}_{i},i=1,2,\cdots,n_{u}\times n_{y}\}$ includes all the selected frequencies for the $i^{th}$ input-output relation.

For the next step of the algorithm, let $\boldsymbol{x}^{(ref)}$ be selected randomly among the sample points in the ED to have its associated trajectory as the reference, i.e. :

Then, the other trajectories are transformed in the frequency axis so as to have the peaks and valleys as close to the corresponding locations in the reference trajectory as possible i.e. :

where $i=1,2,\cdots,n_{u}\times n_{y}$, $k=1,2,\cdots,N_{ED}$ and $\nu$ is the transformed frequency axis called scaled frequency. The transform $\mathcal{T}^{k}_{i}$ consists of a continuous piecewise linear transform of the intervals between the identified selected frequencies that align them to the corresponding ones of the reference trajectory as follows,

where

$j=1,2,\cdots,n_{sf}-1$ and $l=1,2,\cdots,n_{\omega}$. This transformation results in the FRFs which are similar to the reference one in the scaled frequency domain:

where the set ${\mathcal{N}}^{(k)}_{i}=\{\nu_{i,1}^{(k)},\nu_{i,2}^{(k)},\cdots,\nu_{i,n_{\omega}}^{(k)}\}$ consists of the discretized scaled frequencies which are non-equidistantly spread over the frequency range of interest.

Figures 2(a) and 2(b) illustrate the FRFs of a 2-DOF system versus frequency and scaled frequency, respectively. An example of such a transform used for transforming the FRFs of a 2-DOF is presented in Figure 3.

The non-smooth behavior of the FRFs makes their surrogation by polynomials a problematic task. To solve this issue, one PCE could be calculated for each scaled frequency. This means that to compute PCEs for the FRFs, two sets of PCEs are required. The first set is to predict the selected frequencies, collected in the matrix $\boldsymbol{\mathcal{F}}$ (16), which are required for performing the stochastic transformation as explained in Section 3.3. This matrix includes eigenfrequencies of the system, therefore by obtaining this set of PCEs, the problem of random eigenvalue calculation is solved by the use of PCE as a byproduct. This problem has been addressed in some recent works, e.g. Pichler et al. (2012). Since the number of random outputs for this set is not very large, PCE can be applied to each of the selected frequencies separately, i.e. for $i=1,2,\cdots,n_{u}\times n_{y}$ and $j=1,2,\cdots\times n_{sf}$

The second set of PCE is for the system response at each individual scaled frequency. To this end, let $\tilde{\bf H}({\mathcal{N}}^{ref})\in\mathbb{C}^{N_{ED}\times(n_{\omega}\times n_{u}\times n_{y})}$, defined in Algorithm 1, be a matrix of trajectories at scaled frequencies. Since the FRFs have complex-valued responses, whereas the PCEs are defined for real-valued functions only¹¹1Limited literature is available on the use of PCE for complex-valued functions, see e.g. Soize and Ghanem (2004)., separate PCEs need to be performed for real and imaginary parts of the FRFs. Therefore, the matrix $\boldsymbol{\mathcal{G}}=\{\boldsymbol{\mathcal{G}}^{\mathfrak{R}},\boldsymbol{\mathcal{G}}^{\mathfrak{I}}\}=\{real(\tilde{\bf H}({\mathcal{N}}^{ref})),imag(\tilde{\bf H}({\mathcal{N}}^{ref}))\}\in\mathbb{R}^{N_{ED}\times(2\times n_{\omega}\times n_{u}\times n_{y})}$ is the response matrix for which the PCE should be built.

The number of random outputs for this set, $N=2\times n_{\omega}\times n_{u}\times n_{y}$, can be extremely large. As discussed in Section 3.2, the PCEs are therefore applied directly to the principal components of $\mathcal{G}$, yielding:

where $u_{\boldsymbol{\alpha}}^{\mathfrak{R}}$ and $u_{\boldsymbol{\alpha}}^{\mathfrak{I}}$ are respectively the vectors of coefficients of the PCEs made for the real and imaginary parts of the FRFs.

To predict the surrogate model response at a new sample point $\boldsymbol{x}^{(0)}$, several steps need to be taken to transform the PCE predictors in Eqs. (21) and (22) from the scaled frequency axis $\nu$ to the original frequency axis $\omega$. The matrices $\hat{\boldsymbol{\mathcal{G}}}^{\mathfrak{R}}$ and $\hat{\boldsymbol{\mathcal{G}}}^{\mathfrak{I}}$ are obtained by evaluating the second set of PCEs in Eqs. (21) and (22), respectively. Then, the FRFs at the scaled frequencies can be obtained at the new sample point by the inverse vectorization of ${\hat{\tilde{\bf{H}}}}({\mathcal{N}}^{ref})=\hat{\boldsymbol{\mathcal{G}}}^{\mathfrak{R}}+j\hat{\boldsymbol{\mathcal{G}}}^{\mathfrak{I}}$ where $j=\sqrt{-1}$. To obtain the FRF at the original frequency $\omega$ the following transformation is used,

where $\mathcal{T}^{0}_{i}$ is obtained by evaluating Eq. (18) at $\hat{\boldsymbol{\mathcal{F}}}^{(0)}=\hat{\boldsymbol{\mathcal{F}}}(\boldsymbol{x}^{(0)};\omega)$ which is the matrix of selected frequencies at the new sample point $\boldsymbol{x}^{(0)}$ evaluated by Eq. (20). Besides, $\Omega_{i}^{(0)}$ is a set of discretized frequencies which are non-equidistantly spread over the frequency rage of interest. In order to provide the frequency response at the desired frequencies $\Omega_{d}$, interpolation is inevitable. The algorithm for predicting the system response at a new sample point is briefly presented in Algorithm 2.

In this section, the proposed method will be applied to two case studies. The first one is a simple 2-DOF system to illustrate how the method works. The second one is a 6-DOF system with a relatively large (16-dimensional), parameter space. For the sake of readability, only results for one output (the $1^{st}$ output for the 2-DOF and the $6^{th}$ output for the 6-DOF system) are shown for each case while the results for the other outputs are reported for completeness in the Appendices.

To assess the accuracy of the surrogate models quantitatively, the following measure based on the root mean square (rms) error of the vectors is defined.

in which $(\bullet)$ is the vector of interest. $(\cdot)^{ex}$ and $(\cdot)^{approx}$ represent results obtained by running the true and surrogate models, respectively. This error aims at measuring the relative difference between vectorial data, such as one FRF or the mean and standard deviation of several FRFs. For the mean and standard deviation of the data, the reference results are obtained by evaluating the true model at 10,000 Monte-Carlo samples and the approximations are calculated by the PCE surrogate at the same 10,000 points.

As the first example, the simple 2-DOF system shown in Figure 4 is selected to highlight the steps of the proposed method. In this system, stiffness is assumed to be uncertain

where $\xi$ is a standard normal random variable. Other properties of the system are listed in Table 1. The system has one input force $f$ at mass 1, two physical outputs $q_{1}$ and $q_{2}$ and thus, two FRFs. The FRFs of the system are obtained in the range of 10 to 35 Hz with a frequency step of 0.01 Hz, as shown in Figure 1. The selected frequencies are also shown in the figure with red asterisks.

40 points are sampled in the parameter space using Latin Hypercube Sampling (LHS) to form an experimental design (ED) ${\boldsymbol{{\mathcal{X}}}}$ and the model is evaluated at these points to find the system responses of interest, namely the FRFs, $\boldsymbol{\mathcal{H}}$, and the selected frequencies $\boldsymbol{\mathcal{F}}$.

Figure 2(a) shows the FRFs of the system evaluated at ${\boldsymbol{{\mathcal{X}}}}$. To find the transformed FRFs, i.e. FRFs at scaled frequencies $\nu$, one trajectory was selected randomly as the reference and the others were scaled such that their peaks and valleys were at the same scaled frequencies as that of the reference trajectory. The transformed FRFs are shown in Figure 2(b) and the corresponding continuous piecewise-linear transformations in Figure 3.

The next step is to find a suitable basis and the associated coefficients for the polynomial chaos expansion. In this case, since the random variable is Gaussian, the basis of the polynomial chaos consists of Hermite polynomials. The LARS algorithm (Blatman and Sudret, 2011b) is employed here to calculated a sparse PCE with adaptive degree.

The first set of expansion consists of 10 PCEs to surrogate the selected frequencies $\boldsymbol{\mathcal{F}}$, shown in Figure 1 by red asterisks. As the second set of expansions, PCE is made for the dominant components of $\boldsymbol{\mathcal{G}}$ as explained in Section 3.2.1. To do so, the $\hat{N}$ largest principal components are selected such that the sum of their associated eigenvalues amounts to $99\%$ of the sum of all the eigenvalues, i.e. $\sum_{i=1}^{\hat{N}}\lambda_{i}=0.99\sum_{i=1}^{{N}}\lambda_{i}$ in which $\lambda_{i}$’s are the eigenvalues of the covariance matrix of either $\boldsymbol{\mathcal{G}}^{\mathfrak{R}}$ or $\boldsymbol{\mathcal{G}}^{\mathfrak{I}}$. By this truncation, the number of random outputs is reduced from 2501 $\times$ 2 $\times$ 2 to 6 components, namely 3 components for the real part and 3 components for the imaginary part. The spectra of the eigenvalues of the covariance matrix of the $\boldsymbol{\mathcal{G}}^{\mathfrak{R}}$ and $\boldsymbol{\mathcal{G}}^{\mathfrak{I}}$ are displayed in Figure 5. All the PCEs used to make this surrogate model, including those for the selected frequencies and those for the dominant components have orders between 3 and 8.

The efficiency of the proposed approach is assessed by comparing the prediction accuracy of the surrogate model on a large reference validation set (10,000 samples calculated with the full model). PCE estimate of the mean and standard deviation of the surrogate model are compared to their Monte-Carlo estimators on experimental designs of increasing size. The resulting convergence curves are given in Figure 6.

They indicate that the PCE converges faster to the reference results for the mean and standard deviation. Their estimates are approximately two and one order of magnitude more accurate than those from the MC estimators. In addition, one can conclude that 40 points are enough for the ED in this example, since for larger sizes the accuracy does not improve significantly.

It is worth mentioning that the coefficient of variation (COV) of the parameters and the level of damping are among the criteria that can affect the size of the ED. Therefore, larger COV and lower levels of damping are not obstacles for the proposed method provided a sufficiently large ED is used.

The 10,000 model evaluations used to produce the convergence curves in Figure 6 are also used to provide a detailed validation of the performance of the PCE surrogates on various quantities of interests. The results are presented in the following figures. In Figure 7, the selected frequencies obtained by the PCE are shown versus the true ones. The results show that the PCE model accurately predicts the selected frequencies. Since the amplitudes at the resonant frequencies are the most sensitive parts of the FRFs and contains most of crucial information about the system, it is one of the most interesting parts of the FRFs for the researcher. Figure 8 shows the histograms of these amplitudes obtained by the true and the surrogate models at at all 10,000 validation points. Moreover, in Figure 9, the whole FRFs at the validation points are depicted and compared. Their associated means are also shown by black lines. The results indicate accurate prediction of the amplitude at the resonant frequencies as well as the whole FRFs. A quantitative accuracy analysis for the whole FRF was done by using Eq. (24) and its histogram presented in Figure 10 confirms the high accuracy of the surrogate model.

Another accuracy test is given by the comparison between the first two moments of the FRFs. The mean and standard deviation of the trajectories obtained by the true model and the surrogate model are compared in Figure 11 for the $1^{st}$ output and in A.1 for the $2^{nd}$ output. They reveal the accuracy of the proposed surrogate model in predicting the first two moments of the FRFs.

To show the feasibility of the proposed method to estimate the statistics of the FRFs, the results obtained here are compared to their counterparts in two of the most recent works available in the literature. The first study (Jacquelin et al., 2015a) directly uses high-order PCE for estimating the first two moments of the FRF, whereas the second method (Jacquelin et al., 2015b) proposes to use Aitken’s transformation in conjunction with PCEs, Both methods use PCEs of order 50 and tend to produce spurious peaks around the resonance region. The use of Aitken’s transformation slightly improves convergence. Their results for the mean and standard deviation are shown in Figures 12(a) and 12(b), respectively. For comparison, the results from our approach in Figures 11(a) and 11(c) are reproduced in Figure 12 with a scaling similar to the other panels. They indicate that the stochastic frequency transformation approach proposed here, significantly improves the estimation accuracy of the PCE surrogate, as no spurious peaks are visible in this case.

As far as individual comparison between the true and the predicted FRF are concerned, the worst predicted FRF, the one with the maximum error, is presented in Figure 13. They indicate that even for the worst-case, the presented approach results in prediction of the FRFs with excellent accuracy.

The second example is chosen to illustrate the application of the proposed method to a problem with a relatively large parameter space. The system, shown in Figure 14, consists of 10 springs and 6 masses which are modeled by random variables with lognormal distributions. Their mean values are listed in Table 2. The uncertainty on the springs (resp. masses) has a COV = 10% (resp. COV = 5%). The damping matrix is $\boldsymbol{V}=0.1\widehat{\boldsymbol{M}}$, where $\widehat{\boldsymbol{M}}$ is the matrix of the mean value of the system masses. Table 3 provides its corresponding mean of modal dampings evaluated over 10,000 samples.

The system has one input force at mass 6 and 6 system outputs, one for each mass. The FRF of the system is evaluated at a frequency range from 1 to 25 rad/s with the step of 0.01$\pi$ rad/s.

In this example, the ED consists of 400 points sampled from the parameter space using LHS. The marginal distributions of the input vectors $\boldsymbol{X}$ consists of lognormal distributions. Therefore, the chosen PCE basis consists of Hermite polynomials on the reduced variable ${\boldsymbol{Z}}=\ln({\boldsymbol{X}})$. Eq. (6) thus, can be written as

The LARS algorithm has been employed to build sparse PCEs with adaptive degree for both the selected frequencies and the principal components of the scaled FRF.

For the second set of PCEs, PCA has been performed and the dominant components are selected such that $\sum_{i=1}^{\hat{N}}\lambda_{i}=0.999\sum_{i=1}^{{N}}\lambda_{i}$. This truncation reduced the number of random outputs from 761 $\times$ 6 $\times$ 2 to 102 components. Since the dimension of the input parameter space is large, to reduce the unknown coefficients of the PCEs and avoid the curse of dimensionality, a hyperbolic truncation with $q$-norm of 0.7 was used before the LARS algorithm. Besides, only polynomials up to rank 2 were selected here (i.e. polynomials that depend at most on 2 of the 16 parameters). It should be mentioned that all the PCEs used for the surrogate model have eventually maximum degrees less than 10.

The efficiency of the proposed method is assessed by comparing the PCE estimates of the first two moments of the surrogate model with the plain Monte-Carlo estimators on experimental designs of increasing size. The reference validation set is obtained by 10,000 points sampled from the parameter space by LHS at which the full model is evaluated. The results are shown in Figure 15 for the mean and standard deviation at $6^{th}$ output. The results for the other outputs are presented in B.1. They indicate that both mean and standard deviation evaluated by the surrogate model converge faster than those of the Monte-Carlo simulations. Besides, it can be inferred that 400 points are enough for the ED.

In order to assess the accuracy of the surrogate model in estimating various quantities of interests, the same 10,000 points used as the reference validation set to study the convergence are used here. Figure 16 illustrates two of the predicted selected frequencies versus the true ones, namely the best and worst predicted eigenfrequency, so that the accuracy of the surrogate model in this step can be inferred. While the overall accuracy is very good for all frequencies, it tends to degrade somewhat at higher frequencies.

Besides, at all the validation points the FRFs are calculated by both the true model and the surrogate model. The variation of the amplitudes at the first and fifth resonant frequencies are shown as histograms in Figure 17. Plots of the individual FRFs are reported in Figure 18 for $6^{th}$ output. In order to assess the error quantitatively, each response of the surrogate model has been compared with the corresponding one of the true model in the root-mean-square sense. This error is evaluated using Eq. (24) and the corresponding results are presented in Figure 19. They indicate the high accuracy of the proposed surrogate model in predicting the FRFs.

As an individual comparison between the true FRFs and predicted by the surrogate model, two cases are considered: one case with an average and one with the maximum overall error, are selected and their $6^{th}$ output are demonstrated in Figure 20. The other outputs are presented in B.3.

The mean and standard deviations of the FRFs were compared with the reference ones. The results for $6^{th}$ output are plotted in Figure 21. The other outputs are plotted in B.2. They indicate that although both the mean and the standard deviation match with their reference, the standard deviation presents minor mismatch at the peaks.