Synthesis of Shaker Input Signals with High Kurtosis for Vibration Qualification Testing

1. Introduction

In the operating lifetime of mechanical systems working at high dynamics, vibrations are a common source of fatigue damage that could lead components to premature failure. To test whether the components could withstand the elastodynamic loads occurring for a particular application, just the latter is used to tailor qualification tests performed by means of shakers/slip tables. The so -<called Test Tailoring procedure aims at a proper definition of the vibratory profiles (Mission Synthesis) to be used as the excitation for the device under test (DUT), based on the processing of reference environmental data – measured in those working conditions that are particularly significant to represent real-life excitations expected for the DUT – in order to reproduce their most important characteristics in laboratory tests.

Instead of replicating the recorded environmental data as shaker input signals, which would result in the major drawback of losing stochasticity and limiting the qualification tests to represent a given working condition only [1], the conventional Mission Synthesis procedures implemented so far [2] provide the synthesis of a test profile in terms of a Power Spectral Density (PSD). In Random Control tests, the shaker controller then generates the actual vibration after applying the Inverse Fast Fourier Transform (IFFT) to the spectral domain. Amplitudes and phases characterize this domain, with the former being obtained from the PSD and the latter being randomly generated as uniformly distributed random variables.

Owing to their random generation, the probability distribution of the values of the obtained time-series is Gaussian. This could compromise the reliability of tests since most excitation signals encountered in real applications feature non-Gaussian distributions. A parameter that accounts for deviations from Gaussianity is kurtosis [3]: for Gaussian signals, its theoretical value amounts to 3.0, whereas if peaks and bursts or deterministic components are preponderant, it becomes higher or lower than 3.0, respectively. Signals with kurtosis value greater than 3.0 are named Leptokurtic, and so is their distribution.

Over the last decades, the so-called kurtosis control methods have gained increasing interest: the aim is to generate shaker input signals by controlling both the PSD and the kurtosis parameter, in order to carry out more realistic tests by preserving the nature of real-life random excitations (i.e. the values probability distribution due to peaks and bursts). Several kurtosis-control algorithms can be found in the literature [1, 4-17], which are essentially based on three different main approaches (or their possible combinations): Phase Manipulation [5, 8, 12, 16, 17], Amplitude Modulation [10, 11], and Polynomial Transformation [13-15]. The methods proposed feature different characteristics and performance depending on the corresponding targeted goals (not limited to durability testing applications, which is instead the general framework of the present work). Overall, however, one common critical issue exists that needs to be addressed by any procedure of kurtosis-control, as it significantly affects the effectiveness of transferring peaks and bursts from the excited base to the DUT. In fact, it is known that the response of a lightly damped system to a generic excitation may tend to a Gaussian probability distribution due to a filtering effect [18]. Thus, if this occurs, the tests would be no different from the standard procedure that generates directly Gaussian signals. Hereinafter, this filtering effect will be referred to as Papoulis’ Rule, which stems from the mathematician who proved it based on the Central Limit Theorem [18, 19].

In this work, two novel algorithms are proposed for the synthesis of random excitation signals with prescribed PSD and kurtosis value, able to robustly circumvent the Papoulis’ Rule and to excite a generic DUT effectively irrespective of the signal features of the reference environmental data. Indeed, most kurtosis control algorithms prove to be effective only if certain conditions are met (e.g.: the DUT natural frequencies must be known; the reference environmental data to be processed need to be narrowband signals…). The main difference between the two algorithms regards the predicted kurtosis value of the DUT response: constant and equal to the excitation one for the first algorithm, variable and significantly larger for the second one. Moreover, it will also be shown that the proposed methods can be used for accelerated fatigue-life tests with the help of the filter proposed by Kihm et al. [20]. These tests seek to preserve the damage potential of a signal measured from an application and replicate it in a shorter time. The damage potential associated with a vibratory excitation is estimated via a spectral function called Fatigue Damage Spectrum (FDS) [21] and current procedures already available in some commercial software permit the synthesis of a PSD from prescribed FDS and duration of the test [2]. Since the output is the only PSD, the standard procedure leads to accelerated tests featuring Gaussian excitation signals, hence the same problems with not preserving the nature of reference signals – possibly non-Gaussian – may occur in this case as well. The two kurtosis-control algorithms do not directly control the FDS, but by convolving the synthesized signal with the filter mentioned above, it is possible to match a prescribed FDS for accelerated fatigue-life tests while preserving the non-Gaussian features of the environmental data.

The paper structure is as follows: Section 2 reports the theoretical background on the popular procedures of kurtosis-control and the analytical formulations of the two original algorithms proposed in this work; the results of the application of the algorithms for a case study are reported and discussed in Section 3 and Section 4, respectively; finally, Section 5 is devoted to some concluding remarks. In addition, details about some computational issues related to the algorithms implementation are reported in the Appendices: in particular, computationally efficient procedures for the implementation of kurtosis-control algorithms based on Phase Modulation and for the kurtosis value computation of multiple/concatenated time-series are proposed in Appendix A and Appendix B, respectively, whereas the effectiveness of the special filter proposed by Kihm et al. in [20] is analytically proven in Appendix C.

2. Materials and Methods

2.1. Theoretical background

In the field of random vibration testing, the shaker controller typically generates a stationary vibratory signal in the time domain from a prescribed PSD input by applying the IFFT to the spectral domain. Amplitudes and phases characterize the spectral domain; the former are related to the PSD via the equation:

(1)

where the amplitude A_n of the n^th harmonic is related to the n^th PSD element G_n and frequency resolution Δf. The randomness required for the generation of the time-series is guaranteed by the phases, which are defined as uniformly distributed random variables in the interval [0, 2π[. The time histories generated in such a way are always characterized by a probability distribution of their values close to Gaussian. The statistical parameter known as kurtosis is the 4^th statistical central moment of a signal, normalized by the 4^th power of its standard deviation. For Gaussian signals, the theoretical value of kurtosis is 3.0, whereas Leptokurtic signals feature greater values (e.g. due to high peaks caused by micro-collisions). In practical applications, discretized formulations of the signal statistical quantities are typically used:

(2)

(3)

(4)

(5)

where $\overline{x}$, σ, σ², k are the mean value, the standard deviation, the variance, and the kurtosis of the signal, respectively, x_n its n^th sample, and N its total number of samples. Current kurtosis control algorithms are generally based on three approaches: Phase Manipulation (PM), Polynomial Transformation (PT), and Amplitude Modulation (AM).

PM methods [5, 8, 12, 16, 17] employ the formulation of kurtosis written in terms of the amplitudes and phases of a time-series. Steinwolf et al. [4-9, 17] described some analytical methods that control kurtosis by changing only the phases while keeping the amplitudes constant. This would imply preserving the PSD as well since it only depends on amplitudes. In Appendix A, computational cost issues are discussed and an original procedure is proposed to efficiently implement the PM algorithms based on Steinwolf’s studies.

PT methods consist of generating a Gaussian signal x(t) with the prescribed PSD first and then applying an analytical transformation of the form:

(6)

with the coefficients α₁ and α₃ being functions of the target kurtosis value [13-15]. The transformation in Eq. (6) presents two disadvantages: (i) the PSD of the signal is affected by an unwanted disturbance and (ii) Papoulis’ Rule is very likely to occur, thus causing the peaks present in the input signal to be filtered out and leading the response to approaching a Gaussian probability distribution.

The second problem is similar to the case of PM methods. In fact, the effectiveness of kurtosis-control methods lies in the fact that the bursts of the signal must have either of two characteristics: (i) a long enough duration in order to appear also in the unavoidably time-delayed response of the system, (ii) a narrow-band frequency content containing the natural frequency of the system. Typically, the bursts generated by both PM and PT methods do not have a long enough duration and this is the reason why they are likely to be filtered out by the system.

Methods that generate bursts of desired duration employ the AM approach, which consists in modulating a Gaussian signal x(t) having the desired PSD with an appropriate function w(t) in order to obtain a signal with a desired kurtosis value [10, 11]:

(7)

This method is effective in transferring the kurtosis value to the response of the DUT if the signal bursts of the modulating signal have greater duration than the inverse of the bandwidth of the lightly-damped system [11]. The carrier waveform w(t) introduces low frequency components in the spectrum of y(t) compared to that of x(t), albeit negligible if w(t) is appropriately designed.

Both the algorithms presented in this paper, which aim to synthesize vibratory profiles whose kurtosis and PSD match the reference input ones, are based on two different approaches that can hardly be classified into the above-mentioned categories. The following Sections are devoted to the detailed description of their implementation.

2.2. Multi-Level Variance (MLV) algorithm

The algorithm named Multi-Level Variance (MLV) attains the synthesis of a signal by dividing the signal duration T_tot into n_b blocks of the same size and duration T_b (T_b = T_tot/n_b). The generated blocks are not overlapped and have different variance, which is closely related to a modulation procedure, although no modulating function is explicitly used. As shown in the following, the different levels of variance σ_i² (i = 1,…, n_b) are produced in such a way that the synthesized signal complies with the kurtosis and PSD constraints. In general, the PSD of a signal is computed by calculating the Fast Fourier Transforms (FFT) over small-sized blocks and squaring their magnitude to obtain the so-called periodograms. More specifically, the periodogram could be thought of as some sort of PSD computed only for the generic block of the signal. After obtaining the periodograms, the last step is to average them in order to calculate the PSD of the signal. In the algorithm, the PSD of the i^th block of the signal, G_i, is defined to be proportional to the PSD of the reference signal, G, as in the following expression:

(8)

In Eq. (8) σ_tot² and σ_i² are the variance values of the overall signal and the i^th block, respectively. It is to be highlighted that the σ_i² parameters are the unknowns, whereas σ_tot² can be calculated from either the reference signal, Eq. (4), or directly from the reference PSD:

(9)

where G(f) is the PSD amplitude at the generic frequency f of a continuous signal. However, since the processed signal is discrete in practice, the PSD is also discrete and the theoretical computation of Eq. (9) must be discretized. The unknown parameters σ_i² are also related to the overall variance σ_tot² through the following equation (cf. Appendix B for proof):

(10)

From Eqs. (8, 10) the following relation must hold:

(11)

Since the PSD G is computed by averaging the PSDs of the blocks, Eq. (11) is automatically satisfied. Hence, the constraint on the PSD spectrum is fulfilled if n_b values σ_i² that comply with Eq. (10) are found. In addition to Eqs. (10, 11), there is also a relation between the kurtosis values of the overall signal, k_tot, and of the single blocks, k_i (cf. Appendix B):

(12)

The first step of the algorithm is to randomly generate the σ_i² values such that σ_i² ∈ [σ_min², σ_max²] and σ_min² < σ_tot² < σ_max². The ratio r_σ² = σ_min²/σ_tot² ∈ ]0, 1[ must be selected by the user as an input for the algorithm: the closer to 0, the more the variance of the synthesized signal will vary in time, whereas the opposite is true if closer to 1. The parameter σ_max² is not selectable directly by the user because it is adjusted throughout the algorithm iterations in order to converge to the target kurtosis value. The user can also choose the number of bursts with high amplitude excursion of the synthesized signal (parameter n_p = 0,…, n_b); this is typically done by referring to the number of bursts that can be observed in the reference measured data. The algorithm generates n_p blocks with a variance equal to σ_max², which is greater than the variance of the other blocks. The kurtosis k_i and variance σ_i² of the blocks are calculated via Eqs. (4, 5) only on the first iteration, whereas the discrete signal x_n is obtained by the IFFT of the PSD of the single blocks with randomly generated phases and then by concatenating the blocks. A random integer s, representing a single block, is then automatically generated in the interval [1, n_b]. After the σ_i² levels are generated (i = 1,…, n_b, i ≠ s), the variance and the kurtosis of the randomly selected s^th block are calculated via the following relations, which stem from Eqs. (9, 12):

(13)

(14)

Eqs. (13, 14) are used in order to verify if the prescribed PSD and target kurtosis can be achieved: indeed, σ_s² is required to be not null and k_s greater than a lower threshold, k_{s_min}, and smaller than an upper threshold, k_{s_max}. The upper threshold should not be set excessively high because it would lead to generating unrealistic peaks in the s^th block exceeding by far the amplitude of those of the other blocks. Afterward, in order to obtain the k_s given by Eq. (14), the harmonics phases of the s^th block are adjusted using a phase manipulation procedure (cf. Section 2). With very few iterations, where the parameter σ_max² is changed in order to converge towards the target kurtosis value, Eqs. (13, 14) are usually satisfied. The final step of the algorithm is to concatenate the generated blocks by smoothing them through proper interpolation of the values close to the edges of the adjoining blocks, in order to avoid unrealistic discontinuities among them.

In conclusion, the user has to insert:

• the reference input signal or, alternatively, reference PSD and kurtosis value;
• the duration T_tot of the signal to be synthesized;
• the sampling frequency F_s of the synthesized signal (usually the same as the reference signal);
• the signal blocks duration T_b;
• the ratio r_σ²;
• the number of bursts n_p;
• the lower and upper thresholds for kurtosis value k_{s_min} and k_{s_max}, respectively.

2.3. Variable Spectral Density (VSD) algorithm

The algorithm hereinafter called Variable Spectral Density (VSD) splits the signal to be synthesized into n_b blocks of the same duration T_b as the previous algorithm. The major difference is that the PSD G_i of the i^th block is randomly generated.

Let the PSD matrix [G’’_ij] be defined as:

(15)

where N_h is the number of spectral lines. This matrix has N_h rows and n_b columns, with the j^th column being the PSD vector of the j^th block of the reference signal. Eq. (15), where all the columns have the same elements (harmonic amplitudes), refers implicitly to a signal having a stationary PSD. On the other hand, the PSD matrix [G’’_ij] corresponding to the signal synthesized through the MLV algorithm has the following form:

(16)

Both the matrices in Eqs. (15) and (16) satisfy Eq. (11) that may be rewritten in this case, in conformity with the notation used in this Section, as:

(17)

The VSD algorithm synthesizes a signal with a PSD variable over time, corresponding to a PSD matrix having the most general form:

(18)

where the elements must comply with Eq. (17).

The generic matrix in Eq. (18) can be derived from Eqs. (15) and (17), by changing one row, the i^th for instance, as in the following equation, where p ∈ [0, 1] and l is a positive integer such that l ≤ n_b:

(19)

Equation (17) is still complied with if the terms of the type pG_i are l – 1. If similar operations were done not only on the term G_ij but also on other terms and in a random manner, then the PSD of each block could be varied still preserving the overall PSD.

In conclusion, the user is required to perform the following operations:

• insert a reference input signal or, alternatively, reference PSD and kurtosis value;
• set the duration T_tot and the sampling frequency F_s of the signal to be synthesized;
• choose p ∈ [0, 1].

Then, the workflow of the algorithm is based on the following operations that are automatically performed:

• set T_b and i = 1;
• set s = 0;
• choose a random element j of the i^th row;
• generate a positive random integer l ≤ n_b – s;
• set G_ij = [1 + (l – 1)(1 – p)G_i and s = s + l;
• repeat 3 – 4 – 5 with another value for j (different from the values generated in the previous loops) and another value of l, until s ≥ n_b – 1;
• set G_im = pG_i with m ranging over all the elements of the i^th row which have not been modified in point 5;
• if i < N_h, set i = i+1 and repeat from point 2, otherwise proceed to point 9;
• terminate if the kurtosis of the synthesized signal matches the target value (within a certain tolerance, to be preliminarily set), otherwise repeat from point 1 where a different T_b is automatically generated. Decreasing T_b makes the kurtosis value increase and vice versa (this is how the algorithm converges towards the target kurtosis).

As in the case of the MLV algorithm, the very last step is to smooth the signal at the edges of the blocks in order to avoid unrealistic discontinuities.

It is worth noting that a limitation of the algorithm could be that the variation of the PSD over time (i.e. over the blocks) is not controlled but randomly generated.

3. Results

Real field data sampled at 8192Hz and having a duration of 660s are used as reference signal for an example of application (Figure 1). The synthesis by the MLV and VSD algorithms was performed in order to obtain signals with the same sampling frequency and duration for the sake of a direct and more straightforward comparison. Two signals synthesized by the MLV algorithm are shown in Figure 2, hereinafter referred to as MLV1 and MLV2, respectively. The main difference between the two implementations of the algorithm was in the choice of the signal blocks durations and number of bursts: {T_b = 3s, n_p = 10} and {T_b = 0.1s, n_p = 300}, for MLV1 and MLV2 signals, respectively. The parameter r_σ² was set equal to 0.25 for both. Statistical parameters and PSDs of the signals are respectively reported in Table 1 and Figure 3, over the frequency range 0÷4096Hz. The acceleration response of a series of single-degree-of-freedom (SDOF) linear systems, with natural frequency in the range 0÷4000Hz and a constant damping ratio equal to 2.5%, was calculated by convolution between the impulse response of the system and its input. The response kurtosis of the synthesized and reference signals is shown in Figure 4, where the kurtosis parameter is plotted versus the natural frequencies of the SDOF systems. It can be observed that the kurtosis of the synthesized signals remains approximately constant over the frequency range; however, the signal with the lower T_b (MLV2) suffers more from the Papoulis’ Rule, whereas the higher T_b (MLV1) permits the generation of bursts that last enough to be transferred to the response, thus almost preserving the input kurtosis value. By inspecting the plot, it should be evident that the kurtosis of the response to the reference signal is not constant over the frequency range and in particular, for some frequencies it is much higher.

The VSD algorithm only requires the parameter p to be specified, and two synthesized signals with different p (respectively 0.1 and 0.6) are shown in Figure 5, hereinafter referred to as VSD1 and VSD2, respectively. Their corresponding statistical parameters are listed in Table 2 and their PSD in Figure 6. The kurtosis of the SDOF systems responses is shown in Figure 7. It can be observed that the kurtosis changes over the frequency range and proves significantly higher than the input kurtosis, to which extent depending on parameter p: the lower p, the higher the output kurtosis.

In the perspective of using the proposed algorithms to generate shaker input signals for fatigue life tests, further processing of the synthesized signals is carried out. For a straightforward comparative illustration of results, the duration of the hypothetical durability tests is ideally assumed to be equal to the duration of the reference field data. In particular, the FDS curves over the frequency range 0÷4000Hz are calculated and plotted with a resolution of 1/12^th octave. For brevity, only the FDS of the reference, MLV1, and VSD1 signals are shown. The most influent parameters in the FDS calculation are the Wohler’s curve slope (i.e. the exponent b in the well-known Basquin’s law NS^b = const.) and the damping ratio of the SDOF systems, which were set to 7 and 2.5%, respectively. The plots of the reference and the synthesized signals FDSs are reported in Figure 8. It should be evident that the Fatigue Damage curves are not overlapped. In order to adjust the FDS of the signals to match the reference one, the filter described in [20] is used (cf. Appendix C for an original mathematical proof of the filter effectiveness). The steps of the procedure are:

• calculate the FDSs of the reference and synthesized signals D_r(f_n) and D_s(f_n), respectively;
• define the spectral function: ;
• calculate the Inverse Fourier Transform of W(f_n) to obtain the impulse response of the filter;
• convolve the so-obtained impulse response with the synthesized signal.

Now, the synthesized and the reference profiles have the same FDS. It must be highlighted that the newly obtained signals are different from those synthesized previously and their statistical parameters and PSD change to some extent (Figure 9 and Figure 10, Table 3). However, their nature remains Leptokurtic and similar to the signals MLV1 and VSD1. The plots of the signals FDSs obtained after convolution with the filter are reported in Figure 11 and show a precise matching. It is worth noting that the signals in Figure 9 have a sampling frequency of 10240Hz since an up-sampling was necessary to correctly compute the FDS, as highlighted in [22, 23]. Though, the PSDs in Figure 10 are plotted in the range 0÷4096Hz, consistently with Figure 1b, Figure 3, Figure 6.

4. Discussion

The results show that the MLV algorithm synthesizes shaker input signals that can transfer their kurtosis value to the output of any SDOF linear system, irrespective of the latter’s natural frequency. This aspect eliminates the main flaw of Phase Manipulation and Polynomial Transformation methods. In particular, the kurtosis of the output features small fluctuations around a value close to the excitation one, due to the following reason: all the blocks constituting the synthesized signal have a wideband PSD, thus resonance effects cannot occur. The kurtosis transfer is only due to the energy of the bursts appearing in the signal, and the high excursions in the output cannot physically contain more energy. However, the input parameters of the algorithm are relevant and may affect the results. In fact, if the bursts duration (represented by parameter T_b) is not long enough, the system may not have time to respond and high excursions may not appear in the response. In addition, the MLV algorithm only deterministically imposes the phases of one block, thus achieving a high randomness. Consequently, different runs of the algorithm – with unchanged setup parameters – to process the same reference signal provide the synthesis of different time-series (all complying with the target PSD and kurtosis). Therefore, random vibration tests that require a sizable duration can be carried out starting from short environmental measurements by concatenating many profiles synthesized by the MLV algorithm.

The VSD algorithm can generate narrowband blocks, which can exasperate resonance effects on the DUT: the lower the value of parameter p, the more evident the resonance occurrence. The advantages of this algorithm are two-fold: (i) the synthesized signals are capable of generating SDOF systems responses with high kurtosis, though the excitation features relatively small value of the crest factor, thus subjecting the shaker to non-critical loads; (ii) kurtosis of the response is not constant over the natural frequency and is generally broadly higher than the excitation one (this scenario being quite common in many practical applications, due to the variable spectrogram of real-life vibrations). It must be highlighted, however, that a significant limitation of the algorithm could be that the PSD variation cannot be controlled. On the other side, its randomness contributes to the synthesized signal stochasticity, which is further increased by the complete aleatory generation of the phases in the IFFT transform for each block.

Both the algorithms could lead to signals with a very different damage potential from the excitation one, since kurtosis control does not directly control the FDS. The filter described in [20] was used to correct the FDSs of the signals generated by both the MLV and the VSD algorithms. The new signals showed similar Leptokurtic distributions as their corresponding synthesized ones, with slightly different PSD and statistical parameters. The procedure was used to synthesize signals having the same duration as that of the reference input, for the sake of direct comparison of results, but is actually intended to be used for accelerated fatigue-life testing.

5. Conclusions

In this work, two novel kurtosis control algorithms were proposed. The first algorithm, named Multi-Level Variance (MLV), splits the signal to be synthesized into blocks of the same duration and featuring variable variance. Namely, the shape of the PSD of each single block remains the same but scaled by an appropriate factor. The overall PSD and kurtosis approach closely to the reference input ones within a certain tolerance. The Papoulis’ Rule proves not to apply if the parameter that controls the duration of the bursts in the synthesized signal is appropriately chosen. In particular, the kurtosis of the output features the tendency to be constant over the natural frequency range of the DUT.

The second algorithm, named Variable Spectral Density (VSD), splits the signal into blocks of the same duration, not only with variable variance but in general with variable PSD. The overall PSD and kurtosis approach closely to the reference input ones, within a certain tolerance. The Papoulis’ Rule proves to not apply for every value of the only algorithm parameter that depends on the user’s choice. The kurtosis of the output tends to vary conspicuously over the natural frequency range of the DUT, with broadly higher values than those determined by the reference field data.

The damage potential associated with the shaker input signals synthesized by both algorithms proves different from that of the reference field data. The difference in the FDSs functions can be corrected by applying a special filter found in the literature. This leads to the generation of new profiles that match the reference FDS still featuring PSD and kurtosis values very similar to those of the signals synthesized by the algorithms. Hence, the overall procedure can be effectively employed to generate shaker input signals for accelerated fatigue-life testing.

The proposed algorithms have been validated through the application to many case studies, not reported in this paper for the sake of brevity (but a few examples can be partially considered in references [24, 25], where the formulation of the algorithms was not still optimized), starting from different reference field data featuring different characteristics – e.g. kurtosis value, bandwidth, spectral content, duration – thus proving their robustness for different scenarios.