Evaluating internal condition of hardwood logs based on AR-minimum entropy deconvolution combined with wavelet based spectral kurtosis approach

: The aim of this research was to explore the potential of acoustic impact test to evaluate the condition of hardwood logs in regard to internal decay, void, crack and defect ratio using an acoustic signal separation and enhancement algorithm. Longitudinal acoustic signals were obtained from 15 logs of four hardwood species through acoustic impact testing. The defect components were separated from the acoustic response signals and enhanced based on the autoregressive minimum entropy deconvolution (AR-MED) method, and from which the kurtosis was derived and used as the global feature parameter for evaluating the internal condition of logs. Compared with the acoustic velocity obtained directly from the original signal, the kurtosis was deemed to be a more powerful predictor of log defect ratio with higher coefficient of determination ( R 2 = 0.89) and was not affected by log species. To identify the type of defects, a complex Morlet wavelet-based spectral kurtosis (SK) method was proposed. The research results indicated that the SK can not only determine the type and primary and secondary major defects, but also be able to identify those that were not detectable by global acoustic parameters.


Introduction
Hardwood logs are ideal raw materials for making highvalue wood products such as furniture, flooring, and veneer because of its rich texture and colors and high mechanical properties. As a biological material, hardwood logs typically contain some internal defects such as decay, holes, knots, wounds, and other growth defects. The severity and frequency of defects vary widely within species, and even within the same tree. Early detection of these internal defects in hardwood logs could provide significant benefits to wood industry in terms of saving transportation cost, reducing wood waste, and increasing production efficiency; and ultimately to achieve optimal wood utilization and maximize the profits (Edlund et al. 2006;Fischer et al. 2015;Wang et al. 2009;Xu et al. 2018).
Another potential market for hardwood materials is building industry. Currently, mass timber building is on the rise world-wide and creates increasing demands in wood supply. However, despite its superior mechanical performance, hardwood is still less used in timber structural applications (Schlotzhauer et al. 2019). For example, surveys (Aicher et al. 2018;Kobel et al. 2014) showed that most of the glulam products are made of softwood worldwide. One of the contributing factors is the large quality variations existed in hardwood logs and the significant cost associated with processing the logs. Studies have indicated that, if information of internal defects in hardwood logs (such as location, type, and size of defects) is known and utilized effectively during the sawing process, the lumber value can be increased by 10-21% (Steele et al. 1994;Thomas 2010).
The use of acoustic techniques for detecting internal defects in hardwood logs has been investigated and the preliminary results are promising. However, the accuracy of defect diagnosis is far from satisfactory (Koca et al. 2018;Wang et al. 2004Wang et al. , 2009. Xu et al. (2018) examined an acoustic impact testing combined with advanced waveform analysis to assess the internal condition of yellowpoplar logs. They reported that acoustic velocity did not have a clear correlation with the board grade yields; on the other hand, time centroid and damping ratio of the acoustic signals were found somewhat effective in predicting log quality. In a laboratory investigation, Xu et al. (2019) also found that the acoustic velocity was not sensitive to small defects, but it had a negative nonlinear correlation with the defect ratio, an overall quality parameter for hardwood logs. They concluded that time centroid and the second-order damping are better indicators of log quality than acoustic velocity.
Although several acoustic parameters have been proved to be valid quality indicators of hardwood logs, each has its strength and weakness because the physics and mechanism that relates the acoustic measures to various defects are different. It is desirable to extract new acoustic features as complementary evaluation means to further improve the accuracy of defect detection in hardwood logs.
The difficulty for internal defect detection or quality sorting of hardwood logs lies in the non-stationary and spectrum overlapping of defect components obtained from the acoustic response signals. The measured signal is usually processed directly in time-or frequency-domain to obtain acoustic parameters in traditional acoustic sorting. Thus the parameters reflect only the overall characteristics of a signal, not defect-induced instantaneous features that are usually masked by the dominant stationary components. As a result, the quality sorting of hardwood logs using the traditional method is not very accurate, especially for the logs with minor or early defects. It is assumed that if the defect components could be separated from the response acoustic signals, then a more accurate quality evaluation of hardwood logs could be achieved by extracting the features of the defect components.
Autoregressive (AR) model is one of the most basic and widely used time series analysis methods. Research indicate AR model is sensitive to system state variations (Al-Bugharbee and Trendafilova 2016; Cao and Wang 2016). Wang and Wong (2002) established an AR-based filter to successfully separate the crack-induced transient components from the gear meshing vibration signal. Tian et al. (2013) constructed an AR model to isolate mutation signal caused by tower crane damage under impact load. The results showed there was a positive correlation between the variance of residual signal and damage degree. Although removing well the stationary components from the measured signal, it is difficult for AR model to distinguish noise from singular signal as the minimum phase condition assumed in the AR model (Endo and Randall 2007).
Minimum entropy deconvolution (MED), a kind of blind deconvolution technique in the time domain without any priori assumptions, was originally presented by Wiggins (1978) and successfully applied in seismic signal processing. The method can improve the signal-to-noise ratio, enhance mutation components of the analyzed signal, and make the output signal approximate defect components as possible, so it is suitable for the early defect detection.
For a weak transient signal caused by defects, it can be further enhanced by high-order statistics through restraining the noise. Kurtosis is a numerical statistic reflecting the characteristics of signal distribution and sensitive to the variation of signal variance and amplitude, which has now gradually applied to feature extraction of dynamic signals (De la Rosa and Moreno Muñoz 2008; Hayashida et al. 2014). The amplitude of acoustic waves transmitted in sound wood is the approximately normal distribution, and its kurtosis value fluctuates less. But if there are some defects, such as void, rot or crack, etc., in the log, the amplitude distribution will deviate from the normal curve due to the superposition of defect components, as becoming sharp or flat, and then the kurtosis will change accordingly (De la Rosa and Moreno Muñoz 2008; Mucciardi et al. 2011). Therefore, the change of kurtosis could be used to determine whether random components are introduced into the signal.
So the accuracy of quality evaluation of hardwood logs could be improved through achieving firstly the separation and enhancement of defect signal based on the autoregressive minimum entropy deconvolution (AR-MED) method and then extracting kurtosis from the separated components as the evaluation parameter.
Global variables such as time centroid, damping ratio, etc. can determine whether defects exist in a log but not identify the defect types for the fact that no discrimination being made among the frequency bands (FBs) of the acoustic signal (Xu et al. 2018(Xu et al. , 2019. The application of spectral kurtosis (SK) makes possible for solving the problem (De la Rosa and Moreno Muñoz 2008; De la Rosa et al. 2015). In essence, SK is to calculate the kurtosis at each frequency line for revealing the hidden nonstationary components in signal and determine in which FBs they occur.
This study was based on the acoustic evaluation of 15 hardwood logs reported in a previous paper (Xu et al. 2019). The aim of the research was to further explore the potential of acoustic impact test to evaluate the condition of hardwood logs in regard to internal decay, void, crack and defect ratio using an acoustic signal separation and enhancement algorithm. This paper specifically focused on 1) extracting kurtosis from the acoustic response signals using the AR-MED method coupled with the wavelet-based spectra kurtosis approach; and 2) Examining the relationships between kurtosis and spectra kurtosis and the internal quality of hardwood logs in terms of defect ratio and defect types.
2 Feature extraction of acoustic signals 2.1 Kurtosis Kurtosis K represents the peak height at the average value of the probability density distribution curve. It is the normalized fourth-order central moment, which is defined as where, x i andx is the discrete value and mean value of signal x respectively, σ is the standard deviation. K is a dimensionless parameter and suitable for the early detection of structural defects.

Wavelet-based SK
The SK generated from time-frequency analysis and highorder statistics was designed for detecting the transient components in a signal. Wavelet transform has been widely used in transient analysis of vibration signal for its uniform resolution on logarithmic frequency scale. Among the available wavelet functions, the complex Morlet wavelet is thought to be the most suitable for the purpose. Assuming y(t) is the response of non-stationary signal x(t), its frequency domain expression decomposed by Wold-Cramer stochastic process can be described as follows where, H(t, f) is the complex envelope of y(t) at frequency f, and dx(f) is an orthonormal spectral process related to input x(t). So the SK of y(t) can be defined as the normalized fourth-order cumulant where, <⋅> is the time-average operator. When the complex Morlet wavelet was used to calculate the SK, the modulus square of the wavelet coefficients resembles the series of power spectrum H 2 (t, f) of Eq. (3).

AR-MED based kurtosis extraction 2.3.1 AR model
Assuming that x k is the zero-mean stationary signal sequence, the pth-order AR model of x k can be defined as follows: where, φ i are the AR model coefficients, and ε k is the residual error which represents the difference between the actual signal and the AR prediction. For non-stationary signal, the residual error is a mixture of noise and transient components of signal. The order selection of the AR model is the key. The lower an order is, the smoother the estimated spectral peaks are, which causes a significant prediction deviation, however, the higher the order is, the more intense the estimated spectrum fluctuation are, which leads to pseudo spectrum peaks and statistical instability (Wang and Wong 2000). In this paper, the Akaike Information Criterion (AIC) is used to determine the optimal order of the AR model. When AIC gets the minimum, the optimal order of the model can be obtained (Bengtsson and Cavanaugh 2006).

Minimum entropy deconvolution (MED) technique
Assuming that ε k is an acoustic sequence filtered out stationary components, according to the theory of signal and system it can be described as where, g k is the defect signal component, n k is the noise, h k is the impulse response function of the transmission path, and symbol * represents convolution. The basic purpose of MED is to seek a filter that eliminates the influence of transmission path and noise and make the output approximate the defect signal by solving an optimum set of the filter coefficients. The realization process of MED filter is as follows: 1. Assuming L, ε k and y k is the length, input and output of a filter f k respectively, then 2. Calculating the kurtosis of the output y k 3. Solving the filter coefficients by maximizing the kurtosis K f , i.e., taking the first partial derivative with respect to f k for Eq. (7) and set to zero 5. Substituting Eqs. (6) and (9) into Eq. (8), the Eq. (8) can be rewritten as follows: Eq. (10) can be written in matrix form: c Tf (11) where, the column vector c is a weighted sum of crosscorrelation of the input and output of filter f k , T is the Toeplitz autocorrelation matrix of the input, and f is the coefficient vector of the filter f k .
6. Using an iterative algorithm for solving the filter coefficients of the Eq. (11) as follows: (a) Calculating the autocorrelation matrix T.
The iteration ends if the expected value E(δe) is less than the set tolerance limit. Otherwise, update the filter coefficients f (i) =f (i + 1) and repeat the next iteration loop from step (c).

Feature parameter extraction based on AR-MED
When the cross-section of a log is subject to unit longitudinal impact force, the acoustic signal generated reflects the internal quality of the log. If a log is in homogeneous and continuous condition, the response signal can be represented as the superposition of an infinite number of principal vibration modes according to the vibration theory, whereas if there exist internal defects (such as void, decay, crack, etc.) in the log, mutation, distortion or attenuation components caused by defects inevitably skew or modulate the response signal. Therefore, the acoustic signal x k collected from a log can be expressed as the convolution between the response function of transmission path h k with the combined signal of periodic components d k (deterministic signal), defect components g k and background noise n k , namely Figure 1 illustrates AR-MED based extraction process of feature parameters. The steps are as follows: (1) Predicting and subtracting the deterministic components using the AR model from the response signal to obtain the residual signal containing the defect components.
(2) Using the residual signal as the input of theMED filter, and making the filtered output signals approximate to the defect components through multiple iterative operations. (3) Computing the kurtosis of the output signal, and applying them as the global parameter for quality evaluation or internal defect prediction of logs. (4) Performing complex Morlet wavelet transform on the output signal and computing the SK, and then using the SK for approximate identification of the internal defect types of logs.

Log samples
Fifteen hardwood logs, including two black cherry (Prunus serotina Ehrh.), five white oak (Quercus alba L.), six red oak (Quercus rubra L.), and two cottonwood (Populus deltoides Bartr. ex Marsh.), were obtained from a local wood producer in Madison, Wisconsin for this study. The selection of log specimens was based on visual appearances so as to include both sound and unsound logs. The unsound logs had a range of defects (such as voids, cracks, ring shake, and internal decay) and were in different deterioration levels (Xu et al. 2019).

Experimental procedures
All logs were first assigned log numbers and their basic dimensions were measured. Then a longitudinal acoustic impact test was performed on each log by impacting at its one end using a medium-sized hammer. The acoustic response signals were acquired through two sensor probes connected to a data acquisition card inserted into a laptop, with a sampling frequency of 20 kHz, and the sampling length of 2000 points. Following the impact test, a Fakopp Microsecond Timer was used to conduct a series of radial stress wave transmission tests on each log. Stress wave transmission time measured was applied to evaluate the physical condition at the tested locations of each log to locate the cutting positions as so to disclose further the internal defects. Log defect ratio, defined as the percentage of the deteriorated wood in each log, was calculated by dividing the volume of defect areas by the log volume. i.e., based on the measurements of defect areas on dissected cross sections and their extended lengths, a ground-truth defect map was created for each log, from which dimensions of all the defects and the log could be obtained to calculate their volumes respectively. Detailed experimental conditions and procedures can be found in Xu et al. (2019). The species and physical characteristics of the logs are listed in Table 1.

Results and discussion
The moisture content (MC) of the log samples ranged from 31.3 to 95.6% (Table 1), higher than the fiber saturation point (30% MC). As a result, the logs were considered in green condition, and the acoustic parameters discussed in this paper were regarded as "green wood" parameters.

AR-MED based kurtosis extraction
The kurtosis for all hardwood logs are summarized in Table 2. To demonstrate the effectiveness of the AR-MED algorithm in enhancing the mutation signal, the kurtosis of the residual signal filtered by the AR model and that of the output filtered by the AR-MED were compared using acoustic signal of log no.14 (red oak) as an example ( Figure 2). The response signal of log no. 14 has almost regular variation of waveform with no clutter as shown in Figure 2a, and its corresponding spectrum (Figure 2b) does not appear to contain any anomaly frequencies except for the fundamental frequency and the corresponding higherorder multiple frequency, so it is difficult to determine whether defects exist in the log or not. But from the dissected results, there do exist some defects such as relatively severe decay, void, ring shakes, etc. in the log ( Figure 3) and the defect ratio calculated is about 10% according to the ground-true map (details can be found in Xu et al. 2019). Comparing to the main vibration components, mutation components caused by defects are relatively weak and almost submerged in the dominant components. AR predictive filtering was performed on the signal to remove the stationary components and obtain the residual signal containing defect components. The distribution of the order of the AR model with respect to the AIC values and the residual signal are shown in Figures 2c, d, respectively. The order of 400 was the optimal order for the AR model that separated the stationary and mutation components. However, the AR filter could not isolate noise from mutation components due to the minimum phase condition assumed in the AR model. The kurtosis of 4.64 was small and could not reflect the extent of log deterioration. To suppress noise and enhance defect-induced mutation components, the residual signal was processed further by a MED filter, and the kurtosis of the output reached 18.00 (Figure 2e). Observing its spectrum shown in Figure 2f, the mutation components were separated effectively from the response signal and enhanced significantly.

Relationship between defect ratio and kurtosis
The kurtosis of 15 hardwood logs ranged from 10.8 to 27.5, with a coefficient of variation (COV) of 26.6%. Although a meaningful within species statistical comparison cannot be obtained from the small sample size of the logs in each species, the trend of decreasing kurtosis with the increase of defect ratio is evident. The reason could be attributed to the fact that defects of the logs were mainly rot and void. When these defects are present in logs, the occurrence probability of relatively small signal amplitudes increases, and thus the kurtosis decrease accordingly, and the more severe the defects are, the smaller the kurtosis values are (Mucciardi et al. 2011). Figure 4 shows the data plot of defect ratio and kurtosis of the hardwood logs. An excellent negative relationship was observed between log defect ratio and kurtosis, with the coefficient of determination (R 2 ) of 0.89 for a second-order polynomial regression, which is a significant improvement comparing with the R 2 of 0.17 for the regression result between defect ratio and acoustic velocity (Xu et al. 2019). If the difference of the MC was ignored, the distribution of the kurtosis with respect to defect ratio seemed to be less affected by log species compared with acoustic velocity. As a result, the kurtosis extracted from defect components isolated could represent internal condition of a log more effectively than the acoustic velocity directly derived from the original signal.

SK based defect type identification 4.2.1 Complex Morlet wavelet-based SK extraction
The SK expands the concept of the kurtosis, as a global parameter, to that of a frequency function, revealing how the mutation characteristics of a signal are distributed in the frequency domain. Therefore, the SK is a powerful tool for detecting the presence of transients in a signal. Even if the transients were hidden in deterministic components and noise, the SK can identify them by indicating which FBs they appear in.    Since the SK is related to the central frequency and bandwidth of wavelet, the Kurtogram (Antonia and Randall 2006), a fourth-order spectral analysis tool, was first made to select the optimal central frequency (for simplicity and real-time requirement, the bandwidth is fixed of 300 Hz in this paper). A variety of defect types coexist in a log, but in this study, defect types identified were mainly decay, void, crack, or other structural defects. Also, because large SK values located in some FBs were caused by the dominant components (associated with defects), while small SK values were induced by noise, and thus the small SK can be omitted for more effective analysis (Antoni and Randall 2006;De la Rosa et al. 2015). The FBs with the larger SK values were extracted and contrasted to the actual defects from the cutting-open log one by one to determine the corresponding relationship between them.
The complex Morlet wavelet transform was performed further on the AR-MED based output signal of log no.14 to calculate the SK. The distribution map of the SK (Kurtogram) of log no.14 is shown in Figure 5. Obviously, the SK values calculated are different for wavelets with different center frequencies (namely 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, and 5) for the signal decomposition. The maximum of the SK is used as an indicator for selecting the optimal center frequency. For this signal, when the central frequency of wavelet was of 2, the maximum of SK of 19.68, corresponding to the FB of 1850-2150 Hz, can be obtained, so the number of two is the optimum frequency of wavelet for the signal decomposition, that is, the SK values were obtained based on the optimal wavelet. Because the internal defects of the hardwood log were not single, other FBs with some SK values larger than 40% of the maximum should also be the characteristic intervals of useful signal components. The SK values were sorted and then FBs corresponding to the large values were selected for analysis. It is assumed that the larger the SK value of a certain FB is, the larger the proportion of material components corresponding to the FB in the log is. In addition, some research showed that if the internal defects of a log were dominated by decay or void, the frequency of the acoustic signal transmitted in the log is relatively low (Dackermann et al. 2014;Dunlop 1983;Mucciardi et al. 2011); and if the defects were mainly cracks or splits, the high frequency components would appear in the acoustic signal (Hsiao et al. 2008). Therefore, it is speculated that there exist some defects such as decay, void, etc. in the log if larger SK values were located in low FBs of a signal. On the contrary, defects such as cracks or splits do exist if larger SK values were found in high FBs of the signal. As seen in Figure 5, the largest three SK values were located in FBs of 1850-2150, 1550-1850, and 2450, accounting for about 53.7% of the sum of the SK values. Still, the remaining values were significantly smaller, and thus they could be ignored in analysis. Observing the internal condition of log no. 14 ( Figure 3), it was found that internal rot, void, and ring  shake were the dominant defects, other defects, such as knots, etc., were very small. The predicted defects basically conformed to the actual internal condition of log no. 14, which indicated that the conjecture mentioned above between distribution of the FBs and the corresponding defects should be correct. Therefore, the corresponding relation between the FBs of the SK and defect types could be assumed as listed in Table 3 from which the main types and ranking of severity of defects could be roughly estimated.
According to the above method, the SK values of all 15 hardwood logs were extracted, and the Kurtogram was drawn together and shown in Figure 6. The large SK values in Figure 6 are mainly distributed in the FBs of 650-4550 Hz. To determine the defect types and overcome the observation inconvenience caused by extensively large difference among the SK values of each log, the SK values were normalized by dividing SK located in each FB by the sum of the SK values of all FBs, and in the meanwhile the frequency range of Kurtogram was focused on the interval of 650-4550 Hz. Then the Kurtogram was redrawn in Figure 7 with the smaller SK values being set to zero and only the larger ones kept.

SK based defect identification
The SK of log no.1 was mainly concentrated in the frequency range of 3350-4550 Hz, and the sum of SK from four FBs accounted for about 64.7% of that of all SK values. According to the hypothesis of the corresponding relation between FBs and defects listed in Table 3, there should be no internal defects in the log except end checks, as evidenced by observing the actual sawing results (Figure 8a). For log no.3, there were no apparent differences among the SK values from FBs of 650-3650 Hz, with a COV of 18.2%. It Figure 5: Kurtogram for selecting the optimal center frequency of complex Morlet wavelet (log no. 14-red oak).   was found that there existed a variety of defects, such as void, decay, ring shake, etc., in the log, and the defects were all significant, without obvious primary and secondary by examining internal condition of the log (Figure 8b). The observed results further validated the speculation in Table 3. The major SK distribution of log no. 4 was located in frequency ranges of 1250-1850and 2750-3650 Hz, accounting for about 62.8% of the total SK values, and the SK of each frequency range was not significantly different. It could be deduced from the Table 3 there existed some defects such as rot, voids, and cracks in the log. The basic consistence between the predicted defects with the actual sawing results (Figure 8c) further verified the correctness of the assumptions in Table 3. Similar to the analysis of the above examples, the feature FBs and predicted defects from the remaining logs were listed in Table 2, no more tautology here. By comparing Table 2 with Table 1, the types and primary and secondary of the main defects predicted were roughly consistent with the sawing results, which confirmed the validity of the relation between FBs and defects presented in Table 3. It should be noticed that log no.8 (cottonwood) was a particular case, i.e., global parameters of the longitudinal acoustic signal, whether kurtosis or acoustic velocity, time centroid and damping ratio presented in Xu et al. 2019, implied that the log could be a sound log, but actually one with abnormal internal condition (Xu et al. 2019). Further examining the distribution of the larger SK of the log, it was found that they were located in FBs of 2750-3050 and 1250-1850 Hz, respectively. According to the assumption presented in Table 3, the log should contain the defects such as cracks and decay. It was demonstrated that these defects did exist in the log by checking the sawing results (Figure 8d). The research results further indicated that the SK could determine the existence and type of major defects which not be detected even by the global acoustic parameters.

Conclusions
With respect to the non-stationary and characteristic overlap of the response signals resulting from the acoustic impact testing of hardwood logs, a method of kurtosis extraction based on AR-MED was proposed to evaluate the internal quality of hardwood logs in terms of defect ratios and kurtosis. To further identify the types of defects, a complex wavelet-based SK method was presented. Based on the results of this study, the following conclusions were drawn: (1) The AR-MED process could effectively separate the defect components from the main vibration signals and significantly enhance the kurtosis of defect components, which provided a solid basis for evaluating the internal condition of hardwood logs using the kurtosis values.
(2) Compared with acoustic velocity, kurtosis could offer more accurate assessment of the internal condition of hardwood logs. A strong quadratic regression relationship existed between the log defect ratio and the kurtosis with a coefficient of determination of 0.89. (3) Compared with some global acoustic parameters, the SK can not only determine the types and severity of major defects, but also can identify the defects that are not detectable to global acoustic parameters.
The hardwood logs used in this study contained a wide range of defects that were difficult to quantify with precision. The response signals from impact testing of the logs were complex and variable. Although the results showed a great potential for implementing the acoustic impact testing method combined with advanced signal analysis to predict the types and extent of defects in hardwood logs, large trials with mill-length logs are needed to test the applicability of the method presented in the paper.