The application of carbon fiber reinforced composites (CFRP) is rapidly increasing due to its high specific strength and stiffness. Although CFRP components are either autoclaved or 3D printed to a near net shape, secondary machining is often required for edge trimming and drilling. Machining of CFRP is challenging due to the difference in interaction behavior of the individual constituents with the cutting tool, high strength of carbon fiber, poor thermal conductivity of matrix material and low interlaminar shear strength [1] . In addition, the material removal process is dependent on the fiber orientation [2] , which governs the damage induced during machining process. This damage is in the form of fiber pullouts, delamination, thermal distortion, swelling and splintering. This contributes towards excessive tool wear and prove to be uneconomical, deter the production rate and impair the surface integrity [3] . The sensitivity of CFRP to the machining parameters necessitate monitoring of the machining process to ensure tight tolerances and superior surface quality. Several authors have used different techniques to monitor the milling process using force and acoustic emission signals for surface and tool condition monitoring [4] [5] [6] [7] . It has been well established that the machining quality and surface damage can be best correlated with simultaneous time-frequency analysis of the force signals. Li et. al [8] studied the force monitoring of UD CFRP using discrete wavelet analysis with Daubechies (db3) wavelet. Sudden peaks in detailed coefficient was correlated with fiber debonding (at high fiber angles). Plaza et. al [9] implemented real time monitoring of surface quality in CNC turning of CFRP. Wavelet packet transform was utilized and identified parameters to monitor average roughness Ra at low computation cost. Recently, Pahuja et al. [10] implemented continuous wavelet analysis to find the effect of fiber orientation. They also used wavelet packet transform in monitoring abrasive water jet machining of CFRP-Titanium stacks using acoustic emission signals [11] . Wavelet Packet transform is applicable to milling of CFRP since the milling process is non-stationary in nature, especially for CFRP where the presence of different media results in varying burst-like signals, depending on the fiber direction. Wavelet packet transform (WPT) is known to be superior to continuous wavelet transform (CWT) which is a redundant transform resulting in high computation time. Discrete wavelet transform (DWT) on the other hand only decomposes detailed coefficients at each level. Knowing the potential of WPT in process monitoring, a dearth of literature motivated this study. No effort has been made in identifying optimal WPT parameters in monitoring surface quality (and machining induced damage). In this study, an algorithm is proposed where wavelet parameters (mother wavelet and decomposition level), signal features, and type of force signal was identified. This resulted in high correlation of those identified features with process parameters (feed and speed), and surface quality (roughness Rz) at optimal parameters with R 2 >91%. Overall, this approach can be effectively used to monitor the surface roughness and resulting machining induced damage in milling of CFRP.
The material used in this study is (UD) carbon reinforced plastic (CFRP). A 500x50x6.35 mm unidirectional Graphite Epoxy laminates was used. The laminates were autoclave cured, and manufactured from prepreg unidirectional tape. The fiber used to fabricate CFRP is T800H high strength, intermediate modulus yarn, and toughened epoxy resin 3900-2. The fiber diameter was 5 µm with volume fraction about 60% and average ply thickness of 190 µm. Four different fiber orientations (with respect to tool feed directions) were used: 0°, 45°, 90° and 135°. For a given laminate, all laminae had the same fiber orientation. These panels of different fiber orientations were machined out into rectangular specimens using abrasive water jet.
The edge trimming experiments were conducted on aforementioned materials, using 3 axis Haas CNC milling machine. Upmilling was employed using a 2 flute carbide tool with diameter and flute helix angle of 6.35 mm and 30° respectively. The radial depth of cut was 1 mm with a through-the-thickness edge trimming operation. All the experiments resulted in edge trimming with 50 mm long straight cutting path. Fig. 1 shows the experimental set up and data acquisition system. The experimental conditions and process parameters are tabulated in Table 1 . Angle between fiber axis and tool feed direction A rotating type dynamometer (Kistler 5223B) was used to record the cutting forces F x , F y , F z and M z (moment) at 7999 Hz sampling rate. The surface was inspected using optical microscopy, Scanning electron microscopy (SEM) and 2D profilometry. The instantaneous angle between the cutting edge and fiber direction varies with each interaction of the cutting edge and with the height due to helical shape of the cutting edge. Since this instantaneous interaction angle govern the damage mechanism and resulting cutting forces [12] , the 2D roughness profiles were measured at an interval of 1 mm starting from the top edge. Fig.  2 illustrates the roughness measurement on the machined surface. The roughness profile was measured using a contact type (MarSurf XR 20) 2D
roughness measurement system. A skidless contact type probe of diameter 2 μm was used to record measurement of 5.6 mm evaluation length 0.8 mm cut- Roughness profi les off length. The resolution in vertical and horizontal direction was 7 nm and 0.5 μm respectively. Roughness parameters-Average roughness (R a ), tenpoint average (Rz) and maximum profile height (Rt) were calculated upon gaussian filtering as per ISO 16610 standard. FEI XL30 Sirion FEG microscope was used for SEM observation. Refer to [10] for more details about the experimental setup and measurement equipment.
Wavelet packet tranform (WPT) is a simuntaneous time-frequency tranform in which a mother wavelet function is used to decompose a signal into two subbands of approximation and detailed coefficents. Successive decompsosition(s) is applied by translating and scaling the wavelet function, on both approximation and detailed subbands [11] . In this study, four different signals were captured -Fx, F y , F z and M z during the milling process, as shown in Fig. 3 . Due to the excessive data volume, only 1.5 seconds of force signal for each condition was analyzed. Wavelet packet transform was calculated for the following force vectors: F x , F y , F z , M z and � � � � � � � � � � � . For the sake of brevity, only the resultant force F a is discussed in this article.
A set of 35 different mother wavelets (belonging to four different classes -Daubechies, Coiflet, Symlet and Biorthogonal) was used. The coefficients of wavelet packets (i) were generated at each decomposition level (j) with maximum decomposition level of Level 10. A total of 2 j wavelet packets were resulted at j th decomposition level, with a maximum of 1022 packets at Level 10. Selected features were extracted from each packet at the subsequent steps. A set of 7 different features was explored and calculated for each i th wavelet packet at j th decomposition level. These features and their definition is tabulated in Table within Fig. 3 . These features include arithmetic average, standard deviation, skewness, kurtosis, energy, Shannon entropy and energy-entropy coefficient calculated for each i th packet with K coefficients (K=N/2 j ) where N is the total number of Fig. 3 . Wavelet packet analysis algorithm data points in the analyzed signal. The definition of Mean, standard deviation, skewness and kurtosis as same as conventionally used in statistics. Energy of each wavelet packet at the given decomposition level is given by equation 1 where is C k the k th coefficient.
The total energy at each decomposition level is given by the summation of energy of all packets at that level, as given by equation 2.
Next, the probability function was determined which is the ratio of Energy at each packet to the energy at all packets at that decomposition level. This is given by equation 3.
This probability is further used to calculate Shannon entropy, given by equation 4 and the total Shannon entropy of all the packets at the given level is given by
Further, the ratio of packet energy and total Shannon entropy was determined using equation 4 and 5, as given by equation 6 .
The calculated features were then plotted against the frequency ordered wavelet number, and dominant peaks were identified. The summation of the absolute values of these peaks identified as a unique classifier at a given decomposition level. This WPT indicator is given by equation 7. � � ∑ η ����� (7) Further, regression models using ANOVA were developed between the signal features, and process and workpiece parameter(s). Each regression was a 4 th degree polynomial with insignificant terms (p>0.05) were removed using backward elimination (α=0.05).
A regression model of 39 experimental runs was modeled for each of the 5 force components and each of the 35 wavelets, with 7 features calculated using all the terminal wavelet packets at a decomposition level among 10 levels. In total 12250 regression models were developed and analyzed for coefficient of determination R 2 , adjusted R 2 and RMSE (root mean square error). The mother wavelet, decomposition level and signal feature was selected which resulted in maximum R 2 , adjusted R 2 and minimum RMSE. Similar methodology was followed to correlate the surface quality parameter R z with signal features, and a total of 2450 models were analyzed. Once the optimum parameters for wavelet packet decomposition were identified, the signal could be effectively used in process and quality characterization, and quality control by monitoring the signal features that best represent the surface quality. Fig. 3 shows the flowchart for the aforementioned algorithm.
The maximum Arithmetic mean Ra, ten point average R z and maximum peak-to-valley height R t observed was 23.46 µm, 104.79 µm and 131.37 µm respectively. Fig. 4 shows a positive trend between feed and roughness. In general, high feed and low spindle speed resulted in rougher surfaces. The roughness increased with fiber orientation angle from 0° to 90°. However, at 135° fiber angle, the surfaces were extremely rough with excessive damage such as pullouts and fiber spring back. Fig. 5 shows the surface roughness profiles observed in this study for 5 mm/s and 25 mm/s feed rate conditions. It is worthwhile to note that unlike homogeneous materials, the surface roughness was inconsistent when measured at different locations on the surface generated with a given set of process conditions. This is due to (1) The cutting edge helical geometry that changes the instantaneous interaction angle, (2) The damage is different near the top ply which is unsupported. This variance in roughness measurement was calculated. At 135°, the specimens had the highest variance for Rz (1153.3 µm), while 90° specimens had variance of 95.6 µm. The 0° and 45° specimens showed least variance (4.6-8.32 µm). It can be concluded that 135° and 90° specimens were more sensitive to the machining conditions. To better represent the severity of surface anomalies, the maximum Rz value was considered henceforth. Fig. 6 shows the skewness vs kurtosis plots of surface roughness profiles for specimens machined at different fiber orientations at 6000 rpm. As can be seen, most of the profiles are leptokutuic (kurtosis >3). At θ=0° fiber orientation, the skewness ranged between -1 and 1 and kurtosis between 2 and 6. At θ=45°, skewness for 1000 rpm was slightly positive but mostly concentrated around zero for 3000 and 6000 rpm. Kurtosis was low with a few exceptional cases where it went up to 6. At θ=90° and θ=135° highly negatively skewed profiles were observed with highest kurtosis among all fiber orientation angles. Kurtosis up to 10 was observed with an exceptional case of 17.64 at θ=90°. A few other exceptions were θ=0° at 3000 rpm and 15 and 25 mm/s feed. A high kurtosis was observed at two profiles near the top edge which is likely to have captured the damaged topography due to pullouts and sprung back fibers at the entry side.
This section discusses the correlation between the signal features and process parameters using wavelet packet transform. Fig. 7 shows the R 2 values of regression models of wavelet packets computed using different mother wavelets and decomposition levels for different signal features as mentioned in Fig. 3 except feature -skewness (f3) due to poor regression (R 2 <50%). In general, regression models for most of the features at lower levels were missing, especially at L2. This is due to the lack of distinct peaks detected at these levels. However, at higher levels of decomposition, the signal content is distributed among several wavelet packets. Distinct peaks of coefficient features were identified at this high frequency resolution.
The standard deviation (f2) of the wavelet packets is preferable over the arithmetic mean (f1) as the features correlated better at lower level of decomposition. This means that the correlation between signal features and process parameters can be established at low computation power. This is an important criteria to select the wavelet transform parameters. For standard deviation, among the Daubechies, db9 was the mother wavelet with highest R 2 and adjusted R 2 of 97.67% and 96.31% respectively. The RMSE was 5.97. A few Biorthogonal -Bior 4.4, Bior 5.5 and Bior 6.8 regressed well with 97.76%, 95.89% and 97.56% respectively. The Adjusted R 2 was greater than 93.2% and RMSE less than 8.95. Coiflet 4 and Symlet 5-8 were equally competitive with R 2 greater than 97.14% and RMSE under 6.32. The features skewness and kurtosis lacked dominant peaks while calculating the indicator for each WPT. For kurtosis, Coiflet 4 at Level 5 was reasonable, with R 2 = 81.35%. Among Energy (f5), Entropy (f6) and Energy-Entropy coefficient (f7, R 2 was highest for Entropy, followed by Energy and Energy-Entropy coefficient. For Shannon Entropy, level 4 correlated the signal and process parameters better for Biorthogonal wavelet. Daubechies (db2-db4) and Symlets were second in performance. Using Energy -Entropy coefficient, Bior 5.5 and Bior6.8 at Level 8 was optimum with 89.05% and 89.85% respectively. Fig. 8 shows the results of regression model between process parameters and feature f2 (standard deviation) calculated with Bior 5.5 mother wavelet and decomposition level 8. The R 2 , adjusted R 2 and predicted R 2 for this model was 97.77%, 96.61% and 94.74% respectively. At 6000 rpm (Fig. 8(c) ), a positive correlation between feed and f2 was observed. This means high standard deviation of dominant wavelet coefficients at high feed. This trend was also observed with increasing fiber orientation up to θ=90°. Beyond 90°, a drop in standard deviation was observed. At 3000 rpm (Fig. 8(b) ), a more prominent increase in f2 was observed. The trend between f2 and feed, f2 and fiber orientation was positive. A 286.4% increment from θ=0° to θ=135° was observed at 5 mm/s, which went up to 562.88% at 25mm/s feed rate.
The trend was most accentuated at 1000 rpm, where 187.7% increment was observed at 5 mm/s while 485.1% at 25 mm/s. Overall, the feature f2 increased with increasing resistance to cutting, which also corresponds to higher resultant cutting force.
This section discusses the correlation between surface quality (maximum Rz) and signal indicators calculated using aforementioned features, wavelet type and decomposition level for five different cutting force signals. The regression between signal features and surface quality resulted in poor predictive models. However, when supplemented with the information about workpiece parameter-fiber orientation, the prediction was sufficient. Fig. 9 shows the R 2 values of regression between surface roughness, fiber orientation and signal features extracted using WPT. As evident, indicator f1 (Arithmetic mean) and f2 (Standard deviation) resulted in more acceptable results than skewness and kurtosis. The arithmetic mean resulted in higher R 2 at decomposition level 8, 9 and 10. The maximum R 2 was observed with db3 Daubechies wavelet at Level 8 with R 2 and adjusted R 2 of 97.09% and 95.09% respectively. Symlet 4-6 at Level 9 was also somewhat satisfactory. For standard deviation (f2), Bior 5.5 resulted in R 2 and adjusted R 2 of 95.4% and 93.5% respectively at decomposition Level 8. Next, Wavelet Packet energy indicator resulted in consistent results from level 8 onwards. As shown in Fig. 9 , R 2 at Level 8 was 94.8% calculated using Bior5.5 wavelet. Except for the Biorthogonal 5.5, none of the Biorthogonal wavelets correlated well. 
As evident from Fig. 10(a) and (b), both f2 and f7 signal features resulted in positive trend with surface quality. The minimum roughness was observed at θ=0° and θ=45° specimens. For feature f2, about 453% increment in R z was observed when f2 was spanned within the extremes of the experimental range. When compared to θ=0° orientation, the predicted R z for θ=45° was about 48.38% lower, 43.45% higher for θ=90° and 147.46% higher for θ=135° when compared to θ=0° orientation. Fig. 10 (c) and (d) shows the predicted vs experimental results for f2 and f7 indicators. The predicted R 2 for the model represents the variation in new data explained by the model. The predicted R 2 for f2 and f7 was 88.4% and 58.46%, which makes f2 a better candidate for predicting the surface quality.
In this study, process monitoring of CFRP edge trimming was conducted using wavelet packet transform of force signals. Wavelet Packet transform (WPT) was used to identify the key characteristics of the signals. A novel method was proposed in which a signal is decomposed using wavelet transform for up to a certain decomposition level using a mother wavelet. About 35 mother wavelets and up to 10 decomposition levels were used. The coefficients of the terminal wavelet packets at each decomposition level were used to calculate a spectrum of packet features -arithmetic mean, standard deviation, energy, entropy, and energy-entropy coefficient. Regression models were developed and optimum WPT parameters were identified to better correlate the process parameters (speed, and feed), workpiece parameters (fiber orientation) and signal features. Similar methodology was adopted to correlate the signal features, and surface roughness, and made the roughness predictive model independent of the process parameters. Following key conclusions can be drawn from this study.
1. The force captured in feed direction Fy and resultant cutting force F a were two the force signals that provided satisfactory results. The resultant force was used in this study as it captures the effect of force in z-direction as well. 2. Feature f2 (standard deviation) and f7 (Energy-Entropy) coefficient were found to be the key characteristic features, which can effectively serve as signal classifiers. The machining induced damage as well as the signals generating the anomalous surface topography can be related with high standard deviation, high energy, and low entropy. 3. Among several wavelets and decomposition levels, Biorthogonal 5.5 and decomposition level 8 was found to be optimum in correlating process parameters, workpiece parameters and signal features. 4. The correlation between feed, speed, fiber orientation and the signal feature f2 was determined with coefficient of determination 97.8%. 5. The signal features were insufficient to independently predict the surface quality. When combined with the workpiece parameter (fiber orientation), f2 and f7 was found with coefficient of determination >91%. In practical, this information can be deduced from the cutting geometry, known priori such as in applications like contour trimming. An effective control mechanism can be derived by monitoring the signals.
A thresholding technique can also be employed for predicting the surface quality while machining. For example when using the feature f2, if Rz<20 µm is desired, the threshold values (with 95% confidence interval) would be 3218.5, 5068.4, 2209 and 683.7 for 0°, 45°, 90° and 135° respectively. Similarly, the threshold values can be determined for feature f7. Feature f7 appears to be a saddle at R z <20 μm. The contour maps are shown in Fig. 11 . This work suggests the effectivity of WPT in quality monitoring of CFRP machining. However a different material and tooling conditions are likely to result in different threshold values. In future, the same framework would be extended to neural network and machine learning to adapt to the system specific environment.