The paper analyzes in detail the influence of the Welch method on the estimation of fatigue life in the frequency domain. The aim is to examine how the setting of parameters, such as window size, weighting and overlap, affects the calculation of power spectral density (PSD) and thus lifetime. The paper uses knowledge from the theory of signal processing, statistics and fatigue life estimation in the frequency domain. In order to demonstrate the general validity of the conclusions, four different reference PSD spectra are generated and examined. The spectra are converted to time histories using an inverse fast Fourier transform and then used as input (measured) data for the Welch method. The PSD thus obtained are subsequently subjected to lifetime estimation using the Dirlik and Tovo–Benasciutti models. Based on the obtained results, it was found that the best configuration of the parameters of the Welch method for fatigue calculations is the selection of the smallest possible time window (maximum 60% of the observation length) and the largest possible overlap (minimum 70% of the window size).
Keywords: power spectral density; Welch method; spectral fatigue; probability density function; sensitivity analysis
Digital signal processing is one of the most powerful technologies that has shaped science and engineering for the last five decades. The interfacing of measurement instrumentation and computers has brought significant changes in procedures such as data acquisition, digitization, storage, processing and evaluation [[
In practice, many mechanical structures and components (such as vehicles, wind turbines, PCBs, etc.) are exposed to random vibration during operation. In these cases, the assessment of fatigue life using frequency domain methods is advantageous because it provides a direct connection between the power spectral density and the damage intensity and significantly speeds up calculations compared to time domain analysis [[
If the time-varying load or response is in the nature of a stationary random Gaussian process, the power spectral density can be used to described it in the frequency domain [[
In parametric methods, the signal is assumed to be an output of a linear system driven by white noise, i.e., parametric methods estimate the PSD by initially estimating the parameters of the system assumed to generate the signal. Autoregressive (AR), moving-average (MA), or autoregressive moving average (ARMA) time-series models can be used to describe the system [[
The non-parametric methods refer to the methods of estimating the spectral density of a random signal without pre-parameter modelling and are therefore robust and much less sensitive to noise. They are based on the relationship between the PSD and the autocorrelation function, which are a Fourier pair [[
A third type of method used to estimate PSD are subspace methods. They are based on eigen analysis or singular value decomposition of the autocorrelation matrix. The subspace methods are most suitable for line spectra and are effectively used to detect sinusoids or to reduce noise in a signal. However, they produce only the so-called pseudospectrum. These methods include, among others, the multiple signal classification (MUSIC) method and Pisarenko harmonic decomposition method.
One way to estimate the life under variable amplitude or random loading in the time domain is by the rainflow counting algorithm, which allows to determine load cycles (amplitudes and mean values) from a random signal, first introduced by Matsuishi and Endo [[
Although the rainflow method is considered one of the most accurate for estimating fatigue life, several methods have been derived for estimating high-cycle fatigue in the frequency domain [[
The Dirlik method [[
The paper focuses on the correct estimation of the power spectral density in the sense of fatigue calculations. In the literature, it can be seen that the Welch method is the most common and even the only method used to determine the PSD for this application case. For this reason, only this method was chosen for deeper analysis. The Welch method requires several parameters that may affect the fatigue calculation process. One of the key parameters is the window type, which refers to the shape of the window function that is applied to each segment of the signal. Common window types include the rectangular window, the Hanning window, the Hamming window, and the Blackman window. Each window type has its own characteristics and may be more or less suitable for a given application, depending on the desired trade-off between spectral resolution and spectral leakage reduction. The window length is another important parameter, as it determines the length of each segment and therefore the frequency resolution of the resulting PSD estimate. A shorter window length will provide higher frequency resolution but may also result in a noisier estimate. A longer window length will provide a smoother estimate, but at the cost of lower frequency resolution. The overlap between successive segments is also a key parameter in the Welch method. A larger overlap will result in a smoother PSD estimate but will also require more computation time. A smaller overlap will provide a less smooth estimate but will be faster to compute. The optimal overlap will depend on the specific characteristics of the signal and the desired trade-off between smoothness and computational efficiency.
The main result of the article is to show which parameter of the Welch method affects the calculated fatigue life, also providing guidance on how to set these parameters to obtain reliable fatigue life estimation.
The PSD expresses the power contained in the narrow band
The basis of analytical calculation of PSD can be an autocorrelation function using Wiener–Chinchin relations. The power spectral density
(
where
The Welch method uses a modified version of the Bartlett method in which the segments of the series contributing to each periodogram may overlap each other. This means that the discrete time signal
(
where
In addition, the Welch method allows the data segment to be weighted by a time window function before calculating the periodogram. Therefore, the Welch method is also called the weighted overlapped segment averaging method or the periodogram averaging method. The result is a modified periodogram given by:
(
where
(
The Welch estimate of the power spectral density is an average of the modified periodograms of all segments:
(
where
Because the segments usually overlap, data values at the beginning and end of the segment tapered by the window in one segment, occur away from the ends of adjacent segments. This guards against the loss of information caused by windowing.
A weighting window using a trigonometric cosine function, which has a value close to zero at the edges, is named after mathematician Julius von Hanns as the Hanning window. It is one of the most-used time windows. This window is defined as follows:
(
The next well-known window designed by Richard W. Hamming is a modification of the Hanning window. It has an even greater effect on the dampening of the side lobes. It is defined by two coefficients
(
It most effectively suppresses the first sub-lobe when
The power spectral densities obtained in this way are prepared for the calculation of parameters, such as spectral moments
(
where
Another frequently used probabilistic model is the Tovo–Bennasciutti model, whose probability density function is [[
(
where variables
It can be noted that the Tovo–Benasciutti model uses a linear combination of the narrow-band
(
where
(
The fatigue-lifetime estimate in the frequency domain is obtained from the damage intensity
(
The aim of the experiment was to investigate how window size, overlap, and weighting affect the estimation of fatigue life using the Welch method to calculate PSD. The range of the parameters and types of the weighting function used in the calculations are listed in Figure 1.
In addition, to check the influence of the wave spectrum width on the power spectrum density estimates, 6 different PSD functions were used as waveform sources. The shapes of these PSD functions were selected on the basis of literature and standard recommendations, see Figure 2.
PSD was selected as the source of the time signals to ensure the same PSD shape for each time waveform during the simulation. Time courses were obtained by inverse Fourier transform with random phase shift. All calculations were performed in MATLAB.
For analysis purposes, six different PSDs were created and used in fatigue life calculations. The first three of them were (user-defined) synthesized functions with constant, linear and quadratic character, respectively. They are suitable to examine the sensitivity of PSD to changes in the parameters of the Welch method. As fatigue analyses are relatively widespread in engineering practice, the authors considered it necessary to give a practical example of PSD. To generalize the findings, calculations were performed for six different types of PSD functions. Therefore, three other functions were taken from the technical standards ISO 16750-3 and EN 61373. ISO 16750-3 deals with mechanical loads of road vehicles with a focus on environmental conditions and electrical testing for electrical and electronic equipment. EN 61373 deals with railway applications, rolling stock equipment, shock and vibration tests. All of the above PSDs are shown in Figure 2. The fatigue life results calculated based on them are considered the exact reference solution. The last two PSDs, shown in Figure 2e,f, are defined in logarithmic scale, so they had to be converted to a linear scale to be usable for further calculations.
Since the Welch method requires time signals, it was necessary to convert the reference PSDs to the time domain. In engineering practice, signals that contain a large number of samples are usually used to describe the time history, so the histories are smooth and well defined. To achieve the same effect, the length of the PSD corresponding to the number of frequency lines (bins) was set equal to the next higher power of 2, and then upsampled several times. This extended the acquisition interval of the time signal, which was subsequently obtained by inverse fast Fourier transform (IFFT), see Figure 3.
The sampling frequency was set so that the resolution of the frequency spectrum was 1 Hz in all cases. Since the PSDs were defined as one-sided, several adjustments had to be made before using the IFFT algorithm. First, a two-sided spectrum was calculated. A random number generator implemented in MATLAB was used to define the phase values. The generator uses the Gaussian distribution, so the results of the time history took the form of a Gaussian load. To avoid introducing inaccuracies into the calculations, the same set of random numbers was used for each signal. The conjugated symmetric vector thus obtained ensures that the IFFT provides a real value output. This approach also speeds up the discrete IFFT algorithm:
(
where
The kurtosis function was used to ensure the highest accuracy of normal distribution. The time histories obtained by IFFT are considered to be the input (measured) time signals. An example of the time history corresponding to PSD A is shown in Figure 3.
The MATLAB function "PWELCH" was used to calculate the PSDs. The DOE algorithm was based on loops, which ensured the automatic change of selected parameters of Welch method in the range given in Figure 1. This allowed a quick comparison of the results obtained. Within the DOE analysis, the "Full factorial" technique was used, the design matrix of which is a combination of all inputs in all specified values. All PSDs were further processed to assess the impact of the individual parameters. First, the effect on spectral moments was analyzed. The spectral moments can be calculated as [[
(
where
For
As can be seen, in general, the Welch method provides a sufficiently accurate estimate of the variance for a window size of up to 60% of all samples regardless of the character of the PSD. When larger segments are used, accuracy decreases exponentially. The plots in Figure 4 also show that increasing the overlap improves the estimate, with the highest accuracy being obtained with an overlap above 50% of the window size. Optimal adjustment of the window size and overlap reduces the variance of the periodograms. This effect can be also seen in the PSD plots shown in Figure 5, where the blue line represents the reference PSD.
The PSD of window size 20%, overlap 90% and Hamming window is shown by a red line. It is clear from the figure that the PSD calculated for these Welch method parameters shows a very good match with the reference PSD, in contrast to the PSD whose window size was 85% and the overlap was 25%. This PSD function (shown in yellow) shows significant local variations in value compared to the periodogram method, i.e. line smoothness is lost, which is undesirable when calculating the lifetime by spectral methods.
To confirm the general validity of the above findings, the same analysis was performed for the PSDs obtained using the Hanning window. The procedure was identical to the previous case with Hamming. The dependence of the zero spectral moments on the size of the window and the overlap is shown in Figure 6. A comparison of the PSDs obtained with different settings of the Welch method parameters and the reference PSD functions is shown in Figure 7. Based on the results, the same conclusions can be drawn.
The fatigue life calculation was performed for two probabilistic models—Dirlik and Tovo–Benasciutti. In both models, the S-N curve of a carbon steel given in [[
(
where
Estimates of total fatigue life determined by the Dirlik probability model for each type of PSDs are shown in Figure 9 and Figure 10. The plots show that if the window size is larger than 60%, the fatigue life increases exponentially for all cases. The accuracy of the estimate also decreases in areas where the overlap is less than 50% of the window size. It should be noted that the lowest value of the window size was 5% of all samples. However, testing has shown that it is not appropriate to reduce the window size further due to the loss of calculation accuracy. When taking a very small number of samples that enter the FFT calculation, the PSD values will be distorted and will not correspond to reality [[
The presented results of the life calculations were performed for the Dirlik probability model. Similar behavior was found when using the Tovo–Benasciutti model. The difference in the extent of the damage using the Tovo–Benasciutti model and the Dirlik model are on average 0.35% to 8.43%, which can be considered acceptable. The differences between the models depending on the changing parameters of the Welch method are graphically shown in Figure 11. Note that when calculating the fatigue life, it is necessary to emphasize the correct choice of the probabilistic model. If the calculation used, for example, Rayleigh's model, the results would be inaccurate, as this model is more suitable for narrowband signals. Therefore, the two most popular broadband load models were selected for fatigue analysis.
In order to verify the findings resulting from the sensitivity analysis described in the previous section, a FEM simulation of a simple component was performed. The analyzed object was a clamped beam of rectangular cross-section with two notches in the middle. A point mass of 50 g was applied at the free end of the beam, as shown in Figure 12. The beam was excited at the bond site. PSD type D from the previous analysis was used for the excitation. Since the lifetime estimation was performed in the frequency domain, i.e., using the spectral method, it is assumed that the excitation signal is stationary and ergodic, Gaussian-distributed and has zero mean value.
The lifetime calculation was performed in the FE-SAFE program. The input was data obtained from a preliminary simulation performed in Abaqus/CAE. The pre-simulation consisted of modal analysis and linear dynamic analysis. The output of the modal analysis was mode shapes, eigenfrequencies and modal stresses (Table 2). Three modes were identified in the range up to 500 Hz. The aim of the linear dynamic analysis was to obtain the generalized displacements and the phase angle of the generalized displacements for all modes (Figure 13). Note that unit excitation was considered in the calculation. From the response characteristics and modal stresses, real stress fields corresponding to unit excitation were obtained by modal superposition. If the actual excitation in the form of a PSD is subsequently included, it is possible to obtain stress PSDs, which enter into the fatigue life calculation.
A total of 380 lifetime calculations were performed, in which the parameters of the Welch method were changed in accordance with the previous analysis. For this purpose, a computational loop was created in the Isight program (Figure 14). The power spectral densities entering the lifetime calculation were taken from the previous analysis in MATLAB. The outputs from the FEM pre-simulation were the same in each loop cycle. For this reason, the location of the damage did not change. However, the degree of damage varied depending on the setting of the Welch's parameters. The location of the damage can be seen in Figure 15, which shows the result of the life calculation for the reference PSD. The reference value of the lifetime of the analyzed sample is approximately 43 s (10
The Welch method plays an important role in digital signal processing. It is commonly used to estimate the power spectral density. The Welch method reduces the variance of the periodogram method by averaging. From a practical point of view, the PSD is obtained by dividing the time series into overlapping subsegments, computing a modified periodogram for each subsegment and averaging the periodograms. Power spectral density is also an important characteristic in the analysis of fatigue life of components in the frequency domain. The aim of the paper was therefore to examine how the setting of Welch's parameters affects the results of such an analysis. Three parameters were considered: type of weighting window, size of subsegment, and size of overlap of subsegments. To generalize the findings, calculations were performed for six different types of PSD functions. The results were verified on the example of FEM analysis of the service life of a notched beam.
Based on the obtained results, it can be stated that the best configuration of Welch method parameters for PSD calculation and the subsequent determination of lifetime by the spectral method is the smallest possible time window (maximum 60% of all samples) and the largest window overlap (minimum 70% of window size). The results show that the lifetime values in this area are smoothed and close to the reference value. At this point, it should be noted that the above findings are valid for the probability models used, i.e., Dirlik and Tovo–Bennasciutti.
These findings are particularly important in terms of the research and development of mechanical systems and components, whose fatigue life is crucial to ensure their long-term and safe operation. It has been shown that incorrect parameter settings can lead to significant errors in lifetime estimation in the frequency domain. Spectral methods are often used in engineering practice, as they provide a relatively good estimate and are not time consuming. For this reason, it is important to point out the need for the correct transformation of time data to spectral ones.
Graph: Figure 1 Design of experiment flow chart.
Graph: Figure 2 Synthesized PSD using constants (a), constants with linear connection (b), and constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 3 Time history obtained from the reference PSD A (a) and from the upsampled PSD A (b).
Graph: Figure 4 Zero spectral moments (variance) depending on window size and overlap for PSDs calculated with using Hamming window according to synthesized PSD using constants (a), constants with linear connection (b), and constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 5 Comparison of PSDs calculated with different window size and overlay using Hamming window, blue line represents the reference PSD, red is for window size 20% and overlap 90%, yellow is for window size 85% and the overlap 25%. Comparisons prepared according to PSD using constants (a), constants with linear connection (b), constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 6 Zero spectral moments (variance) depending on window size and overlap for all PSDs computed using Hanning window according to synthesized PSD using constants (a), constants with linear connection (b), and constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 7 Comparison of PSDs calculated with different window sizes and overlay using Hanning window, blue line represents the reference PSD, red is for window size 20% and overlap 90%, yellow is for window size 85% and the overlap 25%. Comparisons prepared according to PSD using constants (a), constants with linear connection (b), constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 8 Idealized S-N curve defined using Basquin's equation.
Graph: Figure 9 Fatigue life estimates determined by Dirlik probability model (Hamming window). Comparisons prepared according to PSD using constants (a), constants with linear connection (b), constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 10 Fatigue life estimates determined by Dirlik probability model (Hanning window). Comparisons prepared according to PSD using constants (a), constants with linear connection (b), constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 11 Percentage difference between Dirlik and Tovo–Bennasciutti models (Hanning window). Comparisons prepared according to PSD using constants (a), constants with linear connection (b), constants with quadratic connection (c), standardized PSD functions according to ISO 16750-3 (d) and EN 61373 (e,f).
Graph: Figure 12 Specimen used in numerical simulation.
Graph: Figure 13 Linear dynamic responses.
Graph: Figure 14 Isight computational loop.
Graph: Figure 15 The results of fatigue life analysis for the reference PSD D.
Graph: Figure 16 Matrix of lifetimes for different Welch's parameters compared with the reference lifetime value.
Table 1 Mechanical properties and material constants of the material.
Tensile Strength Yield Strength Elongation E ν (MPa) (MPa) (%) (MPa) (-) (MPa) (-) 483 375 15 210,000 0.3 118 3.4496
Table 2 Modal analysis results.
Mode Mode Shape Eigenfrequency [Hz] 1. 39.01 2. 120.32 3. 466.72
Conceptualization, A.K. and A.N.; methodology, A.K. and A.N.; software, A.K.; validation, A.K., R.H. and A.N.; formal analysis, R.H.; investigation, A.K.; writing—original draft preparation, A.K.; writing—review and editing, A.N. and R.H.; visualization, A.K.; supervision, A.N. and R.H.; funding acquisition, A.K. All authors have read and agreed to the published version of the manuscript.
Not applicable.
Not applicable.
Not applicable.
The authors declare no conflict of interest.
By Adam Kaľavský; Adam Niesłony and Róbert Huňady
Reported by Author; Author; Author