Simpliﬁed Analysis for Multiple Input Systems: A Toolbox Study Illustrated on F-16 Measurements

: This paper introduces a nonparametric, nonlinear system identiﬁcation toolbox called SAMI (simpliﬁed analysis for multiple input systems) developed for industrial measurements of vibro-acoustic systems with multiple inputs. It addresses the questions related to the user-friendly (semi-)automatic processing of multiple-input, multiple-output measurements with respect to the design of an experiment and the analysis of the measured data. When the proposed toolbox is used, with minimal user interaction, it is easily possible (a) to decide whether the underlying system is linear or not, (b) to decide whether the linear framework is still adequate to be used, and (c) to tell an inexperienced user how much can be gained using an advanced nonlinear framework. The toolbox is illustrated on openly accessible F-16 ground vibration testing measurements.


Introduction
The goal of this paper is to introduce a simple but efficient nonlinear, nonparametric toolbox for the automated analysis of vibro-acoustic industrial measurements. The goal of vibration testing is to obtain experimental data for a whole vibrating structure such as a road or air vehicle. Using these data, it is possible to validate and improve the dynamic models of systems under test. Since the development of advanced digital signal processing, it became possible to use a large variety of user-defined excitation signals to experimentally determine the broadband frequency response functions (FRFs), which are required to obtain parametric models (e.g., resonance frequencies, modes shapes). In this paper, special multisines (also known as pseudo-random noise signals) are considered excitation signals. The main advantage of the proposed signals is that there is no problem with spectral leakage or transience, and they deliver excellent linear models while providing useful information about the measurement.
Many vibro-acoustic structures are inherently nonlinear. As coping with nonlinear systems became increasingly important, numerous modeling methods were proposed. For a general overview of the nonlinear modeling techniques, we refer to [1][2][3][4]. For complexity reasons, nonlinear systems are often approximated with linear systems, because that theory is well-understood, and the cost of nonlinear modelling can be very high. However, the loss of accuracy by this decision is rarely studied. To this end, using the proposed toolbox, it is easy to tell how much can be gained by using an advanced nonlinear modeling technique in order to make a well-balanced decision.
The state-of-the-art modeling tool known as the best linear approximation (BLA) framework is already available for single-input, single-output (SISO) systems [5][6][7][8][9]. Unlike in the SISO case, the design and analysis of a multiple-input, multiple-output (MIMO) experiment are more involved, Vibration 2020, 3, 70-84; doi:10.3390/vibration3020007 www.mdpi.com/journal/vibration and the complexity grows with the number of input and output channels. Further, to properly analyze the estimation results and to troubleshoot, an experienced user is needed. This work extends the BLA methodology to MIMO systems with a large number of input and output channels, and it provides an automated, user-friendly interpretation of the measurement data by extracting the user-relevant information, thereby reducing the interactions from the few days needed to evaluate hundreds of transfer functions manually to a couple of minutes. This automatic processing includes, amongst other things, MIMO excitation signal design, transient analysis, estimation method choice, highlighting strategic information (e.g., nonlinearity detection), and warning the system (e.g., indication of sensory faults). This paper is organized as follows. Section 2 briefly describes the considered systems, assumptions, and BLA framework applied in this work. Section 3 introduces the main components of the toolbox. Section 4 provides an experimental illustration of the ground vibration testing of an F-16 aircraft. Conclusions can be found in Section 5.

Definitions and Assumptions
The considered systems are electronic, mechanical, or civil dynamic vibrating structures. The dynamics of a linear MIMO system can be nonparametrically characterized in the frequency domain by its frequency response matrix (FRM, a matrix whose elements are FRFs [10]) G at frequency index k, which relates n i inputs U, to n o outputs Y obtained from data records with a length of N: where G[k] ∈ C n o ×n i , Y[k] ∈ C n o ×1 , U[k] ∈ C n i ×1 , k = 0 . . . (N − 1)/2 at frequency f k = k f s /N with sampling frequency of f s . In order to make the text more accessible, the frequency indices and dimensionalities will be omitted. The system represented by G is nonlinear when the superposition principle is not satisfied in steady state [11]. In this work an arbitrary number of input and output channels is considered. The underlying systems are damped, bounded-input, bounded-output stable, time-invariant, nonlinear systems where the linear response of the system is still present, and the output of the underlying system has the same period as the excitation signal (i.e., the system has PISPO behavior: period in, same period out [11]). The output is measured with additive, independent and identically distributed Gaussian noise (denoted by E) with zero mean and finite variance, and it can contain (a smooth, decaying) transient response (denoted by T) such that the measurement Y measured is given by:

Measurement and Instrumentation
The toolbox relies on the use of classical instrumentation and measurement setups [12]. The reference (excitation) signal is an ideal multisine (pseudo-random noise) generated by the toolbox. Multisines are generated in the frequency domain: the magnitudes are freely set by the user, but the phases are randomly chosen [13]. As a result, the signals have a Gaussian distribution. The excitation system (e.g., a shaker) receives the reference signal (e.g., voltage signal) and generates the input signal (e.g., force) to the system (e.g., a vehicle). The output signals (e.g., accelerations) contain the waveforms of the system's reaction to the input signal. It is important to stress that the toolbox requires the measured time domain signals to be stored.
From an instrumentation point of view, it is assumed that the measurement system is perfectly synchronized, the sampling frequency is kept constant, the input signal is measured precisely such that the signal-to-noise ratio (SNR) is at least 20 dB; otherwise, the FRM estimate will be biased. However, if the reference signal is available, then the toolbox can compensate for the bias by the appropriate choice of estimator [13][14][15].

Best Linear Approximation Framework
The best linear approximation (BLA) of a nonlinear system is a modelling approach that minimizes the mean square error between the measured output of a nonlinear system and the output of the linear model [13]. It makes use of the knowledge that the excitation signal has both stochastic and deterministic properties. The excitation signal is a random phase multisine signal with M different (random phase) realizations; each of the blocks (periods) is repeated P times. By using multiple realizations (i.e., random phase rotations) of the multisines, the richness of the signal is increased. Instead of directly using the averaged input and output data (i.e., the use of classical H1 framework [12,16]), a partial BLA estimate is calculated for each period of the excitation [11]. A BLA FRM estimate (see Figure 1), for a given signal, is then calculated via the average of partial BLA estimates (see Figure 2). In this case we can easily estimate the noise and nonlinearity levels. Figure 1 shows the theoretical structure of the considered BLA estimator. In this figure, G linear is the linear (transfer function), phase coherent component of the model for which phase rotations at the input result in a proportional phase rotation at the output [11]. For the cases in which the output has phase non-coherent behavior, we can capture the (non-coherent) nonlinearities represented by G s . Increasing the number of random realizations of the multisine signal allows us to tackle the random output phase rotations as a "half-stochastic" (nonlinear) noise source. The "half-stochastic" term refers to the fact that G s does not vary over the repetitions of the same signal segment but only over different realizations. G E represents the classical FRF measurement noise. The usage of periodic excitation reduces the effects of the measurement noise component G E (see Equation (3)). The usage of multiple realizations reduces the impact of non-coherent nonlinearities represented by G s . G Bias represents the nonlinearities remaining after multiple realizations of the excitation signal.

Best Linear Approximation Framework
The best linear approximation (BLA) of a nonlinear system is a modelling approach that minimizes the mean square error between the measured output of a nonlinear system and the output of the linear model [13]. It makes use of the knowledge that the excitation signal has both stochastic and deterministic properties. The excitation signal is a random phase multisine signal with M different (random phase) realizations; each of the blocks (periods) is repeated P times. By using multiple realizations (i.e., random phase rotations) of the multisines, the richness of the signal is increased. Instead of directly using the averaged input and output data (i.e., the use of classical H1 framework [12,16]), a partial BLA estimate is calculated for each period of the excitation [11]. A BLA FRM estimate (see Figure 1), for a given signal, is then calculated via the average of partial BLA estimates (see Figure 2). In this case we can easily estimate the noise and nonlinearity levels. Figure 1 shows the theoretical structure of the considered BLA estimator. In this figure, is the linear (transfer function), phase coherent component of the model for which phase rotations at the input result in a proportional phase rotation at the output [11]. For the cases in which the output has phase non-coherent behavior, we can capture the (non-coherent) nonlinearities represented by . Increasing the number of random realizations of the multisine signal allows us to tackle the random output phase rotations as a "half-stochastic" (nonlinear) noise source. The "half-stochastic" term refers to the fact that does not vary over the repetitions of the same signal segment but only over different realizations.
represents the classical FRF measurement noise. The usage of periodic excitation reduces the effects of the measurement noise component (see Equation (3)). The usage of multiple realizations reduces the impact of non-coherent nonlinearities represented by . represents the nonlinearities remaining after multiple realizations of the excitation signal. The considered steady-state model at period and realization at excited frequency bins is given as a straightforward extension of Figure 1 If there is a transient term present in the measurement, it is either discarded or estimated [16]. For computational optimality, the FRM estimation is solved output-channel-wise (i.e., one output channel per unit time). The BLA estimate is calculated at the excited frequency lines only. In the frequency band of excitation, the BLA estimates at the non-excited bins are obtained by (linear) interpolation.
A partial FRM estimate is obtained via averaging over the periods (blocks) in a realization; see Equation (3). If is sufficiently large, then (considering the law of large numbers and the distribution properties of the measurement noise) the term [ ] in Equation (3) converges to zero.
Because the stochastic nonlinear contribution [ ] does not vary over the repetitions of the same realization, one has to average over the different realizations. If is sufficiently large, then the The considered steady-state model at period p and realization m at excited frequency bins is given as a straightforward extension of Figure 1 by: If there is a transient term present in the measurement, it is either discarded or estimated [16]. For computational optimality, the FRM estimation is solved output-channel-wise (i.e., one output channel per unit time). The BLA estimate is calculated at the excited frequency lines only.
In the frequency band of excitation, the BLA estimates at the non-excited bins are obtained by (linear) interpolation.
A partial FRM estimate is obtained via averaging over the P periods (blocks) in a realization; see Equation (3). If P is sufficiently large, then (considering the law of large numbers and the distribution properties of the measurement noise) the termĜ S does not vary over the P repetitions of the same realization, one has to average over the M different realizations. If M is sufficiently large, then the nonlinear noise sourceĜ s in Equation (3) converges to zero. After the 2D averaging (i.e., averaging over P and M) the BLA estimate is obtained. The method is illustrated in Figure 2. The noise covarianceσ 2Ĝ is estimated from the improved averaged sample varianceσ 2Ĝ [m] of each FRM realization (|X| .2 represents elementwise squaring and absolute value of the elements of matrix X): where inσ 2Ĝ [m] the additional normalization with P is needed to show the improved covariance (noise) estimate (this term corresponds to G E /P). In other words, averaging over repeated blocks results in an improvement of the SNR. Similarly, the additional normalization with M inσ 2Ĝ is needed to show the noise estimate improvement over different realizations (this term corresponds to G E /MP). If the user wants to see the covariance (noise) with respect to one period (block), one has to multiplyσ 2Ĝ with MP (this normalization is used in Figure 1). Vibration 2020, 3 FOR PEER REVIEW 4 nonlinear noise source in Equation (3) converges to zero. After the 2D averaging (i.e., averaging over and ) the BLA estimate is obtained. The method is illustrated in Figure 2. The noise covariance is estimated from the improved averaged sample variance [ ] of each FRM realization (| | . represents elementwise squaring and absolute value of the elements of matrix ): where in [ ] the additional normalization with is needed to show the improved covariance (noise) estimate (this term corresponds to / ). In other words, averaging over repeated blocks results in an improvement of the SNR. Similarly, the additional normalization with in is needed to show the noise estimate improvement over different realizations (this term corresponds to / ). If the user wants to see the covariance (noise) with respect to one period (block), one has to multiply with (this normalization is used in Figure 1). and stand for the improved (experiment-wise) stochastic nonlinearity and noise estimates. and stand for the period-wise stochastic nonlinearity and noise estimates (see Figure 1).
The total variance of the FRM is estimated from the improved variance (i.e., with extra normalization factor ) of each partial BLA estimate [ ] : The difference between the total variance and the noise variance is an estimate of the variance of the stochastic nonlinear contributions such that = − . If the user wants to see nonlinear contributions with respect to one period (block), he/she has to multiply with (this normalization mode is used in Figure 1). For the interpretation of different normalization modes, see Sections 4.1.1 and 4.1.2. For the detailed calculation, we refer to [11,14].

Introduction
Simplified analysis for multiple input systems (SAMI) is a user-friendly, Matlab-based toolbox, at its present state available for academic usage per request. It supports command line and graphical [m] are obtained via the period-wise estimates. The BLA estimateĜ BLA and its improved variance estimatesσ 2Ĝ BLA are obtained via partial BLA estimates (vertical direction).σ 2 NL andσ 2Ĝ stand for the improved (experiment-wise) stochastic nonlinearity and noise estimates. G S and G E stand for the period-wise stochastic nonlinearity and noise estimates (see Figure 1).

The total variance of the FRMσ 2Ĝ
BLA is estimated from the improved variance (i.e., with extra normalization factor M) of each partial BLA estimateĜ [m] : The difference between the total variance and the noise variance is an estimate of the variance of the stochastic nonlinear contributions such thatσ 2 NL =σ 2Ĝ BLA −σ 2Ĝ . If the user wants to see nonlinear contributions with respect to one period (block), he/she has to multiplyσ 2 NL with M (this normalization mode is used in Figure 1). For the interpretation of different normalization modes, see Sections 4.5.1 and 4.5.2. For the detailed calculation, we refer to [11,14].

Introduction
Simplified analysis for multiple input systems (SAMI) is a user-friendly, Matlab-based toolbox, at its present state available for academic usage per request. It supports command line and graphical user interfaces. The toolbox has been tested with many simulations and real-life industrial measurements, and it is optimized for MIMO experiments. The key idea is to make use of the statistical properties of the proposed excitation signal such that it becomes possible to split up the classical coherence function into noise, nonlinearity (and transient) components.
The toolbox provides a user-friendly interpretation of the nonlinear MIMO measurement data by extracting the user-relevant information and it addresses the questions related to the automatic processing of MIMO measurements with respect to the whole process. The warning-notification system was developed by experiences learned from typical mistakes which popped up during the measurement campaigns. The internal software architecture of SAMI consists of the following interconnected layers: • The design of experiment layer addresses the signal design, and the choice of measurement parameters.

•
The pre-processing layer considers a check-up of the input (reference) channels and provides an early warning to the user when the inputs are too strongly correlated. Furthermore, it segments the measurement data over the periods and realizations, and it sets up the processing parameters for the BLA estimation.

•
The BLA estimation layer provides the BLA FRF estimation, calculates advanced statistics of the noise and nonlinearity, and when it applies, estimates the transient term.

•
The post-processing layer makes the estimation results and warnings accessible in a condensed form. It provides users with the FRF, noise, and nonlinearity estimates. It is possible to automatically highlight the FRF (or input, output, reference) channels that have significant nonlinearity or noise levels based on the user-defined profiles. Furthermore, channels with sensory faults and/or imperfections, and correlated inputs are detected as well.
Vibration 2020, 3 FOR PEER REVIEW 5 processing of MIMO measurements with respect to the whole process. The warning-notification system was developed by experiences learned from typical mistakes which popped up during the measurement campaigns. The internal software architecture of SAMI consists of the following interconnected layers: • The design of experiment layer addresses the signal design, and the choice of measurement parameters.

•
The pre-processing layer considers a check-up of the input (reference) channels and provides an early warning to the user when the inputs are too strongly correlated. Furthermore, it segments the measurement data over the periods and realizations, and it sets up the processing parameters for the BLA estimation.

•
The BLA estimation layer provides the BLA FRF estimation, calculates advanced statistics of the noise and nonlinearity, and when it applies, estimates the transient term.

•
The post-processing layer makes the estimation results and warnings accessible in a condensed form. It provides users with the FRF, noise, and nonlinearity estimates. It is possible to automatically highlight the FRF (or input, output, reference) channels that have significant nonlinearity or noise levels based on the user-defined profiles. Furthermore, channels with sensory faults and/or imperfections, and correlated inputs are detected as well.
The communication between internal layers is via the all-in-one format of the toolbox: a single structure containing all the measurement, processing and estimation information. Figure 3 shows an overview of the SAMI's GUI. The displayed content is dynamic; under normal circumstances there is minimal user-interaction required (the details of the measurement and processing parameters are hidden); however, the user can overwrite any of the parameters.

Multisine Signals
The toolbox uses multisine signals to assess the underlying systems in a time-efficient way [13]. In order to avoid any spectral leakage, to reach full nonparametric characterization of the noise, and to be able to detect nonlinearities, a periodic signal is needed. Many users prefer noise excitations, The communication between internal layers is via the all-in-one format of the toolbox: a single structure containing all the measurement, processing and estimation information. Figure 3 shows an overview of the SAMI's GUI. The displayed content is dynamic; under normal circumstances there is minimal user-interaction required (the details of the measurement and processing parameters are hidden); however, the user can overwrite any of the parameters.

Multisine Signals
The toolbox uses multisine signals to assess the underlying systems in a time-efficient way [13]. In order to avoid any spectral leakage, to reach full nonparametric characterization of the noise, and to be able to detect nonlinearities, a periodic signal is needed. Many users prefer noise excitations, because they seem easier to implement, but in this case, nonlinearities are not directly identifiable, and there is a possible leakage error. The best signal that satisfies the desired properties is the easy to implement and apply a (random phase) multisine signal which looks like white noise and behaves like it, but it is not noise. If the multisine contains all/only odd or even harmonics, then it is called a full (band)/odd or even multisine. The usage of an odd multisine with randomly missing harmonics allows one to distinguish between even and odd nonlinear distortions on the non-excited bins (called the detection lines).
Basically, there are three modes supported in the toolbox: (1) combined (default, multiple random realizations and periods of missing odd); (2) fast (one realization, multiple periods of missing odd); and (3) robust (multiple realization of full) multisines. For a detailed overview of excitation signals we refer to [10,12]. The multisine module of SAMI with a signal illustration is shown in Figure 4.
Vibration 2020, 3 FOR PEER REVIEW 6 and there is a possible leakage error. The best signal that satisfies the desired properties is the easy to implement and apply a (random phase) multisine signal which looks like white noise and behaves like it, but it is not noise. If the multisine contains all/only odd or even harmonics, then it is called a full (band)/odd or even multisine. The usage of an odd multisine with randomly missing harmonics allows one to distinguish between even and odd nonlinear distortions on the non-excited bins (called the detection lines).
Basically, there are three modes supported in the toolbox: (1) combined (default, multiple random realizations and periods of missing odd); (2) fast (one realization, multiple periods of missing odd); and (3) robust (multiple realization of full) multisines. For a detailed overview of excitation signals we refer to [10,12]. The multisine module of SAMI with a signal illustration is shown in Figure  4.

Measurement Processing
Once the measurement is available, it can be easily processed with the toolbox. In most cases, automated data segmentation is used to retrieve the blocks (i.e., the individual periods and realizations). If meta-information is available, then the automated segmentation parameters are compared with the meta-information (or with the manually given information). If there is a significant difference between these, a warning is given. The automated warning system is based on a user-profile and can be changed any time. The toolbox works best with the generated signals; however, other type of excitation signals (noise, sine sweep, etc.) can be processed as well. Once the measurements are processed, a new window appears with the detailed analysis of the measurement. An example for a (MIMO) measurement of the F-16 aircraft is shown in Figure 5. The SISO F16 benchmark data is elaborated on in the next section. Figure 6 shows another interesting case with which to highlight the capabilities of the toolbox. A lowly damped car structure has been tested using 100 blocks of MIMO multisines; each block is 2 s long, the transient is 6 s (i.e., 3 blocks) long. The figure illustrates the situation which occurs if only one measured block is available: the classical H1 [16], windowed H1 [18,19], and the automated choice (in this case a special implementation of local rational method [15]) are compared with each

Measurement Processing
Once the measurement is available, it can be easily processed with the toolbox. In most cases, automated data segmentation is used to retrieve the blocks (i.e., the individual periods and realizations). If meta-information is available, then the automated segmentation parameters are compared with the meta-information (or with the manually given information). If there is a significant difference between these, a warning is given. The automated warning system is based on a user-profile and can be changed any time. The toolbox works best with the generated signals; however, other type of excitation signals (noise, sine sweep, etc.) can be processed as well. Once the measurements are processed, a new window appears with the detailed analysis of the measurement. An example for a (MIMO) measurement of the F-16 aircraft is shown in Figure 5. The SISO F16 benchmark data is elaborated on in the next section. Figure 6 shows another interesting case with which to highlight the capabilities of the toolbox. A lowly damped car structure has been tested using 100 blocks of MIMO multisines; each block is 2 s long, the transient is 6 s (i.e., 3 blocks) long. The figure illustrates the situation which occurs if only one measured block is available: the classical H1 [16], windowed H1 [18,19], and the automated choice (in this case a special implementation of local rational method [15]) are compared with each other.

F16 Measurement
This section concerns the ground vibration testing measurement campaign of a decommissioned F-16 aircraft with two dummy payloads mounted at the wing tips; see Figure 7. The detailed description of the measurement and benchmark data are openly accessible [20].
The right wing is excited by a shaker using combined multisines: odd multisines with skipping one random bin within each group of four successively excited odd lines. This sparse grid is used to detect even and odd nonlinear contributions. The sampling frequency is 400 Hz, and the period length 16,384 samples. The reference (voltage), input (force), and output (accelerations) signals are measured. The range of excitation is between 1 and 60 Hz; there are three periods and nine realizations per excitation level. There are three different input levels measured at 12.2, 49.0, and 97.1

F16 Measurement
This section concerns the ground vibration testing measurement campaign of a decommissioned F-16 aircraft with two dummy payloads mounted at the wing tips; see Figure 7. The detailed description of the measurement and benchmark data are openly accessible [20].
The right wing is excited by a shaker using combined multisines: odd multisines with skipping one random bin within each group of four successively excited odd lines. This sparse grid is used to detect even and odd nonlinear contributions. The sampling frequency is 400 Hz, and the period length 16,384 samples. The reference (voltage), input (force), and output (accelerations) signals are measured. The range of excitation is between 1 and 60 Hz; there are three periods and nine realizations per excitation level. There are three different input levels measured at 12.2, 49.0, and 97.1 N RMS.

F16 Measurement
This section concerns the ground vibration testing measurement campaign of a decommissioned F-16 aircraft with two dummy payloads mounted at the wing tips; see Figure 7. The detailed description of the measurement and benchmark data are openly accessible [20].
The right wing is excited by a shaker using combined multisines: odd multisines with skipping one random bin within each group of four successively excited odd lines. This sparse grid is used to detect even and odd nonlinear contributions. The sampling frequency is 400 Hz, and the period length 16,384 samples. The reference (voltage), input (force), and output (accelerations) signals are measured. The range of excitation is between 1 and 60 Hz; there are three periods and nine realizations per excitation level. There are three different input levels measured at 12.2, 49.0, and 97.1 N RMS.
In the 1-15 Hz band, the aircraft possesses about 10 resonance modes. The first few modes below 5 Hz correspond to rigid body motions of the structure. The first flexible mode around 5.2 Hz corresponds to wing bending deformations. The mode involving the most substantial nonlinear distortions is the wing torsion mode located around 7.3 Hz: the mounting interface of the payload features nonlinearities in stiffness and damping, due to clearance and friction. Therefore, the analysis is shown for the 1-15 Hz domain at the payload connection.
Vibration 2020, 3 FOR PEER REVIEW 8 corresponds to wing bending deformations. The mode involving the most substantial nonlinear distortions is the wing torsion mode located around 7.3 Hz: the mounting interface of the payload features nonlinearities in stiffness and damping, due to clearance and friction. Therefore, the analysis is shown for the 1-15 Hz domain at the payload connection.

Data Processing
The data processing is fully automated by the toolbox. In short, the first step is the automated segmentation of data: retrieving the periods and realizations. Next, the trends (such as the mean/offset values) from the individual segments are removed [13] and the transient is analyzed. Last, an early quality assessment of experimental data is done, and the parameters for the estimation procedure are set. Based on the SNRs of the input-output data, the existence of a reference signal, the number of periods, and the length of the transient, an appropriate estimator is automatically set.
Currently, the automated possibilities include BLA, H1, indirect BLA [13], and special implementations of the LPM (local polynomial method [10]), LRM (local rational method [15]), indirect LPM, and LRM. Note that with manual overriding more estimation methods are available. This section contains the results of the following single-line Matlab code directly used on the benchmark data: CreateAnalyzePlotAIO(Force(:), Acceleration(2,:),Voltage(:),Fs).
As can be seen, the toolbox requires only the signals in vector/matrix form, and optionally the sampling frequency (if it is not set, then a normalized frequency scale is used). The figures shown in this paper were obtained from the toolbox (for accessibility resized, with no warnings). Figure 8 shows the visualization of toolbox transient check-up routine. The left side of the figure shows the acceleration (output) measurement at the payload connection. In order to determine the length of the transient (i.e., the number of delay blocks), the last block (period)-assumed to be nearly in steady-state is subtracted from every preceding block. Because the transient decays as an exponential function, the differences are shown in logarithmic scale. Using automated statistical analysis of all available signals, the transient term is estimated as one block. Because there are more than two transient-free periods available, the input-output signals are measured sufficiently good (see next subsection). Because the reference signal is available, the indirect BLA estimator is automatically selected [13].

Data Processing
The data processing is fully automated by the toolbox. In short, the first step is the automated segmentation of data: retrieving the periods and realizations. Next, the trends (such as the mean/offset values) from the individual segments are removed [13] and the transient is analyzed. Last, an early quality assessment of experimental data is done, and the parameters for the estimation procedure are set. Based on the SNRs of the input-output data, the existence of a reference signal, the number of periods, and the length of the transient, an appropriate estimator is automatically set.
Currently, the automated possibilities include BLA, H1, indirect BLA [13], and special implementations of the LPM (local polynomial method [10]), LRM (local rational method [15]), indirect LPM, and LRM. Note that with manual overriding more estimation methods are available. This section contains the results of the following single-line Matlab code directly used on the benchmark data: CreateAnalyzePlotAIO(Force(:), Acceleration(2,:),Voltage(:),Fs).
As can be seen, the toolbox requires only the signals in vector/matrix form, and optionally the sampling frequency (if it is not set, then a normalized frequency scale is used). The figures shown in this paper were obtained from the toolbox (for accessibility resized, with no warnings). Figure 8 shows the visualization of toolbox transient check-up routine. The left side of the figure shows the acceleration (output) measurement at the payload connection. In order to determine the length of the transient (i.e., the number of delay blocks), the last block (period)-assumed to be nearly in steady-state is subtracted from every preceding block. Because the transient decays as an exponential function, the differences are shown in logarithmic scale. Using automated statistical analysis of all available signals, the transient term is estimated as one block. Because there are more than two transient-free periods available, the input-output signals are measured sufficiently good (see next subsection). Because the reference signal is available, the indirect BLA estimator is automatically selected [13].

Reference and Input Signal
It is crucial to analyze the excitation system and to quantify its accuracy and linearity. For that reason, the reference signals (generated multisine signals) are measured as well. From a system identification point of view, the actuator (shaker) can be seen as a SISO system (with voltage signal as input, the force signal as output). The reference, input signals, and the BLA estimate of the shaker used in this experiment are shown in Figure 9, at the low, medium, and high level of excitation. The left figure shows the generated reference (voltage) signals and their noise estimates. This measurement has very good quality (SNR is greater than 75 dB); the even and odd nonlinear distortions are hidden in the noise. The excitation forces (see middle figure) are measured with approximately 50-80 dB SNR. It is interesting to point out that the highest input signal is 18 dB higher than the lowest one, but the SNR is decreased with around 8 dB. This indicates the presence of (weak) nonlinearities at the excitation system. Similarly, this information can be seen by looking at the even and odd distortions: the higher the excitation, the higher the nonlinearities with an odd dominance (even distortions usually pop up as noise in the measurement, odd distortions usually influence the shapes of the spectrum).
The BLA analysis of the shakers is shown in the right figure. As can be seen, nonlinear distortions are present which are most likely due to (1) imperfections of the shaker; (2) the interactions between the aircraft and the shaker. Overall, the measurement quality (the gap between the actuator transfer

Reference and Input Signal
It is crucial to analyze the excitation system and to quantify its accuracy and linearity. For that reason, the reference signals (generated multisine signals) are measured as well. From a system identification point of view, the actuator (shaker) can be seen as a SISO system (with voltage signal as input, the force signal as output). The reference, input signals, and the BLA estimate of the shaker used in this experiment are shown in Figure 9, at the low, medium, and high level of excitation. The left figure shows the generated reference (voltage) signals and their noise estimates. This measurement has very good quality (SNR is greater than 75 dB); the even and odd nonlinear distortions are hidden in the noise.

Reference and Input Signal
It is crucial to analyze the excitation system and to quantify its accuracy and linearity. For that reason, the reference signals (generated multisine signals) are measured as well. From a system identification point of view, the actuator (shaker) can be seen as a SISO system (with voltage signal as input, the force signal as output). The reference, input signals, and the BLA estimate of the shaker used in this experiment are shown in Figure 9, at the low, medium, and high level of excitation. The left figure shows the generated reference (voltage) signals and their noise estimates. This measurement has very good quality (SNR is greater than 75 dB); the even and odd nonlinear distortions are hidden in the noise. Figure 9. The measured reference (voltage) and input (force) signals, and the FRFs of the actuator (shaker), are shown. Darker shades refer to higher excitation levels. Thick grey shades refer to the signals, FRFs. Thin grey shades refer to noise estimates. Orange shades (on signal measurements) refer to the odd distortions. Blue shades (on signal measurement) refer to the even distortions. Red shades (on FRF estimations) refer to (aggregated) nonlinear distortions. Observe that the higher the excitation, the more the nonlinear distortion on the output and FRF measurements.
The excitation forces (see middle figure) are measured with approximately 50-80 dB SNR. It is interesting to point out that the highest input signal is 18 dB higher than the lowest one, but the SNR is decreased with around 8 dB. This indicates the presence of (weak) nonlinearities at the excitation system. Similarly, this information can be seen by looking at the even and odd distortions: the higher the excitation, the higher the nonlinearities with an odd dominance (even distortions usually pop up as noise in the measurement, odd distortions usually influence the shapes of the spectrum).
The BLA analysis of the shakers is shown in the right figure. As can be seen, nonlinear distortions are present which are most likely due to (1) imperfections of the shaker; (2) the interactions between the aircraft and the shaker. Overall, the measurement quality (the gap between the actuator transfer The excitation forces (see middle figure) are measured with approximately 50-80 dB SNR. It is interesting to point out that the highest input signal is 18 dB higher than the lowest one, but the SNR is decreased with around 8 dB. This indicates the presence of (weak) nonlinearities at the excitation system. Similarly, this information can be seen by looking at the even and odd distortions: the higher the excitation, the higher the nonlinearities with an odd dominance (even distortions usually pop up as noise in the measurement, odd distortions usually influence the shapes of the spectrum).
The BLA analysis of the shakers is shown in the right figure. As can be seen, nonlinear distortions are present which are most likely due to (1) imperfections of the shaker; (2) the interactions between the aircraft and the shaker. Overall, the measurement quality (the gap between the actuator transfer function and its nonlinear estimates) is good, much higher than 20 dB needed to fulfil the assumption on the precise excitation signal measurement.

Payload Measurement
Further, in order to simplify the analysis, the output and FRF are shown at the payload connection only. The output (acceleration) measurements are shown in Figure 10. As can be seen, the SNR is around 40-50 dB at the resonances. At the higher excitation level, the SNR decreased with approximately 10 dB with respect to the lowest level excitation.
Vibration 2020, 3 FOR PEER REVIEW 10 function and its nonlinear estimates) is good, much higher than 20 dB needed to fulfil the assumption on the precise excitation signal measurement.

Payload Measurement
Further, in order to simplify the analysis, the output and FRF are shown at the payload connection only. The output (acceleration) measurements are shown in Figure 10. As can be seen, the SNR is around 40-50 dB at the resonances. At the higher excitation level, the SNR decreased with approximately 10 dB with respect to the lowest level excitation. Figure 10. The output (acceleration) measurement shown at the payload connection. Darker shades refer to higher excitation levels. Thick grey shades refer to the signals. Thin grey shades refer to noise estimates. Orange shades (on signal measurements) refer to the odd distortions. Blue shades (on signal measurement) refer to the even distortions. Observe that the higher the excitation, the more the dominance of odd nonlinear distortions.
It can also be observed that the resonances were shifted, which is also a further indication of nonlinearities. This is due to the fact that at higher level of excitation we have dominant odd distortions, which usually manifest in changing resonance locations and shapes (it is the so-called hardening or softening stiffness nonlinearity effect).
Please note that even distortions usually manifest as excessive noise on the measurement (see Figure 9).
When the actuator is imperfect, the interpretation of the levels on the detection lines at the output can be jeopardized. If there is power on the (non-excited) detection lines at the input signal power, then the output signals will contain not only the nonlinearities due the system but also the nonlinear contributions of the actuator. However, the toolbox automatically cleans up the output spectrum with the impurities of the actuator [13]. Figure 11 shows the FRFs at low, medium, and high levels of excitation. This is the classical approach; FRFs at multiple levels of excitation are compared with each other. It is interesting to point out that despite the fact that the highest excitation is only 18 dB higher than the lowest level excitation, it can be clearly observed that FRFs at different levels differ a lot from each other by, for instance, shifting resonances and varying damping. This clearly indicates the presence of nonlinearities. It can also be observed that the resonances were shifted, which is also a further indication of nonlinearities. This is due to the fact that at higher level of excitation we have dominant odd distortions, which usually manifest in changing resonance locations and shapes (it is the so-called hardening or softening stiffness nonlinearity effect).

FRFs at Different Excitation Levels
Please note that even distortions usually manifest as excessive noise on the measurement (see Figure 9).
When the actuator is imperfect, the interpretation of the levels on the detection lines at the output can be jeopardized. If there is power on the (non-excited) detection lines at the input signal power, then the output signals will contain not only the nonlinearities due the system but also the nonlinear contributions of the actuator. However, the toolbox automatically cleans up the output spectrum with the impurities of the actuator [13]. Figure 11 shows the FRFs at low, medium, and high levels of excitation. This is the classical approach; FRFs at multiple levels of excitation are compared with each other. It is interesting to point out that despite the fact that the highest excitation is only 18 dB higher than the lowest level excitation, it can be clearly observed that FRFs at different levels differ a lot from each other by, for instance, shifting resonances and varying damping. This clearly indicates the presence of nonlinearities. However, for a given level of excitation, an inexperienced user is not able to determine if there are nonlinearities present using the classical approach. The usage of the proposed multisines allows us to obtain noise and nonlinearity level estimations, as explained in Section 2.2. With the help of these curves one can distinguish between the effects of noise and nonlinearity on the coherence. For instance, when looking at the high-level excitation case, the most dominant resonance (around 7.3 Hz) has an approximate SNR of 43 dB, and an SNLR (signal-to-nonlinearity ratio) of 15 dB. This means that at that resonance the main error source is the nonlinearity. If a linear model is used, then the expected error level will be in the order of the SNLR. If an appropriate nonlinear model is used, the expected error level could drop to the SNR. This kind of extra information would have been impossible to derive from the H1 framework.

FRFs at Different Excitation Levels
Furthermore, the toolbox provides an automated warning system for the FRF, actuator, reference, input, and output signals when significant nonlinearities, noise, sensory faults, and correlation issues are discovered. This is partly illustrated in Figure 5 Figure 11 above shows the noise and nonlinearity estimates normalized with respect to one block (period). This normalization mode is very important for understanding, modeling, simulation, and control. Of course, it is possible to show the improved noise and nonlinearity estimates with respect to the whole experiment (measurement). The toolbox supports four different normalization modes illustrated in Figure 12: block-wise (default), realization-wise, and experiment-wise normalization modes, and nonlinear detection mode. We authors recommend always checking the period-wise (block-wise) noise and nonlinearity levels first (mode 0). Using many realizations and simultaneously SNR SNLR possible improvement Figure 11. The FRF estimation shown at the payload-right wing connection. Darker shades refer to higher excitation levels. Thick grey shades refer to the FRFs. Thin grey shades refer to noise estimates. Red shades refer to the nonlinear distortions. Observe that the higher the excitation: (1) the more the nonlinear distortions on the FRF measurements; (2) the lower the resonances; and (3) the higher the damping. However, for a given level of excitation, an inexperienced user is not able to determine if there are nonlinearities present using the classical approach. The usage of the proposed multisines allows us to obtain noise and nonlinearity level estimations, as explained in Section 2.2. With the help of these curves one can distinguish between the effects of noise and nonlinearity on the coherence. For instance, when looking at the high-level excitation case, the most dominant resonance (around 7.3 Hz) has an approximate SNR of 43 dB, and an SNLR (signal-to-nonlinearity ratio) of 15 dB. This means that at that resonance the main error source is the nonlinearity. If a linear model is used, then the expected error level will be in the order of the SNLR. If an appropriate nonlinear model is used, the expected error level could drop to the SNR. This kind of extra information would have been impossible to derive from the H1 framework.
Furthermore, the toolbox provides an automated warning system for the FRF, actuator, reference, input, and output signals when significant nonlinearities, noise, sensory faults, and correlation issues are discovered. This is partly illustrated in Figure 5 Figure 11 above shows the noise and nonlinearity estimates normalized with respect to one block (period). This normalization mode is very important for understanding, modeling, simulation, and control. Of course, it is possible to show the improved noise and nonlinearity estimates with respect to the whole experiment (measurement). The toolbox supports four different normalization modes illustrated in Figure 12: block-wise (default), realization-wise, and experiment-wise normalization modes, and nonlinear detection mode. We authors recommend always checking the period-wise (block-wise) noise and nonlinearity levels first (mode 0). Using many realizations and simultaneously showing the improved nonlinearity level (experiment-wise normalization, mode 2) might give an incorrect message to the user, since it converges to zero. It is also possible to show the noise improvement over one realization while freezing the nonlinearity level (mode 1). The last normalization mode (nonlinear detection mode, mode 3) can be used to detect the presence of nonlinearities (nonlinear variability mode). The transition between different noise and nonlinear modes can be seen in Figure 2 as well.
Vibration 2020, 3 FOR PEER REVIEW 12 showing the improved nonlinearity level (experiment-wise normalization, mode 2) might give an incorrect message to the user, since it converges to zero. It is also possible to show the noise improvement over one realization while freezing the nonlinearity level (mode 1). The last normalization mode (nonlinear detection mode, mode 3) can be used to detect the presence of nonlinearities (nonlinear variability mode). The transition between different noise and nonlinear modes can be seen in Figure 2 as well. Figure 12. FRFs, noise, and nonlinearity estimates shown at different normalization modes at a high level of excitation. Mode 0/1/2/3 refers to the block-wise/realization-wise/experiment-wise/nonlinear detection mode. Observe that the BLA FRF remains the same estimate. The noise and nonlinear distortion estimates are shown in case of normalization mode 0 with respect to one block; mode 1 with respect to one realization; mode 2 with respect to the whole experiment (all realizations). Mode 3 represents the nonlinear (variability) detection mode (where the nonlinearities are shown with respect to one block; noise is shown with respect to the all blocks). ̂2 and ̂2 stand for the improved (experiment-wise) nonlinearity and noise estimates. and stand for the period-wise nonlinearity and noise estimates (see Figures 1 and 2).

Using One Realization
We authors recommend using as many realizations (and periods) as possible. However, it is sometimes not possible to measure for long. Using the proposed multisines, it is possible to have a rough estimate about the nonlinearity and noise levels using at least one realization and two periods. In the classical literature [13] the one-realization case is called the fast method; the multiple realizations case (with full multisines) is called the robust method, where robust refers to the robustness estimate of the nonlinear and noise quantities (see also Section 3.2).
A comparison of the BLA quantities between 1 and 9-realization cases are shown for low and high levels of excitation in Figure 13. Observe that at the low level of excitation (left figure), the noise and nonlinearity levels are almost identical and the FRF estimate is smooth, whereas at high level of excitation, the one realization FRF looks noisy; its noise estimates are higher. This is why many realizations (and periods) are advised. Figure 12. FRFs, noise, and nonlinearity estimates shown at different normalization modes at a high level of excitation. Mode 0/1/2/3 refers to the block-wise/realization-wise/experiment-wise/nonlinear detection mode. Observe that the BLA FRF G BLA remains the same estimate. The noise and nonlinear distortion estimates are shown in case of normalization mode 0 with respect to one block; mode 1 with respect to one realization; mode 2 with respect to the whole experiment (all realizations). Mode 3 represents the nonlinear (variability) detection mode (where the nonlinearities are shown with respect to one block; noise is shown with respect to the all blocks).σ 2 NL andσ 2Ĝ stand for the improved (experiment-wise) nonlinearity and noise estimates. G S and G E stand for the period-wise nonlinearity and noise estimates (see Figures 1 and 2).

Coherence Function
Industrial practitioners tend to look at the (multiple) coherence function estimation [21] only, since the interpretation of the variances requires hands-on experience. When the coherence function gives one, it means that there is a 100% linear correlation between the measured output and input data. When the coherence is lower than one, it indicates the presence of (among others) high-level noise and/or transient and/or leakage and/or nonlinearities. Using the toolbox, the multiple coherence function estimate can be split up to noise, nonlinearity, and transient components. Figure 14 shows the detailed coherence information. The analysis of this figure is very exciting. The coherence estimate with respect to the total measurement (thick black line) shows how many "imperfections" can be found in the measurement. When this coherence is split into components, one can see that the largest lack of coherence originates from the transient data (pink line). When an appropriate number of delay blocks are cut out from the measurement, this component is eliminated.
The second largest lack of coherence contributor is the nonlinear component (red line). When an adequate nonlinear model is used, this component will be eliminated. The last component is the noise coherence component (thin black lines with asterisks), which tells us that the measurement quality was excellent. The advantage of the variance technique over the coherence component technique is that different normalizations can be used; relative information is shown. For each pair of input-output, different variance levels are assigned. This restriction of the coherence component technique is related to the processing method (which is evaluated measurement-wise, experiment-wise).

Coherence Function
Industrial practitioners tend to look at the (multiple) coherence function estimation [21] only, since the interpretation of the variances requires hands-on experience. When the coherence function gives one, it means that there is a 100% linear correlation between the measured output and input data. When the coherence is lower than one, it indicates the presence of (among others) high-level noise and/or transient and/or leakage and/or nonlinearities. Using the toolbox, the multiple coherence function estimate can be split up to noise, nonlinearity, and transient components. Figure 14 shows the detailed coherence information. The analysis of this figure is very exciting. The coherence estimate with respect to the total measurement (thick black line) shows how many "imperfections" can be found in the measurement. When this coherence is split into components, one can see that the largest lack of coherence originates from the transient data (pink line). When an appropriate number of delay blocks are cut out from the measurement, this component is eliminated.
The second largest lack of coherence contributor is the nonlinear component (red line). When an adequate nonlinear model is used, this component will be eliminated. The last component is the noise coherence component (thin black lines with asterisks), which tells us that the measurement quality was excellent.
Vibration 2020, 3 FOR PEER REVIEW 13 Figure 13. FRF, noise, and nonlinearity estimates are shown at low and high level of excitation using only one and all realizations.

Coherence Function
Industrial practitioners tend to look at the (multiple) coherence function estimation [21] only, since the interpretation of the variances requires hands-on experience. When the coherence function gives one, it means that there is a 100% linear correlation between the measured output and input data. When the coherence is lower than one, it indicates the presence of (among others) high-level noise and/or transient and/or leakage and/or nonlinearities. Using the toolbox, the multiple coherence function estimate can be split up to noise, nonlinearity, and transient components. Figure 14 shows the detailed coherence information. The analysis of this figure is very exciting. The coherence estimate with respect to the total measurement (thick black line) shows how many "imperfections" can be found in the measurement. When this coherence is split into components, one can see that the largest lack of coherence originates from the transient data (pink line). When an appropriate number of delay blocks are cut out from the measurement, this component is eliminated.
The second largest lack of coherence contributor is the nonlinear component (red line). When an adequate nonlinear model is used, this component will be eliminated. The last component is the noise coherence component (thin black lines with asterisks), which tells us that the measurement quality was excellent. The advantage of the variance technique over the coherence component technique is that different normalizations can be used; relative information is shown. For each pair of input-output, different variance levels are assigned. This restriction of the coherence component technique is related to the processing method (which is evaluated measurement-wise, experiment-wise). The advantage of the variance technique over the coherence component technique is that different normalizations can be used; relative information is shown. For each pair of input-output, different variance levels are assigned. This restriction of the coherence component technique is related to the processing method (which is evaluated measurement-wise, experiment-wise).

Conclusions
In this work a novel automated MIMO BLA framework was developed to provide a user-friendly interpretation of the nonlinear behavior of MIMO measurement data by extracting user-relevant information. The toolbox turned out to be useful for modelling nonlinear FRFs because: • It requires minimal user-interaction, and no expert-user • Excitation signals are optimized for structures with multiple inputs • The reference, input and output measurements were nonparametrically characterized • The actuator system is characterized to improve the output measurement quality in case of nonlinear actuator and/or feedback • An advanced frequency response matrix estimation method is used. It is a simple but robust estimation process; at each excitation level of noise, nonlinearity and transient information can be retrieved.