Wavelet analysis to decompose a vibration simulation signal to improve pre-distribution testing of packaging.

In general, vehicle vibration is non-stationary and has a non-Gaussian probability dis- tribution; yet existing testing methods for packaging design employ Gaussian distributions to represent vibration induced by road pro ﬁ les. This frequently results in over- testing and/or over-design of the packaging to meet a speci ﬁ cation and correspondingly leads to wasteful packaging and product waste, which represent $15bn per year in the USA and € 3bn per year in the EU. The purpose of the paper is to enable a measured non-stationary acceleration signal to be replaced by a constructed signal that includes as far as possible any non-stationary characteristics from the original signal. The constructed signal consists of a concatenation of decomposed shorter duration signals, each having its own kurtosis level. Wavelet analysis is used for the decomposition process into inner and outlier signal components. The constructed signal has a similar PSD to the original signal, without incurring excessive acceleration levels. This allows an improved and more representative simulated input signal to be generated that can be used on the current generation of shaker tables. The wavelet decomposition method is also demonstrated experimentally through two correlation studies. It is shown that signi ﬁ cant improvements over current international standards for packaging testing are achievable; hence the potential for more ef ﬁ cient packaging system design is possible. 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC


Introduction
Mechanical damage in packaged products is commonly attributed to vehicle vibration and shock inputs. This damage arises as a consequence of fatigue due to packaged-product resonance or as a result of high level shock events [1]. Both forms of damage are possible due to the random, non-stationary nature and non-Gaussian distribution of vehicle vibration [2]. The pre-distribution testing of packaging plays a crucial role in mitigating the risk of damage by ensuring an appropriate level of protection is provided. In recent years, the method used to simulate vibration during testing has been subject to much debate over its inability to realistically simulate vehicle vibration [3,4]. Such deficiency and corresponding uncertainty may result in: (a) The packaging system underperforming if qualified against an unrepresentative test input; or (b) Overdesigned packaging if the test inputs are above thresholds that are considered appropriate for real road profiles or a particular distribution network. Both of these deficiencies lead to waste in terms of packaging material and product. In the U.S. alone, unsaleables cost retailers and brand owners $15 billion/year [5]. In Europe, a 2004 study by Smithers Pira estimated the cost of damage in Europe to be between €2bn and €3bn per annum [6].
A time signal, xðtÞ, may be decomposed into Fourier harmonics according to which has a periodicity of T ¼2π/ω 0 and harmonics defined at the discrete frequencies, f n ¼ nω 0 /2π. The discrete harmonics, X f n À Á , become continuous Fourier transform coefficients in the limiting case of T-1. In terms of the discrete frequencies, the power spectral density (PSD) of the signal is defined by where Δf ¼ 1=T. If the time signal is sampled at t r ¼ rΔt, where Δt¼T/N, the Discrete Fourier Transform (DFT) may be used to evaluate the harmonics using If an ensemble of Q signals is considered, an average PSD may be evaluated [1]: An essential feature of current testing regimes is the ability to perform accelerated testing. This is important for testing houses who need to maximise utilisation of facilities and offer competitively priced services, and for manufacturers who demand shortened test and development cycles. To reduce long test durations the average PSD may be time compressed (accelerated). This is based on the assumption that, in the elastic region, fatigue life of a material is linearly accumulative and that the relationship between stress cycles and magnitude of stress follows the relation given by [7] m where N f is the fatigue life at stress level S l , γ is the slope factor of the curve for the given material, and m is the material constant.
As previously stated, if it is assumed that fatigue is linearly accumulative, an approximation for the cumulative effect of stress cycles at different stress levels can be given by [8] where n i is the number of stress cycles at stress level i, N f i is the fatigue level (cycles to fatigue) at stress level i, and F is the number of stress levels in a stress cycle history leading to failure. As described in [9], expanding Eq. (6) with substitution of Eq. (5) and the time for a specified number of cycles, a simplified relation between test duration and intensity, and actual journey duration and intensity is possible. This can then be used to compress a vibration test [10]. The increase in test intensity is inversely proportional to the reduction in test duration following the power law, where T a is actual journey duration, T c is compressed test duration, a a is actual journey intensity and a c is compressed test intensity (where intensity is represented as the root mean square with respect to g). For the purpose of testing packaging, γ is generally set to a value between 2 and 5 with 2 being the worst case, and the time compression factor to 5 [11]. The PSD for the compressed test is generally applied to a test specimen using a shaker table with a standard random vibration control strategy. An input vibration signal may be constructed from Fourier coefficients associated with randomly distributed phases. In contrast to actual vehicle vibration, the signal produced by the averaging process is stationary with a probability distribution having a particular level of kurtosis. For the purpose of simulation on shaker tables the level of kurtosis is an important parameter to consider. This is because high levels of kurtosis typically signify the presence of large variations in acceleration levels and correspondingly, given the capabilities of the shaker table, the high level accelerations may not be accurately reproduced. Consequentially, the varying intensity and discrete high amplitude shock events, which are prevalent in vehicle vibration, are typically lost. This is further compounded by the fact that ISTA [12]/ASTM [13] random vibration test spectra are all based on a 3-sigma-limited Gaussian distribution with a kurtosis, β 2 ¼ μ 4 =μ 2 2 ¼ 3, where μ i is the i-th central moment of the distribution. For reasons previously stated, current simulation methods necessitate a need for conservatism in packaging design. While it should be possible to simulate a time replication signal that includes all shock events and varying intensity of vibration, hardware (acceleration) limitations of current shaker tables means that this is not generally possible. An alternative approach is to seek an improved compromise between the relative accuracy of the simulated test with respect to the original (actual) journey, and the ability to simulate the signal through existing vibration table systems. One recently proposed means to overcome this problem is to modify the kurtosis directly. Such an approach is justified because kurtosis has been shown to modify the crest factor and to do so in the same sense i.e. positive change in kurtosis results in a positive change in crest factor and vice versa [14]. In addition to a theoretical justification for modifying kurtosis a number of shaker table controller manufacturers now offer kurtosis control [15], which allows a user to select a target value of kurtosis, typically in the range from 3 to 9. While the potential of kurtosis control is apparent the aforementioned challenge of achieving the best compromise between the accuracy of the simulated test with respect to the original, and the ability to simulate the signal through existing vibration table systems, where it is desirable to reduce kurtosis levels to lie within certain limits or below a maximum threshold, remains unresolved. In general, the use of a test signal having a single kurtosis value will incur higher levels of accumulated input from the effects of multiple transient events. Accordingly, this paper describes the use of wavelet analysis to evaluate time variations of frequency spectra in order to decompose a signal into components with separate kurtosis levels. The highest kurtosis signals are contained within outlier components, which are applied over a shorter sub-period of the full test signal period. The theory of wavelets is first discussed and the selection of a suitable mother wavelet is summarised. The proposed wavelet decomposition method is described in detail with an illustrative example. The proposed method is then verified and validated for two journeys through consideration of damage correlation with time replication.

Wavelet analysis
Although there exist alternate approaches for decomposing time-series signals we have selected Wavelets because of their ability to localise in time and frequency. Importantly, this enables the decomposition, post-processing and then reconstruction of a signal, via segments of the signal that possess particular characteristics (features), In contrast to applying a uniform operation to a signal, such as a transform or threshold, we are able to extract those segments that can be simulated with existing hardware and then post-process the remainder. This approach thus offers the opportunity to achieve the aforementioned compromise of an improved and more representative simulated input signal that can be used on the current generation of shaker tables. 'The wavelet transform is a tool that cuts up data, functions or operators into different frequency components, and then studies each component with a resolution matched to its scale' [16]. Wavelet analysis is used in many applications because of its effectiveness in detecting facets in signals. It is used in the decomposition of both electrocephalagraph (EEG) and electrocardiogram (ECG) signals to identify transients with respect to both the time and frequency domain [17,18]. Further applications include fault diagnostics in mechanical machinery and crack identification in structures [19,20]. The dilation of a mother wavelet, ψðtÞ, is determined by a scale, a, which is inversely proportional to frequency. A translation parameter, b, is then associated with translation (time). In general, the scaled and translated wavelet is normalised: The dilation of the wavelet allows the resolution of the analysis to vary, which is a benefit over use of the Short Term Fourier Transform (STFT). It enables a lower frequency resolution, but a higher temporal resolution when analysing high frequencies (low scales). Conversely, it affords a higher frequency resolution, but lower temporal resolution at low frequencies (high scales). This property of wavelets is useful when analysing signals such as vehicle vibration, where low frequency vibration tends to be predominant and continual throughout the signal, while high frequency events tend to occur over short durations.
The continuous wavelet transform (CWT) is found by calculating the convolution of the wavelet function, ψ a;b ðtÞ, and the time domain signal, xðtÞ, using where Cða; bÞ is the wavelet coefficient at scale a and translation b. The discrete wavelet transform (DWT) follows by setting C jk ¼ Cð2 j ; k2 j Þ, where j ¼ 1; 2; …; J and k ¼ 1; 2; …; K.
The most appropriate mother wavelet for use with vehicle vibration analysis needs to be selected to enable efficient identification of key features in a signal. One candidate is the Morlet wavelet, which has relatively good time and frequency resolution capability [21]. Furthermore, the Morlet wavelet has been used successfully to identify impulsive features within a signal, which is relevant for the capture of shock response events in vehicle vibration [22]. It is commented that the requirement is to characterise the on-board vehicle response, rather than a direct shock input, which would be filtered by the suspension/vehicle dynamics. Alternative forms of mother wavelet such as the Daubechies and Meyer wavelets were also considered but due to the requirement for a continuous wavelet, and given the expected characteristics of the response signal, the complex form of the Morlet wavelet was adopted. The Morlet wavelet is given by Step 8 Calculate the kurtosis, β2, of all signals ( m= 1,2,…, M and )

No
Step 6 Extract internal and outlier components

Vibration signal
Step 2 Construct stationary signal ( xm) Step 3 Calculate the wavelet coefficients of ym and xm ( and ) Step 4 Construct a decomposition envelope defined by Step 5 Identify wavelet coefficient outliers for ym Step  where ω 0 is the central frequency of the mother wavelet. For real vehicle vibration signals only the real part of the Morlet wavelet is required, which may be adapted to whenever ω 0 Z 5. In this case the mean error from zero is negligible.

Wavelet decomposition for pre-distribution testing
Used independently, wavelet analysis provides an appropriate method to identify transient events in a signal. However, it does not directly offer a suitable method of producing a simulation test for application with shaker tables and one that can be compressed (accelerated). In this paper, this limitation is overcome using a combination of time-frequency wavelet analysis and frequency analysis tools. More specifically, the Morlet wavelet is used to identify changes in vehicle vibration, allowing the signal to be decomposed into a number of stationary components. Each component can then be analysed and represented using its PSD and duration. By applying the PSD of each component in terms of a vibration signal having appropriate spectral content with randomised phasing, an overall simulation signal may be constructed that is more equivalent to the non-stationary and non-Gaussian probability distribution of the vehicle vibration.
Practically, the method involves nine steps and two stages of decomposition. The first stage is used to decompose the original signal into a number of components (Mþ1), where components m ¼1 to M have lower kurtosis values than the original signal. The second stage is then used to evaluate and further decompose each component, to further reduce the kurtosis of each component. During each decomposition wavelet analysis is used to separate a vibration signal into components: (1) One with a lower level of kurtosis; and (2) One with a higher level of kurtosis.
In general, both will have non-Gaussian probability distributions (kurtosis a3), while Eq. (2) may be considered as a residual or outlier component. At the end of the iterative process, a simulation signal is constructed from PSDs, to generate randomly phased component signals.
The purpose of the decomposition is to enable a measured non-stationary vibration signal to be replaced by a constructed signal that includes as far as possible any non-stationary features in the original signal. The constructed signal consists of a concatenation of decomposed shorter duration signals, each having its own kurtosis level. A threshold level of kurtosis, β 2T , may be used to terminate the iterations required for the decomposition process. Fig. 1 illustrates the nine step process to achieve decomposition of the original signal and reconstruction for simulation.

Step 1: Time-averaged signal PSD
In principle, the DFT of Eq. (3) may be applied to each signal, y m ðtÞ, to yield the harmonics Y m ðf n Þ. However, the sample length may be excessive. Alternatively, averaging over Q windowed segments, each having sample length N, may be applied in the form The time-averaged signal PSD is then calculated from the average spectrum [1]: where the one-sided Fourier spectrum is used.

Step 2: Construction of a stationary signal
A signal, assumed to be acceleration for vehicle vibration analysis, may be constructed from the PSD of Eq. (13) having the same overall RMS acceleration and duration as the parent signal, y m ðtÞ. As the constructed signal would be derived from averaged data, the RMS acceleration will be constant. Therefore it will not recreate any non-stationary features of the original vehicle vibration, nor will it have a non-Gaussian probability distribution. However, it will form the starting point of an iterative process to decompose it into other components that may be usefully implemented for simulated test signal generation (Steps 3-9).
A stationary signal may be created by inverting Eq. (13). Firstly, the PSD is used to derive the Fourier harmonics where the phases φ n n ¼ 1; 2; …; N ð Þ form a randomly chosen set. A time domain signal segment of period NΔt can then be constructed: since Xðf À n Þ ¼ Xðf n Þ Ã . In order to then obtain a time domain signal with the same duration, L, as the original vibration signal, the process in Step 2 needs to be repeated L=NΔt times. The constructed signal parts are concatenated to form a time domain signal, x m ðtÞ. An example of the stationary signal, x m ðtÞ, derived from a recorded acceleration signal, y m ðtÞ, is illustrated in Fig. 2.

Step 3: Calculation of the wavelet coefficients
For the purpose of decomposition, wavelet analysis is performed on both the original vibration signal, y m ðtÞ, and the constructed stationary signal, x m ðtÞ. The Morlet wavelet coefficients are evaluated according to Eqs. (9) and (11) with discrete scales, a ¼ 2 j , up to level j ¼ 6. These are denoted by C y;m ða; bÞ for y m ðtÞ, and C x;m ða; bÞ for x m ðtÞ. These wavelet coefficients will form the basis of decomposing the original signal into a component that has a lower kurtosis than the original signal, together with a remainder component that possesses characteristics that increase the spread of the original signal's probability distribution, and correspondingly its kurtosis.
By considering all time translations, b, this envelope provides a mechanism to partition the original signal, y m ðtÞ, through assessment of where its wavelet coefficients lie (Step 5). An example of a decomposition envelope, derived from the signal shown in Fig. 2, is given in Fig. 3. The scales used are in the range 0 r a r64.

Step 5: Identification of the outlying coefficients
The decomposition envelope is compared with the wavelet coefficients of the original vibration signal, y m ðtÞ, at all translations, b. Any coefficients that are outside the envelope's maximum or minimum values are identified according to: These coefficients are termed outliers. The corresponding wavelet coefficients that lie inside the envelope are denoted by C I y;m ða; bÞ. Regions without outliers indicate the periods in the original signal that are within the bounds of the constructed stationary signal.  where y I m ðtÞ and y O m ðtÞ denote the components derived from C I y;m ða; bÞ and C O y;m ða; bÞ, respectively. These components are obtained by inverting the wavelet transforms of the partitioned wavelet coefficients from Step 5.
The data windows surrounding locations that exceed the envelope are extracted from the vibration signal. The size of the window is determined by the frequency resolution required. As vehicle vibration is typically analysed within 1-200 Hz, a minimum data window size of 1 s is deemed appropriate.
The extracted portions of the signal, and similarly the remaining portions of the signal, are concatenated to form two components, the outlier component and remaining inner component i.e. falling within (inside) the envelope. An example of this, for the signal in Fig. 2, is given in Fig. 4, where (a) is the remaining component, and (b) is the outlier component.
and Steps 1-7 repeated. This will yield the expression of the original signal in the form Note that the decomposition envelopes will generally differ for each component. Theoretically, following the iterative process, a vibration signal may be decomposed into an infinite number of stationary components. However, as previously mentioned, practicalities and the limited window size constrain this process. The iterative process continues until either:

Step 8: Further decomposition of components
If the kurtosis of particular components remains above a threshold value, then a second level of iteration may be applied to the components. Steps 1-7 may be applied to form For the worked data in Table 1, a second level of decomposition was carried out using a limiting number of S ¼ 1 iterations per component (M ¼ 3). Table 2 shows the pertinent results for the 11 components. It is evident from Table 1  , which indicates that these components again contain the shock-like acceleration inputs.

Step 9: Construction of a simulation signal
In the construction of the simulation signal (termed constructed signal) the PSD for each of the decomposed components is first generated. This results in a stationary signal for each component over its evaluated duration (c.f. Table 2). These components may then be concatenated to form a simulation signal suitable for vibration testing. For the previously considered data relating to Table 2, an example of the constructed simulation signal is shown in Fig. 6. The original signal is given in Fig. 5 for comparison. The concatenation process requires a method of ordering the components in Table 2 to be selected. For the purpose of this study, the components have been ordered according to the sequence of decomposition. However, in practice and depending on the intended controller and shaker table setup, it may be beneficial to order them so as to minimise demands on the system. For example, components may be concatenated in ascending or descending order of RMS. Further, where there are large steps between concatenated components it may be desirable to blend the signals to avoid any transients.
The wavelet decomposition method enables a non-Gaussian and non-stationary vibration signal, to be decomposed in to a series of components with comparatively lower kurtosis. By then constructing a signal on a shaker table, from each component's PSD, a simulation signal can be constructed that more accurately represents the original signal, than the current average PSD method defined earlier. To assess the appropriateness of this method, comparisons between the average PSD and RMS distribution of the original and constructed signals are given in Table 3 and Figs. 7 and 8, respectively.
From these results, it can be deduced that the PSD of the constructed signal shows a good overall correlation with the original signal, which is the aim of the proposed method. For completeness the percentage distributions of RMS levels for the original and constructed signals are given in Fig. 8. The distributions are significantly different, which is to be expected. That is, signals can possess a similar PSD but drastically different RMS distributions and there is no expectation that they should be similar. Of particular note is that the concentration around 0.27 g to 0.29 g is caused by the simulation of part 4.4, which contains all of the high level events from the signal.

A correlation study
Correlation studies are an important part of validating a simulation method as they allow the assessment of a simulation method's ability to recreate actual damage. There are many mechanisms of damage for packaged-products, e.g. scuffing, denting or distortion, each of which can be described qualitatively, but cannot easily be translated into a quantitative  Table 2 of the original signal in Fig. 5.   measure. It is consequentially difficult to quantitatively validate simulation methods. Further there is an added complexity of the inter-dependency of various damage mechanisms, meaning that the isolation and measurement of a single damage mechanism is non-trivial. In order to address this, [9] created a device for measuring the damage mechanism of scuffing. The rig was designed to: 1. Provide a measurable response to vibration input, within the range of vehicle vibration (1-200 Hz). 2. Represent a damage mechanism that is widely observed in the distribution of packaged goods. 3. Provide a repeatable measure with respect to variation in both input intensity and input duration.
This rig was further used in a correlation study by [23] to compare the appropriateness of various simulation methods. While 'scuffing' represents only a single damage mechanism the study demonstrated that, as a comparative test, scuffing was both repeatable and a good indicator of the closeness or representiveness of simulated signals to real-time vibration signals. Thus, in order to evaluate the proposed wavelet decomposition method of this paper, the damage mechanism of scuffing is considered.

Method of measurement of scuffing
A detailed description of the rig used to measure 'scuffing' and its validation are given by [9]. A summary is provided herein for completeness. The rig used is illustrated in Fig. 9.
A sheet of clay-coated board is printed with a flexographic ink, 'blue 10' from Sun Chemicals [9]. This is mounted on a 25 kg cubic mass. An unprinted sheet of Ivorex high white uncoated board is then affixed to the surface of a steel container ( Fig. 9 part a), inclined at 45°to the horizontal. The face of the cubic mass with the printed board on is then placed onto the inclined surface so that the printed and unprinted board are in contact. Foam blocks of polyethylene (PE) with a closed cell density of 75 kg/m 3 are placed below and above the mass, so as to respond as a packaged item with cushion shock protection (Fig. 9, parts a and b). The mass is then fully enclosed within the steel frame container designed such that resonant frequencies are well above 200 Hz (Fig. 9, parts c and d).
Following testing, the Ivorex board is removed and the residual ink transfer is determined by colour scanning the board and measuring the 'blueness' of the corresponding array of pixels (Fig. 10). This results in an overall percentage measure of 'blueness' and hence scuffing.

Journey 1 -Vibration signal from a motorway
In order to provide a means of comparison of the wavelet based decomposition method to other methods, the measured acceleration time history shown in Fig. 11(a) was adopted. This signal is taken from [23] where it was used to compare the established ISTA method [12] with the two-way split spectra method; the three-way split spectra method; and RMS modulation. With reference to Eq. (20), the signal of Fig. 11(a) was decomposed into 21 components with M¼7 and Sr2. The constructed simulation signal is shown in Fig. 11(b). The kurtosis, RMS values and duration of each of the PSDs used to construct the simulation signal are provided in Table 4. Two components (8 and 21) have high kurtosis (13.5 and 13.2). The high kurtosis is due to the relatively long duration of these constructed components, enabling high level events from the  original signal to be included with lower level vibration, resulting in a wide spread distribution, concentrated around zero. The kurtosis of these segments could likely be reduced through further decomposition. A comparison of the overall kurtosis and RMS acceleration of the original and constructed signal is given in Table 5. The difference between the RMS accelerations is 0.002 g, which is relatively small and would therefore suggest good correlation between the two signals. The constructed signal has a kurtosis of 5.6, which is a significant reduction from the original signal that has a kurtosis of 7.4. Given that the aim of the proposed method can be considered to involve reducing the kurtosis while retaining the PSD it is thus important to compare also the PSDs. The PSD of the original signal and the constructed signal is given in Fig. 12.
The overall shape of the constructed signal's PSD correlates strongly with the original vibration signal's PSD. There are small differences in the frequency ranges 7-8 Hz and 300-400 Hz. These differences arise due to the averaging process when calculating the PSDs for each segment. Testing was carried out using the same signal as [23] and also a sample of Ivorex coated board from the same batch used in [23]. Importantly, this enables the results from the current study to be compared directly with other methods tested in [23]. The percentages of scuff damage produced from a time replication simulation and wavelet decomposition method are shown in Table 6 and presented in Fig. 13 alongside the results from [23].
In [23], six test methods were considered including: Time replication [1]; High and low separate two-way split spectra [24]; Three-way split spectra [25]; Shock on random [26] and Root mean square (RMS) modulation [27]. The results in Fig. 13 show a very strong correlation between the two replication simulations with damage levels of 0.59% and 0.60%, respectively. This demonstrates the repeatability of the correlation study. The level of scuffing produced using the wavelet decomposition method correlates well with both time replication tests, with the difference in the level of scuff damage being 5% and 7% (relative) and 3% and 4% (absolute), respectively. When compared with the other simulation approaches, the wavelet decomposition method correlates directly with the modulated RMS approach, each achieving 0.56%.
Although small ( $ 5%) there is some difference between the measured damage levels for time replication and the wavelet decomposition method. This is to be expected given that the scuff rig had a confidence limit of 95% with an error margin of 5% [9]. Consequentially, a tolerance of 5% is applied in Table 6 and Fig. 13, and gives rise to a significant region of agreement between the time replication and wavelet decomposition method. In order to explore this further the wavelet decomposition method was used to simulate an additional, more extreme, journey.

Journey 2 -Vibration signal from minor roads
The time-history for the test vehicle vibration is given in Fig. 14 and was obtained on minor roads in the South West of the UK, in a Ford Luton Box van with the vehicle travelling between 40 and 60 mph.
The condition of the minor roads was severe, with significant road damage and potholes, due to a lack of maintenance. This, combined with the speed, resulted in a time-history vibration containing a large number of shock events (events over 1 g).
The percentage scuff damage produced using wavelet decomposition and time replication of the time-history, is shown in Table 7.
In this test the wavelet decomposition method correlated well with the time replication with only 8% difference in the level of scuff damage.

Discussion
A wavelet decomposition method has been used in the process of constructing a simulation signal from a measured vehicle vibration signal. It has been shown that by simulating a signal as a series of stationary components with varying kurtosis/RMS levels and spectral content, it is possible to reproduce a greater part of the original non-stationary characteristics. That is to say that the effect of averaging is significantly reduced enabling a greater proportion of the high intensity events to be included, and thus more representative simulation signals created.
The proposed wavelet decomposition method has also been evaluated experimentally for two journeys (Sections 4.2 and 4.3). In both cases the induced damage levels for the constructed signal correlates well with time replication simulation with a difference of between 5% and 8%. In the first journey the method is further compared to existing approaches. It is shown that the method performs similarly to Modulated RMS, which outperforms the other methods of high and low separate twoway split spectra, three-way split spectra, shock on random and RMS modulation; significantly with a simulated scuffing level of 0.56% compared to 0.59% for time replication.
One of the key tenets of this approach is that it can be used on the majority of existing shaker tables. While the approach produces a significantly more representative simulation signal, the iterative decomposition can result in at least one final component with a high kurtosis (49). These final components (components 8 and 21 in Table 4) have kurtosis values of 13.5 and 13.2, which for the purpose of simulation, increases the likelihood of the signal being moderated. That is, parts of the signal and in particular high-level shock events will be replicated with much lower amplitude vibration. Further iteration of these components may result in several subcomponents that have a distribution with a lower kurtosis, however, due to the nature of vehicle vibration, the final component(s) will likely always contain discrete high level shock events. Although not the primary motivation of the proposed approach, the residual components can be considered to represent the majority of the high-level shock events. Hence, an alternate application of the proposed approach would be to identify and extract the high-level shock events which could then be simulated separately e.g. through drop test.
Whilst the nature of the vibration limits the accuracy of the stationary components, the size of the data window (Section 3.6) used in the decomposition of the signal also has an effect. To enable appropriate accuracy for frequency analysis of each component, a minimum window size used for decomposition needs to be chosen. For the purpose of analysing vehicle vibration a data window size of 1 s has been used. This allows for frequency analysis from 1 Hz.

Conclusions
Vehicle vibration is inherently non-stationary and has a non-Gaussian probability distribution; yet existing testing methods for packaging design are based on a single Gaussian distribution. This can lead to either over or under conservatism in packaging design, which in turn may lead to excessive material or product damage, respectively.
To address this deficiency, a method has been presented that uses time-frequency (wavelet) analysis to decompose a non-Gaussian and non-stationary signal into a series of stationary components. The approach allows an improved compromise between the relative accuracy of the simulated signal with respect to the original. Further, and of particular relevance to industry, is the ability to execute the simulated signal on existing vibration shaker controllers, and thus allow for test compression (acceleration). The suitability of the Morlet wavelet for vehicle vibration decomposition is discussed. A nine-step decomposition strategy has been presented where a vehicle vibration signal is iteratively decomposed into components that are stationary. These are combined to construct a simulated signal for testing. The suitability of the method demonstrated by comparison of the PSD and overall RMS of the results of a decomposed simulation signal with that of an original journey signal which reveal a good correlation.
In order to validate the utility of the approach the decomposition method was applied experimentally in a correlation study of scuffing for two journeys: a motorway and minor roads. In both cases the measured damage levels were within 8% of those measured for the time replication studies. Additionally, for the first journey, the results were compared with a range of typically applied simulation approaches where it was shown to outperform all, and match the Modulated RMS approach.
Overall, the results show that the method is effective in creating an improved and, in particular, a more representative simulated input signal. Further the utility of the approach in constructing a simulation signal that can be used on existing multi-axis shaking tables has been demonstrated. Satisfying this latter constraint is essential for industrial/commercial organisations to apply regimes for accelerating testing and thus provide commercially viable testing services.