# Accuracy Improvement of In-line Near-Infrared Spectroscopic Moisture Monitoring in a Fluidized Bed Drying Process

Andrey Bogomolov 1,2 \*, Joachim Mannhardt 1 and Oliver Heinzerling3

1 Blue Ocean Nova GmbH, Aalen, Germany, 2 Samara State Technical University, Samara, Russia, 3 Drug Product Development, AbbVie Deutschland GmbH & Co. KG, Ludwigshafen am Rhein, Germany

An exploratory analysis of a large representative dataset obtained in a fluidized bed drying process of a pharmaceutical powder has revealed a significant correlation of spectral intensity with granulate humidity in the whole studied range of 1091.8–2106.5 nm. This effect was explained by the dependence of powder refractive properties, and hence light penetration depth, on the water content. The phenomenon exhibited a close spectral similarity to the well-known stochastic variation of spectral intensities caused by the process turbulence (the so-called "scatter effect"). Therefore, any traditional scattercorrective preprocessing incidentally eliminates moisture-correlated variance from the data. To preserve this additional information for a more precise moisture calibration, a time-domain averaging of spectral variables has been suggested. Its application resulted in a distinct improvement of prediction accuracy, as compared to the scatter-corrected data. Further improvement of the model performance was achieved by the application of a dynamic focusing strategy when adjusting the model to a drying process stage. Probe fouling was shown to have a minor effect on prediction accuracy. The study resulted in a considerable reduction of the root-mean-square error of in-line moisture monitoring to 0.1%, which is close to the reference method's reproducibility and significantly better than previously reported results.

Keywords: fluidized bed drying, moisture monitoring, NIR spectroscopy, light scatter, scatter correction, lighthouse probe, process analytical technology

# INTRODUCTION

Fluidized bed drying is a common unit operation routinely performed in the pharmaceutical production of solid dosage forms. In a typical batch granulation process, the drying stage immediately follows either the fluidized bed or high-shear granulation stage. It is often considered as one of the most critical steps for achieving stable product quality, i.e., for obtaining granules with desired properties at their minimal variability. Therefore, a close monitoring of the residual moisture content in the process medium is necessary for any quality assurance system in granulate production.

In modern industrial practice, moisture is commonly analyzed in isolated samples. Karl Fischer titration is a classic water analysis technique that has been widely used for decades. A viable alternative accepted by pharmacopeias is thermogravimetric analysis with a drying balance that determines moisture content in the sample as percentage weight loss on drying (LOD). At present,

### Edited by:

Federico Marini, Università degli Studi di Roma La Sapienza, Italy

### Reviewed by:

Ludovic Duponchel, Université de Lille, France Huawen Wu, BaySpec, Inc., United States

> \*Correspondence: Andrey Bogomolov<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Analytical Chemistry, a section of the journal Frontiers in Chemistry

Received: 16 April 2018 Accepted: 10 August 2018 Published: 10 October 2018

### Citation:

Bogomolov A, Mannhardt J and Heinzerling O (2018) Accuracy Improvement of In-line Near-Infrared Spectroscopic Moisture Monitoring in a Fluidized Bed Drying Process. Front. Chem. 6:388. doi: 10.3389/fchem.2018.00388 both techniques are realized as compact desktop devices enabling the at-line analysis of samples taken from a running process.

For the process type studied here, the at-line analysis of granulate moisture content typically takes 20–30 min, representing a good alternative to off-line laboratory analysis of the final product. However, such operability is insufficient to carry out real-time process control, for example by generating alarms on abnormal process states and performing timely corrections. For the same reason, at-line analysis is hardly suitable for accurately determining the process end-point the time point at which the product reaches its optimal properties. Therefore, instant in-line monitoring of the moisture content in fluidized bed drying is strongly desired to provide a necessary level of process control and to meet growing quality requirements.

Near-infrared (NIR) spectroscopy is an undoubted favorite among real-time sensor systems for moisture monitoring in the production of solids, specifically, in the drying step (Roggo et al., 2007; Burggraeve et al., 2013; Da Silva et al., 2014). In such systems, the diffuse reflectance spectra of the process material are typically measured through an immersion probe. The key advantages of NIR spectroscopy as an in-line analytical technique include the suitability for measurements in media of highly variable bulk density, nondestructiveness, and the capability to place the probe into an appropriate position within the process space while keeping it connected to a remote spectrometer through a fiber optic cable.

The classic NIR spectroscopic moisture analysis relies on two intensive water absorption bands around 1,440 and 1,930 nm, enabling quantitative determination of the moisture in a wide concentration range. In low-selective NIR spectra, the component bands are essentially overlapped and their quantitative analysis requires the application of multivariate modeling, also known as chemometrics. In particular, the partial least-squares (PLS) regression algorithm (Sjöström et al., 1983) is widely accepted in process chemometrics (Bogomolov, 2011).

Over the last decades, the practical acceptance of NIR spectroscopy for in-line moisture monitoring in fluidized bed processing of powders and solids have been constantly growing. Published works (Frake et al., 1997; Rantanen et al., 2000; Zhou et al., 2003; Green et al., 2005; Nieuwmeyer et al., 2007; Skibsted et al., 2007; Luukkonen et al., 2008; Mantanus et al., 2009; Alcalà et al., 2010; Corredor et al., 2011; Peinado et al., 2011; Burggraeve et al., 2012; Demers et al., 2012; Möltgen et al., 2012; Obregón et al., 2013) have focused on the general feasibility of the analysis or on the investigation of specific experimental or modeling aspects (e.g., important process influences, sampling, control strategy, and model transfer). At the same time, the resulting models are typically built and validated on relatively small sets of samples and batches, which can be accounted for by the technical complexity of industrial experiments. Hence, the accuracy estimates reported for similar process setups and conditions are very diverse (Zhou et al., 2003; Green et al., 2005; Nieuwmeyer et al., 2007; Skibsted et al., 2007; Mantanus et al., 2009; Alcalà et al., 2010; Corredor et al., 2011; Peinado et al., 2011; Burggraeve et al., 2012; Demers et al., 2012; Möltgen et al., 2012) and the "ultimate" moisture determination accuracy by in-line NIR spectroscopy under widely variable process conditions remains unknown. Therefore, despite significant progress, the method can hardly be regarded as completely established yet.

In-depth considerations of NIR spectroscopic analysis in terms of light propagation in the complex fluidized bed process medium are rare (Rantanen et al., 2000; Luukkonen et al., 2008; Burggraeve et al., 2013). One of the main obstacles complicating the NIR spectroscopic monitoring of fluidized bed drying is related to process turbulence. A highly variable density of the material around the probe, and consequently the quantity of light reaching the detector, causes intensive random fluctuations of the overall intensity of in-line spectra that are often referred to as the "scatter effect." The problem is commonly resolved by preprocessing the spectra prior to the modeling step. The three most-used scatter correction methods are multiplicative scatter correction (MSC), standard normal variate (SNV), and spectral derivatives (Rinnan et al., 2009). The application of a scatter correction method to in-line process NIR spectra is ubiquitous; no exception has been found in the literature. In most cases, the choice of the preprocessing method is empirical or arbitrary.

In some publications, it was noticed that the NIR spectra expressed in the logarithmic reflectance units (lg(1/R)) exhibited a significant downward shift of the background as the drying progressed (Frake et al., 1997; Rantanen et al., 2000; Zhou et al., 2003; Luukkonen et al., 2008; Burggraeve et al., 2012). Two plausible explanations were suggested, both related to the altering of light scatter conditions in the course of drying. On one hand, the uniform decrease in spectral intensities could be caused by an increase in scattering particle size; this explanation was given by Burggraeve et al. (2012) and Frake et al. (1997). On the other hand, the presence of water on crystal surfaces affects the reflective properties of the granulated powder, resulting in a deeper light penetration and a subsequent higher absorbance of wetter samples (Rantanen et al., 2000; Luukkonen et al., 2008). Rantanen et al. (2000) provided an experimental evidence of the latter phenomenon by using the pharmaceutical excipient (microcrystalline cellulose) as well as inorganic glass beads ("ballotini") with a known size distribution.

The present work aims at building an accurate and robust functional prediction model for in-line moisture content monitoring in fluidized bed drying based on a large representative set of designed process data. Both experimental and modeling factors have been scrutinized to improve the performance of the prediction model. A thorough exploratory data analysis has been applied to help understand the process multivariate trajectory delivered by in-line diffuse-reflectance NIR spectroscopy better. In this study, we focus on efficiently using of the whole spectral information, including both absorption and scatter-related effects of water, to improve the performance of in-line moisture monitoring.

# MATERIALS AND METHODS

Twenty-five pilot-scale fluidized bed drying batches of a pharmaceutical powder mixture were studied by using a 256 pixel diode-array TIDAS 1121 SSG NIR spectrophotometer with

a wavelength range of 1091.8–2106.5 nm (J&M Analytik AG, Germany) that was equipped with the Lighthouse ProbeTM (LHP) from GEA Pharma Systems nv – Collette, Belgium (Engler et al., 2009) immersed into the process medium. The LHP was periodically cleaned and recalibrated without process interruption (see section S1.4 of **Supplementary Material**). The total number of cleaning cycles in all batches was 19.

The data of each batch included from 396 to 1,213 NIR spectra collected at 5-s intervals (16,303 spectra in total). In the course of the process, 301 samples of about 5 g (between 5 and 26 samples from each batch) were isolated and analyzed for moisture content as weight loss on drying using a HR73 halogen moisture analyzer (Mettler Toledo GmbH, Switzerland). Reproducibility checks for three LOD analyzers performed during the whole study showed that the measurement standard deviation error does not exceed 0.06% (section S1.2 of **Supplementary Material**).

The main process and the sample information are summarized in **Table S-1**. Out of the 301 samples, three were rejected from further analysis as evident outliers (section S2.3.1 of **Supplementary Material**).

Individual batch conditions were set in accordance with a developed experimental design to cover the whole range of practical process variability. Moisture content in the selected samples varied between 2.38 and 25.92%. The active pharmaceutical ingredient (API) was present in four assay levels: 0 (placebo), 0.1, 1.0, and 10.0 mg. The range of process temperatures was 30.5–49.7◦C. Eight batches (88 reference samples) formed a validation subset that was representative of the process conditions and used for model validation; the other 17 batches were used as the calibration set in that case (**Table S-1**).

A subset of 101 experimental samples were additionally analyzed off-line by using an MPA Fourier-transform (FT-) IR spectrometer (Bruker, Germany) with an integrating sphere (section S1.5 of **Supplementary Material**).

Principal component analysis (PCA) and PLS regression are multivariate data analysis algorithms described in the literature (Sjöström et al., 1983; Wold et al., 1987). The multivariate spaces, namely, PCA model principal components (PCs) and PLS latent variables (LVs) represented by their score (**t**) and loading (**p**) vectors, were used for exploratory data analysis. Conventional data preprocessing methods employed were MSC, SNV, and first-derivative using the Savitzky–Golay smoothing filter, as described by Rinnan et al. (2009).

Three validation techniques were applied with each regression model: leave-one(-sample)-out (LOO), a.k.a. full cross-validation (CV), leave-a-batch-out (LBO) CV, and validation by a preselected set (**Table S-1**). The performance of the models was characterized by rootmean-square errors (RMSE) of calibration, validation, and prediction, as well as corresponding determination coefficients R 2 .

A detailed description of data acquisition and analysis is given in section S1 of **Supplementary Material**.

# RESULTS AND DISCUSSION

# Exploratory Analysis of In-line Spectral Data

**Figure 1** presents a set of 1,213 in-line NIR spectra obtained in batch B03 (**Table S-1**). An expected intensity reduction of the main water band in the 1,920–1,940 nm range during the process is clearly observed. Another distinct feature is the high variability of spectral intensities over the whole wavelength range (the so-called "scatter effect"), caused by strong instant density fluctuations of the granulate (and its spatial distribution) around the probe.

At the same time, the overall spectral intensity tends to fall gradually during the process, generally following the dynamics of water reduction. This trend can be illustrated by the time dependencies of the spectral intensity at two separate wavelengths: 1932.0 nm at the maximum of the main water band and 1708.1 nm where no noticeable water absorption is expected. Both intensities strongly correlate with the reference moisture content (**Figure 2A**). Data smoothing along the time scale makes this correlation even more distinct.

The moisture- and time-dependent changes in the batch processes can be effectively visualized by using data animation (section S2.1 and **Video S-1**, **Supplementary Material**). Animated spectral data reveal the same trends, namely water band reduction and stochastic background variation accompanied by a gradual fall of the spectrum intensity in the whole range.

In this situation, preprocessing is desirable, but it should be applied to the data variable vectors, i.e., along the time scale, as shown in **Figure 2A**. As the turbulence effect is supposed to be pure noise, the smoothing of variables is a straightforward way to eliminate it with a minimal loss of the informative variance.

One of the simplest smoothing techniques, the moving window averaging algorithm, has been used to preprocess the matrix of spectral data **X**. In this method, each element xij in **X,** where i and j are respectively the object (spectrum) and variable (wavelength) indices, is replaced by a corrected value x s ij calculated as a mean of the surrounding points within a window having the width defined by an odd number k (Equation 1):

$$\mathbf{x}\_{ij}^{s} = \frac{\sum\_{i=i-(k-1)/2}^{i=i+(k-1)/2} \mathbf{x}\_{ij}}{k} \tag{1}$$

The transformation is performed for each variable in **X**. (k – 1)/2 end-points on each side of the variable vector were smoothed with a reduced window of (l – 1) . 2 + 1 points, where l is the point ordinal number from either spectrum end.

Data averaging within a selected time window is similar to a respective enhancement of the spectrum acquisition time, thus enlarging the virtual sample size captured by a single measurement. However, in contrast to the measurement time adjustment, the mathematical averaging does not place any limit on the time step of data acquisition, i.e., it can be performed with a time window that is much wider than the physical step size. A positive effect of the variable smoothing for the modeling of a fermentation process data has been reported (Skibsted et al., 2001).

Pair-wise correlations between the LOD values and the intensities at individual variables in the corresponding (closest to the sampling times) in-line spectra were analyzed in the whole wavelength range. **Figure 2B** presents linear correlation coefficients (r) as a function of wavelength in B03. All spectral variables exhibit a strong intensity correlation with the moisture content, even in the raw data. Eliminating the process noise using the suggested averaging method (Equation 1) results in a dramatic enhancement of r. It also looks natural that correlation maxima are observed around major water bands. However, even beyond the water absorbance regions, this correlation is very high. Thus, the lowest r observed in B03 at the short-wave end of the spectral region is still greater than 0.8 (**Figure 2B**); after the smoothing, this value increases to 0.98. Similar dependencies were observed for all the 25 studied batches.

A high correlation of lg(1/R) with the moisture content in the whole studied NIR range is in agreement with some published observations. This phenomenon can be explained by altering the refractive properties of the granulate (Rantanen et al., 2000). Indeed, in the course of drying, the liquid bridges holding the primary particles together (Burggraeve et al., 2013) are replaced by air. The crystal–air interface is characterized by a higher difference of refractive indices than the crystal–water pair. Thus, drying leads to a higher scatter and hence an increased quantity of diffusely reflected light reaching the detector—that corresponds to a decrease in the spectral intensity expressed in absorbance type of units. For relatively large particles constituting the granules, this effect should be wavelength-independent. An intuitive illustration of the particle wetting effect and its uncomplicated explanation using the representative layer theory was given by Dahm (2013). A similar correlation of the Raman spectral background with the moisture content was observed in our earlier studies on pellet coating (Bogomolov et al., 2010) and granulation process monitoring (Bogomolov, 2011), and was also explained by the effect of moisture on the light propagation conditions in the process medium. Considering the strength of the spectrum variable correlation with the moisture content observed in the whole range of process conditions studied, an earlier explanation of the phenomenon in terms of changing particle size distribution during the drying course (Frake et al., 1997; Burggraeve et al., 2012) has not been confirmed. This hypothesis does not agree with the complex shape of the correlation curve in **Figure 2B**. Particle size distribution can be a minor watercorrelated factor affecting the spectra of the drying process, though.

The effect of humidity on the light penetration depth in porous materials can be compared to the watermark technique commonly used for banknote authentication. The very name of watermarks comes from the visual similarity of paper thickness variation and its wetting effects, both resulting in a decrease in the back-scattered light. Darkening of wetted powders (e.g., sand) is another manifestation of the same phenomenon that is not limited to the visible light and should be inherent in any material with a highly developed surface. The spectral variance related to the changing refractive properties of the powder is also expected to be present in the in-line process spectra. However, being wavelength-independent, the moisture-related spectral changes are masked by the stochastic "scatter effect" and then eliminated by any scatter correction. Earlier studies on inline moisture analysis by using NIR spectroscopy neither paid any significant attention to the analytical information hidden in the "watermarks" nor attempted to use it in the modeling.

A deeper insight into the data structure and its modification by adopting different preprocessing methods was obtained by the PCA of augmented process data (section S2.2 of **Supplementary Material**) that makes possible the investigation of process trajectories of individual batches in the same multivariate factor space.

As one can see from the scores of batch B10 taken as an example here (**Figure 3** and **Figure S-3**), the first PC (95.49% of **X**-variance) of the raw-data model (**Figure 3A**) is strongly associated with the moisture content, while PC2 (4.23%) basically describes the process turbulence. A remarkable similarity of the first two loadings (**Figure S-4a**) with the correlation coefficient r = 0.998 is a confirmation of a close spectral affinity of these two phenomena. A scatter-driven correlation of spectral intensities with the moisture content is confirmed by the uniformly positive **p**1. A simultaneous presence of the water absorption peaks in this plot implies that PC1 tends to capture the whole variance due to the moisture reduction, related to both absorbance and scatter phenomena.

Although the process noise is basically described by PC2, it strongly pollutes PC1 and all further components in the raw-data model. The suggested smoothing method effectively eliminates this noise from the model scores (**Figures 3B,C** and **Figures S-3b,c**) without any essential change to the loadings (**Figures S-4b,c**). In contrast, the SNV, MSC, and first derivative (**Figures S-4d-f)**) strongly modify the whole factor space; they essentially remove random fluctuations from the first two score vectors (less noisy for the first derivative) but further PCs stay very noisy (**Figures 3D–F**). The smoothed data is suitable for exploring the process trajectories in the PCA factor space. Most of the minor features revealed in the refined scores **t**2-**t**7 (**Figures 3B–F**) can be assigned to certain process events, i.e., to changing process phases or LHP cleaning cycles. The PCA score plots for all batches can be found in **Figure S-3**.

**X**-variances captured by individual PCs (**Table S-3** in section S2.2.2 of **Supplementary Material**) indicate at least six significant factors for all preprocessing methods, while the PC8–PC10 are definitely negligible. The PC7 seems to be a boundary case, and its significance should be proved by using other criteria. Considering spectrum-like loadings (**Figure S-4**) and process-reflecting scores, in particular in the time-wise averaged data (**Figure 3C**), seven PCs are likely to be relevant. Additional considerations helping to deduce a number of PCs in the augmented process data are considered in section S2.2 of **Supplementary Material**.

In general, the low variances captured by minor principal components PC2-PC7 (**Table S-3**) illustrate a much higher sensitivity of NIR spectroscopy to water than to other chemical or physical variability sources in the drying process medium. Nevertheless, a thorough study of the complete PCA model resulted in some practically important observations. Thus, LHP fouling and cleaning during the process has a minor effect on the observed in-line spectra, in particular, at the final process stage (section S2.2.1 of **Supplementary Material**).

An exploratory data analysis performed has revealed an essential correlation of all spectral variables with the moisture content. The PCA analysis of the united dataset (16,303 spectra) has shown that this effect is overlaid with a variation on the stochastic spectrum intensity caused by the process noise. Since both scatter-driven effects have similar spectral signatures, the application of conventional normalization or derivative preprocessing methods of scatter correction incidentally removes useful information contained in the spectrum background. Instead, it was suggested to perform the smoothing of spectral variables along the time domain, e.g., using a moving window average.

# Building an Accurate PLS Regression Model of Moisture Content

For efficiently using the additional moisture-related information contained in the spectral variables, the dependence of model accuracy on averaging window width (WW = k points) has been studied. PLS models for all possible odd k values between 3 and 101 in different moisture ranges were compared (section S2.3.2 in **Supplementary Material**). Since the in-line smoothing of time dependencies results in a delay of 2k – 1 trajectory points (half WW) between the process and analysis times (Bogomolov, 2011), light smoothing is technically preferred. WW = 15 was found to be optimal in all cases as it provided an essential improvement of the model accuracy with a reasonable delay of 35 s. The full WW of 70 s approximately corresponds to a material circulation period in this process and dryer type. Thus, each portion of the granulate has a good chance of being exposed to spectroscopic measurement during this time. Due to the averaging, a virtual sample size captured by spectroscopic measurement, and hence the level of scrutiny of analysis, is extended. From this point of view, an optimal WW should correspond to an averaged spectrum that is representative of the bulk material volume, while remaining a nearly instant measurement compared to the total process time. This principle can be suggested as a rule of thumb for optimal data averaging in the drying process analysis and similar applications. A 47-point averaging was found to be a "global" optimum in our case; stronger smoothing does not lead to any significant gain. Based on these observations, 15- and 47-point smoothing windows have been chosen as benchmarks for model comparison (the respective preprocessing methods are designated as S15 and S47). **Table 1** presents a summary

of full-spectrum modeling results for different moisture ranges, preprocessing techniques, and validation methods.

(violet line); and LHP cleaning start/end points (vertical green lines).

The data covers a wide range of moisture contents from 2 to 26% (**Table S-1**). As the prediction error may be nonuniform depending on the drying stages (Mantanus et al., 2009), several PLS models were built corresponding to moisture LOD ranges <20% (D20), <15% (D15), and <10% (D10), in addition to the full-data (D) models. The abundance of measurement points makes possible the use of this data reduction without a significant impact to the model quality. The upper value of moisture content noticeably reduces the RMSE (e.g., for LBO CV, it falls from 0.21 in D to 0.13 in D10), keeping R 2 at the same high level of 0.997–0.998 (**Table 1**). A strong error dependence on the moisture content can be practically employed to improve the performance of moisture monitoring in general. Thus, prediction software can switch to a more precise model as soon as a certain moisture content level is reached, providing an automatic model "focusing" in the process course. By this way, the most critical final stage of drying can be monitored with the highest accuracy.

A number of LVs to be kept in PLS models was estimated from the RMSE of different validation methods and from the explained **X**- and **y**-variances (**Table S-4**). **Figure 4** compares the LBO CV RMSE dependencies on the number of LVs for the models in different moisture ranges (**Figure 4A**) and data averaging degrees (**Figure 4B**). Their common trend is that the validation error reaches a plateau starting from the seventh LV; faint minima at higher factor numbers do not seem significant. Note that LBO TABLE 1 | PLS regression statistics for in-line moisture content determination: model comparison for different moisture ranges and preprocessing techniques using different validation methods; all models were built with 7 LVs.

aDataset used: D – full dataset, D20, D15, and D10 – datasets limited to LOD moisture content below 20, 15, and 10%, respectively; b the number of samples without outliers (see section S2.3 of Supplementary Information); cpreprocessing applied; dcalibration statistics; e full cross-validation statistics; f leave-a-batch-out cross-validation statistics; g validation set (Table 1) prediction statistics; h variable averaging with 15-point window; i variable averaging with 47-point window; jSavitzky–Golay first derivative with second-order polynomial and 15-point smoothing window.

CV is generally the most conservative (i.e., resulting in the highest errors) validation method in **Table 1**. Data scatter correction does not result in any model simplification as expected. **Figure 4B** shows that the validation RMSE for MSC-preprocessed D15 data is even higher than the RMSEV of the model obtained after moderate (S15) data smoothing. This effect is observed for any number of LVs higher than one. Starting from the sixth LV, the prediction error after MSC becomes even worse than in the raw-data model. This behavior agrees with the earlier PCA-based conclusion that conventional scatter correction refines only the two first factors of the multivariate space, transferring the process noise into higher yet significant model dimensions.

The analysis of the captured **X**- and **y**-variances (**Table S-4**) exhibited similar trends. It was also shown that seemingly insignificant variances captured by the seventh LV in the calibration data are still in agreement with the respective precisions of the NIR spectrometer and the LOD analyzer (section S2.3.3 of **Supplementary Material**).

The first two PLS loadings (**Figure S-6**) are almost identical to those in the augmented PCA (**Figure S-4**); therefore, both multivariate modeling spaces are essentially the same. Meaningful shapes of the first seven loadings, which are similar in PCA (**Figure S-4**) and PLS models (**Figure S-6**) as well as PCA scores (**Figure 3**), provide an additional justification of the chosen model's complexity. The noticeable positive offset of **p**1 in raw and smoothed data models (**Figures S-6a-c**) indicates that PLS regression makes use of both absorbance and scattercorrelated variances for moisture calibration. The loadings **p3** to **p7** still exhibit similar (as in PCA) interpretable spectrum-like features. Therefore, seven LVs were found to be optimal for all moisture ranges and data preprocessing methods, consistent with the earlier PCA result for all spectral data. This number is also reasonable, considering the physical and chemical complexity of the process as well as the anticipated nonlinearity of spectral responses. It is also acceptable from the point of view of calibration set size.

**Table S-4** also confirms the efficiency of variable smoothing. Cumulative **y**-variances grow with the averaging WW, reducing a misbalance between the **X**- and **y**-variances for any number of LVs, in particular, for LV1. Starting from LV3, the **y**-variance captured in smoothed data becomes higher than that in the models preceded by scatter correction (e.g., MSC). In detail, the problem of deducing the optimal number of LVs is considered in section S2.3.3 (**Supplementary Material**).

Validation statistics presented in **Table 1** evidences that the suggested data averaging approach is advantageous as compared to the MSC, SNV, and first derivative using the Savitzky–Golay smoothing filter. It is remarkable that any scatter correction (most essentially, MSC or SNV) leads to higher calibration and validation errors than those for raw spectral data (This comparison is provided for D and D15, but it holds for all datasets). **Figure 5** illustrates the model performance achieved in D15.

FIGURE 4 | RMSE dependencies (LBO CV) on the number of LVs in PLS models: (A) for nonpreprocessed data in different moisture content ranges: D (squares), D20 (diamonds), D15 (circles), and D10 (triangles); and (B) D15 data with different smoothing degrees: none (solid), S15 (dashed), and S47 (dash-dotted), as well as for MSC preprocessing (red dotted, filled markers).

A subset of 101 process samples was additionally analyzed off-line by using a high-resolution FT-NIR spectrometer (section S2.3.5 of **Supplementary Material**). The integration sphere applied in this case excluded any scatter-related stochastic variation of spectral intensities. Nevertheless, all spectral variables (including the background signal) exhibited the same strong correlation with the sample moisture content (**Figures S-9, S-10**), as in the case of in-line spectra (**Figure 2**). This fact confirms our previously given explanation of this effect in terms of changing light propagation conditions. Moreover, the performance of the PLS model built on 96 off-line spectra (samples with LOD > 15% were used) was found to be essentially the same (cross-validation RMSE = 0.108) as in the model built on respective averaged in-line spectra (S15) of the same process samples (**Table S-6**). This remarkable result provides an additional confirmation of the efficiency of the suggested method. For more details on the off-line analysis results, see **Supplementary Material**, section S2.3.5.

The time dependencies of the predicted moisture content in B12 (**Figure S-7**) illustrate the additional advantages of the suggested preprocessing technique. Variable smoothing most efficiently eliminates the noise contained in process trajectories at the beginning of the drying process, when the moisture content is greater than 15%. It also helps avoid prediction artifacts related to probe cleaning during the "wet" process stage. Section S2.3.4 in **Supplementary Material** provides a detailed discussion of the predicted drying trajectories.

In numerous publications on in-line diffuse-reflectance NIR monitoring of fluidized bed drying and similar processes, data analysis is always prefaced by MSC, SNV, or derivatives without exception. A mandatory application of corrective preprocessing may only be justified in preliminary feasibility studies, when the small calibration/validation dataset does not allow for building models of adequate complexity. The results reported here could be used as evidence for the destructiveness of scatter correction for the moisture calibration, as it eliminates a significant portion of the useful variance. Similar ideas have been formulated in the literature (Chen and Thennadil, 2012), where the information content of MSC coefficients was analyzed. The PLS capability of employing the quantitative information delivered by the scatter has earlier been illustrated in other applications, in particular in particle size analysis (Nieuwmeyer et al., 2007) and the quantitative determination of fat and protein in milk (Bogomolov et al., 2012). In these cases, the predictive models built on raw data exhibited a noticeably better performance, as compared to those in which any scatter-correction was applied. For in-line process data, the suggested smoothing approach, performed in a time rather than spectral domain, presents a viable alternative to the classic scatter correction of spectra, to eliminate noise while preserving useful information contained in the spectral variables.

# CONCLUSIONS

In light of our presented results, the following recommendations to practical NIR spectroscopic monitoring of moisture content in fluidized bed drying and similar process types can be formulated. A very common practice of a priori scatter correction of inline process spectra prior to the multivariate calibration is generally discouraged, because it may eliminate an essential part of the water-related variance from the data and thus deteriorate the resulting prediction model. To avoid this, quantitative modeling should be prefaced by an exploratory analysis of the raw data to investigate the relevance of both absorbance and scatter-related effects of moisture by using a sufficiently large representative set of designed samples and process conditions. These considerations are equally valid in cases when water content is not directly determined, but it should be taken into account by an accurate multivariate model as an important process factor. Process noise, i.e., stochastic background and intensity variations of in-line spectra, can be efficiently eliminated with a minimal loss of useful information by means of data smoothing along the time scale. The parameters of smoothing strengths should be adjusted depending on the process scale and dynamics. Building accurate quantitative models should rely on a methodically determined number

REFERENCES

of latent variables. A deliberate application of less LVs than their optimal number following from the model diagnostics sometimes done by researchers to guarantee an avoidance of overfitting—is not always justified. An underfitting may often be more undesirable for model prediction accuracy.

# AUTHOR CONTRIBUTIONS

AB conceived and wrote the paper and analyzed the data. JM conceived the project and organized and planned industrial experiments. OH performed the experiments and analyzed the data.

# ACKNOWLEDGMENTS

The Ministry of Education and Science of the Russian Federation supported this work within the framework of the basic part of state task on the theme Adaptive technologies of analytical control based on optical sensors (Project No. 4.7001.2017/BP). The authors thank Tomas Vermeire (GEA, Belgium) and the colleagues from the previous Pharmaceutical and Analytical Development department (Weesp, The Netherlands) for their support of the experiments. J&M Analytik AG is acknowledged for organization and support. Prof. Dr. Rudolf W. Kessler (Reutlingen University, Germany) is acknowledged for fruitful discussions. Ivan and Petr Bogomolov helped in manuscript preparation.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2018.00388/full#supplementary-material

particle size from these parameters. Anal. Chim. Acta 746, 37– 46. doi: 10.1016/j.aca.2012.08.006

spectroscopy with consideration of sampling effects on method accuracy. Anal. Chem. 77, 4515–4522. doi: 10.1021/ac050272q

**Conflict of Interest Statement:** OH is an AbbVie employee and may own AbbVie stock/options.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Bogomolov, Mannhardt and Heinzerling. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Eliminating Non-linear Raman Shift Displacement Between Spectrometers via Moving Window Fast Fourier Transform Cross-Correlation

### Hui Chen1,2,3†, Yan Liu1†, Feng Lu1 \*, Yongbing Cao2,4 \* and Zhi-Min Zhang5 \*

### Edited by:

Hoang Vu Dang, Hanoi University of Pharmacy, Vietnam

### Reviewed by:

Andreas Borgschulte, Swiss Federal Laboratories for Materials Science and Technology, Switzerland Sebastian Primpke, Alfred Wegener Institut Helmholtz Zentrum für Polar und Meeresforschung, Germany

### \*Correspondence:

Feng Lu<EMAIL_ADDRESS>Yongbing Cao<EMAIL_ADDRESS>Zhi-Min Zhang<EMAIL_ADDRESS>
†These authors have contributed equally to this work and are co-first authors

### Specialty section:

This article was submitted to Analytical Chemistry, a section of the journal Frontiers in Chemistry

Received: 08 February 2018 Accepted: 05 October 2018 Published: 25 October 2018

### Citation:

Chen H, Liu Y, Lu F, Cao Y and Zhang Z-M (2018) Eliminating Non-linear Raman Shift Displacement Between Spectrometers via Moving Window Fast Fourier Transform Cross-Correlation. Front. Chem. 6:515. doi: 10.3389/fchem.2018.00515 1 School of Pharmacy, Second Military Medical University, Shanghai, China, 2 Department of Vascular Disease, Shanghai TCM-Integrated Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai, China, 3 Quality Control Department, Shanghai Diracarta Biomedical Technology Co., Ltd, Shanghai, China, 4 Department of Foundation and New Drug Research, Shanghai TCM-Integrated Institute of Vascular Disease, Shanghai, China, 5 College of Chemistry and Chemical Engineering, Central South University, Changsha, China

Obtaining consistent spectra by using different spectrometers is of critical importance to the fields that rely heavily on Raman spectroscopy. The quality of both qualitative and quantitative analysis depends on the stability of specific Raman peak shifts across instruments. Non-linear drifts in the Raman shifts can, however, introduce additional complexity in model building, potentially even rendering a model impractical. Fortunately, various types of shift correction methods can be applied in data preprocessing in order to address this problem. In this work, a moving window fast Fourier transform cross-correlation is developed to correct non-linear shifts for synchronization of spectra obtained from different Raman instruments. The performance of this method is demonstrated by using a series of Raman spectra of pharmaceuticals as well as comparing with data obtained by using an existing standard Raman shift scattering procedure. The results show that after the removal of shift displacements, the spectral consistency improves significantly, i.e., the spectral correlation coefficient of the two Raman instruments increased from 0.87 to 0.95. The developed standardization method has, to a certain extent, reduced instrumental systematic errors caused by measurement, while enhancing spectral compatibility and consistency through a simple and flexible moving window procedure.

Keywords: Raman instruments, shift correction, cross-correlation, fast fourier transform, moving window

# INTRODUCTION

Over the last few decades, the use of Raman spectroscopy in combination with chemometric methods has increased significantly for analysis of pharmaceutical products (Sacré et al., 2010; Dégardin et al., 2011; Loethen et al., 2015), detection of food adulteration (Zou et al., 2009; Cheng et al., 2010), and other applications (Mrozek et al., 2004; Taleb et al., 2006; Muehlethaler et al., 2011). Raman spectroscopy is a powerful tool for sample analysis and benefits from several advantages such as high speed, simplicity, non-destructive nature, and cost-effectiveness. To date,

**139**

it has been extensively applied in pharmaceutical analysis by constructing multivariate calibration models. However, these models will be invalid if an existing calibration model is applied to spectra that are collected on a different occasion or a separate instrument, or when the response of an old instrument suffers from variations (Du et al., 2011; Brown, 2013). These variations may, if left untreated, dominate the calibration models, thereby making analysis of samples impractical. Consequently, chemometric techniques have been used to circumvent these problems through instrumental transfer or standardization so as to isolate and compensate for any instrumental and environmental variations.

A number of methods, including both instrumental transfer and standardization, have been discussed in the literature (Wang et al., 1991, 1992; De Noord, 1994; Mann and Vickers, 1999; Nguyen Quang et al., 1999; Hutsebaut et al., 2005; Kompany-Zareh and van den Berg, 2010; Rodriguez et al., 2011b; Weatherall et al., 2013). The direct standardization (DS) and piecewise direct standardization (PDS) developed by Wang et al. (1991, 1992) are the most extensively used procedures for spectral response standardization. Using the PDS method, Gryniewicz-Ruzicka (Gryniewicz-Ruzicka et al., 2011) obtained a very low detection limit for diethylene glycol in pharmaceuticalgrade glycerin by using five portable Raman spectrometers. This method, however, requires the user to measure several standards prior to analyzing samples. In addition, both the use of the moving window strategy and the selection of principal components have a noticeable impact on the performance of PDS, which needs to be determined carefully. Furthermore, neither the DS nor the PDS method can deal with different (i.e., non-linear) shifts in the peaks in Raman spectra. It is worth mentioning that in contrast to the various instrumental spectral responses, Raman shift inconsistencies arise mainly from different charge-coupled device (CCD) detectors (Vickers and Mann, 1999). Nonetheless, the use of inconsistent spectra will diminish significantly the predictive power of a calibration model. As a result, the removal of Raman shifts or wavelength inconsistencies for spectra synchronization has become a particularly significant aspect of Raman spectroscopy analysis. In 1996, a mathematical procedure to correct wavelength drifts to synchronize Raman spectra was presented by Booksh et al. (1996). Typically, empirical data are required to select a number of principal components and channels to increase the synchronization precision. Westad and Martens (1999) developed a more general concept of shift determination and tested it on Raman spectra. The results revealed, however, that the spectra were not reproduced exactly after removal of peak drifts exceeding a discrete spectral resolution. Hutsebaut et al. (2005) used a Raman shift standard scattering (SSS for short) method in combination with a linear fitting to determine shift drifts between measured Raman peak and reference positions. A similar approach was used by Rodriguez et al. (2011b) to transfer Raman spectral libraries among instruments. Nevertheless, the use of Raman shift standards is inappropriate for in-line monitoring applications as a result of the difficulties associated with incorporating one or more of the materials proposed as shift standards in a system for in-line measurements. Recently, another approach for the removal of disturbing factors in the CCD responses and instrumental apparatus functions was proposed by Weatherall et al. (2013). Unfortunately, the use of baselineWavelet continuous wavelet transform as a function to identify major peaks' positions accurately requires idealized line profiles of the corresponding peaks, which is not practical for real Raman spectra. In addition, several parameters that influence the final results, such as the width of the window and the choice of the signal-to-noise threshold, need to be specified, mostly by the users.

As a result of the multifarious theoretical and practical limitations of the existing instrument standardization methods (Chen et al., 2015), there is a significant demand for methods that are easier to implement (i.e., fewer or even no tunable parameters required) in order to acquire better analytical performance. Accordingly, we introduced a cross-correlation method in order to address the problems (such as tunable parameters, need idealized line profiles, etc.) discussed above. Generally, in signal processing, cross-correlation is a measure of similarity of two waveforms as a function of a time-lag applied to one of them, and is also known as a sliding dot product or sliding inner-product (Welch, 1974; Goshtasby et al., 1984). When coupled with fast Fourier transform (FFT) algorithms, the efficiency of FFT can be exploited in the numerical computation of cross-correlations, accelerating thus the convolution calculation (Bracewell, 1980). FFT cross-correlation may therefore be the fastest method in signal processing for shift correction (Bergland, 1969), and benefits from many advantages such as high speed and accuracy. Moreover, it also eliminates the requirement for alignment parameters. Previously, two alignment methods were proposed to estimate the shifts between segments in large chromatographic and spectral datasets, namely, peak alignment by FFT (Wong et al., 2005b) and recursive alignment by FFT (Wong et al., 2005a). However, these two methods move segments by insertion and deletion of data points at the start and end of segments, without considering peak information, which may cause changes in the shapes of peaks by introducing artifacts and removing peak points. Zhan et al. (Zhang et al., 2012) developed another method, known as the multi-scale peak alignment (MSPA) method, to synchronize peaks against a reference chromatogram (aligning peaks from large to small scales), which is accelerated by the application of FFT cross-correlation while preserving peak shape during synchronization. Similarly, Li et al. (2013a) developed a moving window FFT cross-correlation (MWFFT) method to effectively synchronize high-throughput chromatograms without segment size optimization. However, the Raman spectra profiles were different from the chromatograms, which required peak fitting to obtain perfect profiles and a precise Raman shift.

In the present work, the MWFFT was improved and subsequently applied to spectral standardization to address the issues associated with spectral drifts in Raman spectrometers. The performance of this method was compared to that of the SSS method (Hutsebaut et al., 2005) by using two Raman datasets from primary and secondary spectrometers. The aim of our study was to make the MWFFT as a powerful and practical method for standardization across

Raman spectrometers, which can be easily implemented and well-suited for solving Raman shifts displacements between spectrometers.

# MATERIALS AND METHODS

# Standards and Samples

Standards (acetaminophen and cyclohexane) were provided by the National Institute for the Control of Pharmaceutical and Biological Products. Pharmaceutical tablets (listed in **Table 1**) from five different manufacturers were provided by the Shanghai Institute for Food and Drug Control.

# Raman Spectrometers

Two Raman instruments with an excitation wavelength of 785 nm were used, and their physical parameters are listed in **Table 2**. In this work, the i-Raman is regarded as the "master" (primary) instrument, while the GemRam is regarded as the "slave" (secondary) instrument.

The integration times of the standards and drugs were of 2 and 3 s, respectively. Unless stated otherwise, six Raman spectra were collected for each drug during the experiment. It is worth noting that the final spectrum of each drug was calculated as the average of spectra collected from a variety of positions. Moreover, only the spectral region containing the most abundant information (i.e., 300–1,700 cm−1 ) was used in subsequent data analysis.

# Cross-Correlation

In signal processing, cross-correlation is a standard technique to calculate the similarity between and estimate the linear shift of two signals as a function of one relative to the other, which is also known as the sliding dot product. It is obvious that any changes involving the shifting of one signal will affect the correlation coefficient calculated for any combination of two signals that includes this shifted signal. For two discrete signals such as those in the Raman spectra, the cross-correlation is defined as:

$$\mathbf{c}(i) = \frac{\sum\_{i} \left( \mathbf{r}(i) - \bar{\mathbf{r}} \right) \left( \mathbf{s}(i+j) - \bar{\mathbf{s}} \right)}{\sqrt{\sum\_{i} \left( \mathbf{r}(i) - \bar{\mathbf{r}} \right)^{2}} \sqrt{\sum\_{i} \left( \mathbf{s}(i+j) - \bar{\mathbf{s}} \right)^{2}}} \tag{1}$$

where r is the reference signal, s is the signal to be synchronized, c is the cross-correlation values for all lags. As a simple example, consider two simulated Raman spectra r and s that differ only by a known displacement of 90 points along the x-axis. We can determine by how much s be shifted along the x-axis in order to maximize its similarity to r by using cross-correlation. The above formula slides s along the x-axis, calculating the sum of their product at each position. When the value of c is maximized, i.e., the signals match well due to peak synchronization, they make the most significant contribution to the sum of their product. A visual description of the calculation procedure of cross-correlation and estimation of shifts between signals via cross-correlation is shown in **Figure 1**.

TABLE 1 | Correlation coefficients of drug tablets before and after shift correction.

uCorrelation coefficient before shift correction; sCorrelation coefficient by SSS; mCorrelation coefficient by MWFFT.

# Moving Window FFT Cross-Correlation

The FFT is typically used to calculate the cross-correlation between 1D and 2D signals (Papoulis, 1962; Cooley et al., 1969; Dutt and Rokhlin, 1993). In the present work, FFT was used to increase the speed of cross-correlation between two datasets, in which one signal may be shifted relative

to another. In addition, and perhaps more significantly for its application to the spectral synchronization problem, FFT cross-correlation is not heuristic and thus can identify consistently the best match between signals by finding the maximum correlation coefficient (Wong et al., 2005b).

Usually, the cross-correlation method can only estimate linear shifts between Raman spectra. However, Raman shift displacements are often non-linear in real samples. Consequently, we adopted the moving window procedure in this work to address this problem. In this procedure, the shifts relative to the reference can be estimated by FFT crosscorrelation, allowing us to obtain the shift profiles of all samples. Furthermore, MWFFT can be implemented and optimized simply and effectively only if a moving window of appropriate size is utilized. With a window moving from the beginning to the end of the two spectra, one can obtain a matrix of shift points. Accordingly, the shift profile can be obtained by calculating the mode value of each column of the shift matrix. **Figure 2A** shows an example Raman shift profile estimated by using the moving window strategy and FFT cross-correlation. It is apparent from the obtained shift profile that non-linear shifts exist across the entire spectral region, while the change points are observed in two regions with different shifts. By moving the continuous region around the change points, the synchronization procedure can be finished smoothly to obtain the synchronized spectrum, which can be seen in **Figure 2B**, with all the non-linear shifts successfully synchronized.

# RESULTS

There are two common ways of correcting the x-axis in Raman spectrometers (McCreery, 2005). The first one is to simply use the SSS method (Hutsebaut et al., 2005); the second one is based on absolute frequency calibration using the emission line spectra of gases. The SSS method, which requires the acquisition of Raman spectra of common materials with well-established Raman shift peak frequencies in order to correct the Raman shift axis directly, is used as a comparative method in this work. Several well-known Raman shift chemical standards, namely, cyclohexane and acetaminophen, are chosen over others for this study since their spectral combination can provide more signals in the region from 300 to 1,700 cm−1 (see **Table 3**). The left panel in **Figure 3** shows the spectra acquired for the used chemical standards on two instruments, while the right panel shows a plot of their differences. It should be noted that when the SSS method was used, the spectra acquired on the primary instrument were regarded as the reference, i.e., the peak positions in these spectra were used for synchronization. The relevant peak positions obtained on the secondary instrument are compared to those obtained on the primary instrument and are subsequently subtracted from the primary peak positions to afford the corresponding shift displacements. Linear fitting is then used to describe the shift displacements between the two instruments. Finally, the shift correction is carried out by linear interpolation.

# Synchronization of Pharmaceutical Datasets

Data synchronization of the raw Raman spectra are presented to evaluate the performance of the MWFFT method (**Figure 2**). In order to gain further insight into the two shift correction algorithms, and the properties and advantages of MWFFT in particular, different batches of pharmaceutical tablets were examined to verify the practicability and effectiveness of MWFFT. **Figure 4** describes the application of MWFFT—each tablet from a total of 40 drugs was analyzed on average six times on two instruments to obtain six different spectra. Subsequently, these spectra were detected for outlier. The average spectrum obtained from six spectra acquired on the primary instrument can be regarded as a reference without outliers. Analogously, we obtained the spectrum of the same tablet on the secondary instrument, and this represents the spectrum to be synchronized. Finally, MWFFT was applied to remove the shift displacements in order to synchronize the spectra across the two instruments.

Prior to data analysis, adaptive iteratively reweighted penalized least squares (airPLS) (Zhang et al., 2010a,b; Li et al.,

2013b) baseline correction and Savitzky–Golay smoothing (Savitzky and Golay, 1964) (a 9-point wide window and a second-order polynomial) were used in the preprocessing of a variety of pharmaceutical datasets. All processing tasks were implemented on a personal computer (CPU: 2.53G, RAM: 8GB) with MATLAB R2013a. Firstly, we demonstrate the effect of MWFFT by using the pharmaceutical datasets (**Figure 5**). The primary instrument spectra (black lines) are used as references for synchronization. **Figure 5** shows the magnified versions of the sample profiles, focusing on a particular set of peaks in order to allow the performance of the MWFFT method to be evaluated by visual inspection. For the acyclovir and captopril datasets, it can be seen that before synchronization (top panel in **Figure 5**), the peaks in the spectrum collected on the secondary instrument are de-synchronized with respect to that obtained on the primary instrument, and vary from sample to sample. After synchronization (middle panel in **Figure 5**) using the MWFFT method, it is apparent that all the spectra are now properly synchronized. This outcome is attributed to the action of the MWFFT method, which appropriately slides the peaks to match the reference spectrum with a window size of 70 points. In addition, for the sake of comparison, the results obtained using the SSS method for the same spectra are displayed in the bottom panel of **Figure 5**.

# Correlation Coefficient After Synchronization

The correlation, or distance, between a signal and the reference point is often used as an optimization objective function—when the signals match, the correlation coefficient is maximized. In this case, correlation coefficient can be used to assess the synchronization problem (Lee Rodgers and Nicewander, 1988). Generally, the correlation coefficient is a good descriptor of similarity, with a value of 1.00 indicating a perfect match, while 0 indicates significant dissimilarity. The correlation coefficient is simple to use and possesses several desirable properties, which we discussed in detail in our previous work (Gao et al., 2014). The correlation coefficient can be calculated by using the following equations:

$$\begin{aligned} r &= \frac{\sum\_{i=1}^{n} \left( X\_i^p - \bar{\mathbf{x}}^p \right) \left( X\_i^s - \bar{\mathbf{x}}^s \right)}{\int \sqrt{\sum\_{i=1}^{n} \left( X\_i^p - \bar{\mathbf{x}}^p \right)^2 \sum\_{i=1}^{n} \left( X\_i^s - \bar{\mathbf{x}}^s \right)^2} \\\\ \sum\_{i=1}^{n} \left( X\_i^p - \bar{\mathbf{x}}^p \right) \left( X\_i^{sa} - \bar{\mathbf{x}}^{sa} \right) \end{aligned} \tag{2}$$

$$R = \frac{\sum\_{l=1}^{n} \left(\mathbf{x}\_{l}^{n} - \hat{\mathbf{x}}\_{l}^{n}\mathbf{x}\_{l}^{n} - \hat{\mathbf{x}}\_{l}^{n}\right)^{2}}{\int \sqrt{\sum\_{l=1}^{n} \left(\mathbf{X}\_{l}^{p} - \bar{\mathbf{x}}^{p}\right)^{2} \sum\_{l=1}^{n} \left(\mathbf{X}\_{l}^{sa} - \bar{\mathbf{x}}^{sa}\right)^{2}}} \tag{3}$$

Here, Xp and Xs represent the spectra of n drugs measured on the primary and secondary instruments, respectively. Parameters x p , x s , and x sa represent the average spectra of Xp , Xs , and X sa, respectively. Xsa indicates the secondary shift corrected spectrum, while r and R denote a similarity between the primary original spectrum and the secondary spectrum (before or after shift correction). The correlation coefficient of each drug's spectrum was calculated, and the results are summarized in **Table 1**. During the preprocessing, linear interpolation was used to re-compute intensity based on the master Raman shift xaxis in order to unify the spectra obtained using the primary and secondary instruments. It is apparent from **Table 1** that the correlation coefficients between the two instruments improved significantly after shift correction.

As can be seen in **Table 1**, the correlation coefficient assessment prior to the shift correction exhibited a slight variation among different batches of a drug. Nevertheless, these variations are within the three-sigma range. Portable spectrometers are often based on the use of library-based spectral correlation methods (Carron and Cox, 2010), which frequently utilize the hit-quality index (HQI) as the figure of merit to characterize the correlation with each other. The

TABLE 3 | Raman shifts (cm−1 ) used to calibrate standard samples.

aValues as reported by ASTM E1840-96.

typical minimum threshold that classifies an unknown sample as a "Pass" is 0.95 (Rodriguez et al., 2011a, 2013), which is similar to the correlation coefficient. Clearly, the MWFFT method makes a significant contribution to the level of similarity for the spectra obtained using the slave instrument. The synchronization increased the similarities for all drugs above the verification threshold of 0.95, while the similarity for one captopril tablet remained under 0.95 when the SSS procedure was used. Consequently, it is obvious that the MWFFT method can correct the non-linear shifts successfully, synchronizing thus the secondary spectra to the reference spectra in a timeeffective manner. In addition, MWFFT can reduce the systematic differences across spectrometers, which can increase the spectral consistency of different instruments as well as the compatibility with library search. Furthermore, this method can be used as an on-line standardization method across Raman instruments in the future.

# DISCUSSION

# Selection of Reference Spectrum

A wide application of the MWFFT method necessitates the selection of an appropriate reference spectrum. When a drug sample is measured on a secondary instrument to obtain an average spectrum for synchronization, its corresponding standard spectrum contained in the existing spectral library can

FIGURE 3 | Spectra of acetaminophen (A) and cyclohexane (C) acquired on two different instruments. Magnified spectral differences in (B,D) correspond to the shaded areas in (A,C), respectively.

be certainly used as the reference to correct shift displacements. However, when the spectral library does not contain the required reference spectrum, it would be preferable to use the reference spectra of existing drugs with the same generic name in the database in order to obtain a new matrix of shift points. As a result, the shift profile of the new drug can be calculated from the mode of each column of the matrix. Through this profile, one can obtain a new reference spectrum by shift correction, which can be subsequently applied. Otherwise, one can regard the new sample spectrum directly as a reference, and save it in the database for subsequent analysis. The entire procedure is depicted in **Figure 6**.

# Avoiding Peak Detection Using the Moving Window Strategy

The existing peak detection methods, e.g., the wavelet and ridge line peak picking method, need idealized line profiles of the corresponding peaks in order to detect the displacements accurately, which is not practical for the spectroscopic analysis of real samples. Moreover, several parameters need to be specified

lines indicate the reference spectra. The inset shows the full Raman spectra, whereas the shaded areas indicate the region magnified in the main panel.

with a priori knowledge, which largely influence the final results and can be difficult to implement in C programming language. By contrast, the use of the moving window strategy can allow an estimation of non-linear shifts between spectra flexibly and without peak detection for peak synchronization. With a window that moves from the beginning to the end of two spectra, one can obtain an N-dimensional matrix of shift points, where the data points of a Raman spectrum are N. In this case, the shift profile can be calculated from the mode of each column of the matrix, while the mean and median of the matrix can outline the paths of the shifts. The Raman shift profile of metronidazole tablet is shown in **Figure 2A** using a green dotted line. It is apparent that the profiles in all regions are corrected by the MWFFT method, meaning this method is sufficiently flexible for estimation of non-linear shifts between spectra.

# Advantages of the MWFFT Method

The MWFFT method has several distinctive advantages when compared to the traditional methods as a result of the continuity and redundancy of the moving window procedure. Usually, the direct evaluation of cross-correlation requires O (N2 ) time complexity for a Raman spectrum of length N, which is time-consuming for spectra with thousands of data points. Fortunately, cross-correlation can be calculated by using FFT much more efficiently since it can significantly decrease the time complexity of crosscorrelation from O (N2 ) to O (NlogN). The use of the moving window strategy with FFT cross-correlation, with a window size w, leads to a time complexity of one window wlogw. Accordingly, the time complexity of MWFFT is Nwlogw, where N represents the number of data points in a Raman spectrum.

The MWFFT method evaluates the shift of each point. In the moving window strategy, only one parameter needs to be taken into account, which makes this procedure simple and practical, as there is no need for chemical standards. By contrast, the SSS method requires the use of some chemical standards in order to locate the position of each peak, which is used in turn to obtain the corresponding shift displacement. After the shift of each point is estimated by MWFFT, the points in the spectrum are shifted according to their shifts by insertion and deletion. The present work introduced a change point, i.e., a discontinuity point in the shift profile. It is possible to see that the change points (**Figure 2A**), around which insertions and deletions occur frequently, are not in the peak region. Consequently, peak distortions can be effectively avoided, allowing the peak shape to be preserved during the synchronization procedure with MWFFT. Overall, the advantages associated with the use of non-linear shift estimation, insertion and deletion around change points, and shape preservation make MWFFT a flexible, rapid, practical, and precise method for correcting shifts in synchronization of Raman datasets.

# Evaluation of the Synchronization Quality

Generally, Raman spectra will become more consistent, exhibit higher correlation coefficients, and be more similar to each other after a successful synchronization. The correlation coefficient can be used as a criterion for assessing the synchronization quality between the primary and secondary spectra. The synchronized spectra are commonly used to perform library-based searches and are further analyzed by chemometric algorithms. Usually, distance and Euclidean distance in particular (Juday, 1993),

acMean Euclidean distances of acyclovir datasets; caMean Euclidean distances of captopril datasets.

can also be a good criterion for evaluating the quality of synchronization. Generally, the more similar the spectra are, the smaller is the Euclidean distance between them, and vice versa. In this work, the mean Euclidean distance (Dmean) is calculated as follows:

$$D\_{mean} = \,^1\!/\_{\mathfrak{h}} \sum\_{i=1}^n \sqrt{\sum\_{j=1}^k \left(X\_{i,j}^p - X\_{i,j}^s\right)^2} \tag{4}$$

where the rows of matrix X correspond to observations (n), while the columns correspond to variables (k). Xp i and X s i are the ith primary (reference) spectrum and secondary spectrum, respectively. It is worth mentioning at this stage that the normalization algorithm (Heraud et al., 2006) is used to scale the spectra within a similar range before calculating the distances. The results are summarized in **Table 4**. It is apparent that the mean Euclidean distance of the pharmaceutical datasets shift-corrected by SSS and MWFFT were considerably reduced when compared to the uncorrected ones. In addition, for the two datasets, MWFFT performed slightly better than the SSS method in terms of non-linear shift correction.

# CONCLUSIONS

Methods for the synchronization of spectra are indispensable for successful applications using different spectrometers. In the present work, we used the moving window strategy in combination with FFT cross-correlation to synchronize Raman spectra. This technique, abbreviated as MWFFT, was shown to eliminate accurately and effectively non-linear shift displacements between Raman spectra. Owing to the continuity of the moving window technique, non-linear shifts are corrected and shift profiles are obtained for each spectrum. In general, the use of the FFT cross-correlation methodology is timesaving and results in a significant improvement in speed. Moreover, this method can reduce or even remove systematic differences between Raman spectrometers (a dramatic increase in similarity from 0.87 to 0.95 after synchronization of the spectra between master (primary) and slave (secondary) spectrometers), as well as the compatibility with Raman spectral library. It is better than the SSS method in terms of correcting non-linear shifts and does not require the use of Raman shift standards. These advantages make MWFFT a promising shift correction method that addresses the demand for automated, flexible, rapid, and reliable data preprocessing, which plays an important role in Raman spectroscopy analysis using different spectrometers. Finally, MWFFT can be easily implemented

REFERENCES

with C and C++ programming languages (available as open source package at http://code.google.com/p/mwfft), which may be well-suited to solving the Raman shift displacements between spectrometers in the fields that rely heavily on the use of Raman spectrometers.

# ETHICS STATEMENT

The experimental protocol was approved by the Research Ethics Committee of The Second Medical University and Shanghai University of Traditional Chinese Medicine, The findings and conclusions in this article have not been formally disseminated by the State Food and Drug Administration and should not be construed to represent any agency determination or policy.

# AUTHOR CONTRIBUTIONS

HC designed and carried out experiments. Z-MZ and YL assisted with analyzing the results and discussions. HC and Z-MZ wrote the manuscript. FL reviewed and edited the manuscript. YC reviewed and checked our manuscript, gave constructive amendments to the text, and also approved the version to be published.

# FUNDING

This work is financially supported by Ministry of Science and Technology of the People's Republic of China (2017YFF0210103, 2012YQ180132). The studies meet with the approval of the university's review board.

**Conflict of Interest Statement:** HC is employed by the company Shanghai Diracarta Biomedical Technology Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Chen, Liu, Lu, Cao and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Wavelet Transform-Based UV Spectroscopy for Pharmaceutical Analysis

### Erdal Dinç1 \* and Zehra Yazan2

1 Department of Analytical Chemistry, Faculty of Pharmacy, Ankara University, Ankara, Turkey, 2 Department of Chemistry, Ankara University Faculty of Science, Ankara, Turkey

In research and development laboratories, chemical or pharmaceutical analysis has been carried out by evaluating sample signals obtained from instruments. However, the qualitative and quantitative determination based on raw signals may not be always possible due to sample complexity. In such cases, there is a need for powerful signal processing methodologies that can effectively process raw signals to get correct results. Wavelet transform is one of the most indispensable and popular signal processing methods currently used for noise removal, background correction, differentiation, data smoothing and filtering, data compression and separation of overlapping signals etc. This review article describes the theoretical aspects of wavelet transform (i.e., discrete, continuous and fractional) and its characteristic applications in UV spectroscopic analysis of pharmaceuticals.

### Edited by:

Hoang Vu Dang, Hanoi University of Pharmacy, Vietnam

### Reviewed by:

Gaetano Ragno, Università della Calabria, Italy Joseph Dubrovkin, Western Galilee College, Israel

> \*Correspondence: Erdal Dinç<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Analytical Chemistry, a section of the journal Frontiers in Chemistry

Received: 27 April 2018 Accepted: 03 October 2018 Published: 26 October 2018

### Citation:

Dinç E and Yazan Z (2018) Wavelet Transform-Based UV Spectroscopy for Pharmaceutical Analysis. Front. Chem. 6:503. doi: 10.3389/fchem.2018.00503 Keywords: discrete wavelet transform, continuous wavelet transform, fractional wavelet transform, UV spectroscopy, pharmaceutical analysis

# INTRODUCTION

In experimental studies, instruments or devices can provide signals (or graphs) in different formats e.g., spectrum, chromatogram, voltammogram, and electroferogram etc. The analysis of chemicals and pharmaceuticals in various samples is based upon the utilization of the measured signals of substances of interest. In practice, such an analysis for a multicomponent mixture may not be determined without a prior separation step due to spectral overlapping. Therefore, high performance liquid chromatography (HPLC) is one of the most commonly used techniques for quantitative estimation in the quality control of raw materials and commercial products in laboratories. In some cases, chromatographic determination could not be possible due to not only similar physicochemical behavior of analytes but also time- and solvent-consumption for optimal experimental conditions.

In practice, UV spectroscopic methods are widely used in chemical and pharmaceutical analysis. As compared to chromatographic ones, the use of spectroscopic methods provides a rapid analysis with low-cost and acceptable results. However, multicomponent analysis may not be possible with a traditional UV spectrophotometric approach due to spectral interferences of both active and inactive ingredients in samples. In some cases, derivative spectrophotometry (O'Haver and Green, 1976; O'Haver, 1979; Levillain and Fompeydie, 1986; Ragno et al., 2006) and its improved versions e.g., ratio spectra-derivative spectrophotometry (Salinas et al., 1990), ratio spectra-derivative spectrophotometry-zero crossing (Berzas Nevado et al., 1992; Dinç and Onur, 1998; Dinç, 1999), and double-divisior-ratio spectra-derivative spectrophotometry (Dinç and Onur, 1998; Dinç, 1999; Gohel et al., 2014; Shokry et al., 2014) could be used in place of conventional UV spectrophotometric method for analysis of binary and ternary mixtures without using a separation step. However, these spectral approaches may not always yield successful data due to severely overlapping spectral bands, spectral noise and baseline variation. Additionally, high-order differentiation of spectra may lead to spectral deterioration i.e., a decrease in signal intensity and signal-to-noise ratio. As a result, a number of mathematical manipulations (or signal processing methods) are often required to make instrumental signals more meaningful for analysis purpose.

Generally speaking, transform (i.e., Fourier, Hilbert, shorttime Fourier, Wigner distribution, Radon, and wavelet) is a very suitable technique in the pre-treatment step to simplify signals. Fourier transform (FT) is the first method to modify chemical signal (Griffiths, 1977; Cooper, 1978; Griffiths and De Haseth, 1986; Ernst, 1989) with the mathematical essence such as filtering, convolution/deconvolution etc. FT analysis can localize signal in frequency domain very well, but not so much in time domain. In contrast, wavelet transform (WT) has the advantage of localizing signals both in time (position) and frequency (scale) domains, making it a preferable mathematical tool to replace FT in the study of the local property of a signal and the removal of the perturbation of measuring error in spectral analysis. Nowadays, WT is one of the most signal analysis algorithms commonly used in the different fields of chemistry and engineering, providing alternative ways or opportunities to resolve complex spectral bands or diverse data types of signals.

For readers interested in learning the general theory of wavelets, more details can be found in the literature (Mallat, 1988; Chui, 1992; Daubechies, 1992; Newland, 1993; Byrnes et al., 1994; Chui et al., 1994; Vetterli and Kovacevi ˇ c, 1995; Strang and ´ Nguyen, 1996).

In the signal smoothing and de-noising of spectral peaks, the elimination of noise requires an application of appropriate filters to the raw spectral data such as some conventional signal filters Savitzky–Golay, Fourier and Kalman (Brown et al., 1994, 1996). The use of WT in signal analysis is two-fold: (i) to detect the singularities of a signal very likely caused by high-frequency noise and (ii) to separate the signal frequencies at different scales (Palavajjhala et al., 1994; Yan-Fang, 2013; Li and Chen, 2014). To illustrate this, Barclay et al. (1997) performed a comparative study in de-noising and smoothing of Gaussian peak by using wavelet, Fourier and Savitzky–Golay filters i.e., smoothing eliminates high-frequency components of the transformed signal irrespective of their amplitudes, while de-noising eliminates small-amplitude components of the transformed signal irrespective of their frequencies.

Historically, WT principal applications in chemistry were first explored by Walczak and Massart (1997a), who presented an approach based on the application of wavelet packet transform (WPT) to the best-basis selection for the compression and denoising of a set of signals in time-frequency domain. In their paper, the proposed technique was compared to Wickerhauser's approach (Wickerhauser, 1994) of fast approximate principal component analysis (PCA). These authors also published two more papers on the application of wavelets for data processing i.e., the introduction of WPT for noise suppression and signal compression (Walczak and Massart, 1997b) and the use of WT for signal compression and denoising, image processing, data compression and multivariate data modeling in analytical chemistry (Walczak and Massart, 1997c). On the other hand, Alsberg et al. (1997) tried to introduce WT to chemometricians by suggesting the short-time FT technique as a resolution to obtain information about frequency changes over time as well as the WT for de-noising, baseline removal, determination of derivative zero crossings and signal compression. In 1997, WT application in chemical analysis was also confirmed by Wang et al. (1997) and Depczynski et al. (1997). Up to date, WT processing of the different types of raw signals has been reported for liquid chromatography (Shao et al., 1997, 1998a,b,c) and NMR spectroscopy (Neue, 1996; Barache et al., 1997), Raman spectra (Cai et al., 2001; Ehrentreich and Summchen, 2001), and voltammetry (Chen et al., 1996; Fang and Chen, 1997; Zheng et al., 1998; Zhong et al., 1998; Aballe et al., 1999; Zheng and Mo, 1999) IR and Raman spectroscopy (Shao and Zhuang, 2004; Hwang et al., 2005; Chalus et al., 2007; Jun-fang et al., 2007; Lai et al., 2011). In this context, as in the various fields of mathematics and engineering, the implementations of WT in analytical chemistry and neighbor disciplines has become increasingly attractive as an alternative way to analyze complex mixtures previously unresolved by traditional analytical techniques.

With reference to the above-mentioned review, the aim of this paper is to describe the fundamentals of WT methodologies and its typical implementations for UV spectroscopic analysis of pharmaceuticals.

# BRIEF HISTORY OF WAVELETS

In the literature, the first study was related to the Haar Wavelet transform. This family was suggested by the mathematician Alfred Haar in 1909. However, the word "wavelet" was not used in the period of Haar. In fact, the word "wavelet" was invented by Morlet and the physicist Alex Grossman in 1984. After the first orthogonal Haar wavelet, the second orthogonal wavelet known as "Meyer wavelet" was formulated by the mathematician Yves Meyer in 1985. In 1988, Stephane Mallat and Meyer elaborated the concept of multiresolution. In the same year, a systematical method to construct compactly supported continuous wavelets was found by Ingrid Daubechies. Afterwards, Mallat proposed the fast wavelet transform. The emergence of this algorithm increased the implementations of the WT in the signal processing field.

In other words, the history of the wavelet families could be given in the following chronological order: Haar families in 1910, Morlet wavelet concept in 1981, Morlet and Grossman, "wavelet" in 1984, Meyer, "orthogonal wavelet" in 1985, Mallat and Meyer, multiresolution analysis in 1988, Daubechies, compact support orthogonal wavelet in 1988 and Mallat, fast wavelet transform in 1989 (c.f. Chun-Lin, 2010).

Basically, WT can be mainly classified into discrete wavelet transform (DWT) and continuous wavelet transform (CWT) in the signal analysis. The theory and implementations of wavelets in chemistry and related fields were well documented as review papers (Leung et al., 1998; Dinç and Baleanu, 2007b; Dinç, 2013; Li and Chen, 2014; Medhat, 2015) and reference books (Walczak and Massart, 2000a,b; Walczak and Radomski, 2000; Brereton, 2003, 2008; Chau et al., 2004; Danzer, 2007; Mark and Workman, 2007; Dubrovkin, 2018).

## WAVELET TRANSFORM ALGORITHMS

FT is based upon the decomposition of a signal into a set of trigonometric (sine and cosine) functions i.e., FT represents a signal in terms of sinusoids. The representation of FT of a signal from time mode to frequency mode is illustrated in **Figure 1**. For the determination of a local information in the FT, it is required to use an analyzing function ψ having localization properties in both frequency and time domains. This ψ function is named as a wavelet and it must be wave of finite duration.

WT contains the decomposition of a signal into a set of basic functions (wavelets). Basis functions of WT are small waves detected in different times. On the contrary to FT, WT gives information on both time and frequency, making it as an alternative approach to eliminate the resolution problem in signal analysis.

By definition, wavelets are the mathematical methods that convert the data into various coefficients and then analyze each coefficient at a resolution corresponding to its scale. Projection of a signal onto wavelet basic functions is called the wavelet transform. In other words, wavelets are mathematical functions generated from a mother wavelet Ψ(x) by the scaling parameter (dilatation) and shifting parameter (translation) i.e., the signal is expanded on a set of the dilatation (scaling parameter) of functions

$$
\psi \left( \frac{x - a}{b} \right) \tag{1}
$$

The scaling parameter has a significant role for the variation of time and frequency resolution when processing the signal.

For a given mother wavelet (Daubechies, 1992) ψ (x) by the scaling parameter and shifting parameter o fψ (x), a set of functions expressed by ψa,b (x) is obtained from the following equation.

$$\psi\_{a,b}\ (\mathbf{x}) = \frac{1}{\sqrt{|a|}} \ \psi \ \left(\frac{\mathbf{x} - b}{a}\right), \mathbf{a} \neq \mathbf{0}, \mathbf{a}, \mathbf{b} \in \mathbb{R} \tag{2}$$

where a is the scaling parameter, b is the shifting parameter and R is domain of real number. The mathematical expression of a CWT on a function f (x) is given below

$$CWT\left\{f\left(\mathbf{x}\right);a,b\right\} = \int\_{-\infty}^{\infty} f(\mathbf{x}) \psi\_{a,b}^\*(\mathbf{x}) d\mathbf{x} = \left\langle f(\mathbf{x}), \psi\_{a,b} \right\rangle \tag{3}$$

here the superscript ∗ is related to the complex conjugate and hf(x), ψa,b i represents the inner product of function f(x) onto the wavelet function ψa,b(x).

The original signal can be completely reconstructed by a sampled version of the CWT. Usually, the exemplar is follows as

$$a = 2^{-m} \\ \text{and } b = n2^{-m} \tag{4}$$

Here a and b denote scale and dilatation parameters, respectively, and R is the real number. The expression of DWT can be given as

$$DWT = \int\_{-\infty}^{+\infty} f \text{ (X)} \ \psi\_{m\_\*n}^\*(\mathbf{x}) \, dt \tag{5}$$

Where ψ∗ m, n (X) = 2 −m ψ (2m x − n) is the dilated and translated version of the mother wavelet. In the application of the DWT, only outputs from the low-pass filter are processed by WT. However, in the wavelet packet decomposition of signals, both outputs from the low-pass and high-pass filters are manipulated by WT (Strang and Nguyen, 1996). Multiresolution decomposition with wavelets is an interesting topic for signal and image analysis (Mallat, 1988; Daubechies, 1992).

Some families of wavelets with names and their coding list are illustrated in **Table 1**.

For signal processing, there is also another WT approach i.e., fractional wavelet transform (FWT) specifically designed for rectification of the limitations of the WT and fractional FT (Blu and Unser, 2000, 2002; Unser and Blu, 2000). FWT is based on the fractional B-splines. As it is already known, the splines play an important role on the early development of the theory of WT.

TABLE 1 | Families of wavelets with names and their coding list.

A B-spline is generalization of the Beziers curve. Let a vector known as the knot be defined by T = {t0, t1, . . . , tm} where T is a non-decreasing sequence with ti ǫ [0, 1], and define control point P0, Pn. The knots t0, t1, . . . , tm is called internal knots. If p = m n - 1 denotes the degree, the basis function is defined as follows:

$$N\_{i,\ 0}(t) = f\left(\mathbf{x}\right) = \begin{cases} 1, \text{ if } t\_i \le t < t\_{i+1} \text{ and } t\_{i+1} \\ \quad 0 \text{ otherwise} \end{cases} \tag{6}$$

and

$$N\_{i,p}^{\\\uparrow}(t) = \frac{t - t\_i}{t\_{i+p} - t\_i} N\_{i,p-1} \ (t) + \frac{t\_{i+p+1} - t}{t\_{i+p} + 1 - t\_{i+1}} N\_{i+1,p-1} \ (t) \tag{7}$$

Therefore, the curve defined by

$$C(t) = \sum\_{i=0}^{n} P\_i \, N\_{i, \, p} \, (t) \tag{8}$$

is a B-spline

Fractional B-spline: The fractional B-spline is defined as

$$\beta\_+^a(\alpha) = \frac{\sum\_{k=0}^{+\infty} (-1)^k \binom{\alpha+1}{k} \binom{x-k}{k}\_+^\alpha}{\Gamma\left(\alpha+1\right)} \tag{9}$$

where Euler's Gamma function is obtained by

$$
\Gamma\_{\pm}(\alpha+1) = \int\_0^{+\alpha} \varkappa^{\alpha} \ e^{-\varkappa} \, d\varkappa \tag{10}
$$

and

$$(\boldsymbol{\chi} - k)^{+}\_{\alpha} = \max \left( \boldsymbol{\chi} - k, \ 0 \right)^{\alpha} \tag{11}$$

The forward fractional finite difference operator of order α is defined as

$$
\Delta\_+^{\alpha} f \begin{pmatrix} \mathbf{x} \end{pmatrix} = \sum\_{k=0}^{+\infty} (-1)^k \begin{pmatrix} \alpha \\ k \end{pmatrix} f \begin{pmatrix} \mathbf{x} - k \end{pmatrix}, \tag{12}
$$

where

$$\left(\frac{\alpha}{k}\right) = \frac{\Gamma\_{\mid}(\alpha+1)}{\Gamma\left(k+1\right)\left(\alpha-k+1\right)}\tag{13}$$

B-splines fulfill the convolution property, namely

$$
\beta\_+^{\alpha 1} \* \beta\_+^{\alpha 2} = \beta\_+^{\alpha 1 + \alpha 2} \tag{14}
$$

The centered fractional B-splines of degree α is defined as

$$\beta\_\*^{\alpha}(\mathbf{x}) = \frac{1}{\Gamma\_-(\alpha+1)} \sum\_{k \in \mathbb{Z}} (-1)^k \left| \frac{\alpha+1}{k} \right| \left| \mathbf{x} - k \right|\_\*^{\alpha} \tag{15}$$

where

$$|\mathbb{x}|\_{-}^{\alpha} = \begin{cases} \frac{|\cdot|^{\alpha}}{-2\sin\left(\frac{\pi}{2}\alpha\right)}, & \alpha \text{ not even} \\\frac{X\log x}{(-1)^{1+n}\pi}, & \alpha \text{ even} \end{cases} \tag{16}$$

The fractional B-spline wavelet is defined as

$$\begin{aligned} \psi\_+^{\alpha} \left( \frac{\chi}{2} \right) &= \sum\_{k \in \mathbb{Z}} \frac{(-1)^k}{2^{\alpha}} \\ &\sum\_{1 \in \mathbb{Z}} \binom{\alpha+1}{1} \beta\_\*^{2\alpha+1} \left( 1+k-1 \right) \beta\_+^{\alpha} \left( \chi - k \right) \text{(17)} \end{aligned}$$

We mention that the fractional splines wavelets of degree obey the following

$$\int\_{-\infty}^{+\infty} X^n \, \psi\_+^{\alpha} \, \left( \alpha \right) d\mathfrak{x} = 0, \dots, \left[ \alpha \right] \tag{18}$$

and the Fourier transform fulfills the following relations

$$
\hat{\psi}\_+^{\alpha}(\varpi) = \mathcal{C} \left( j\varpi \right)^{\alpha+1}, \text{ as } \varpi \to 0 \tag{19}
$$

and

$$\hat{\psi}\_{\*}^{\alpha}(\varpi) = \mathcal{C} \,(j\varpi)^{\alpha+1}, \,\, as \,\, \varpi \to 0 \tag{20}$$

where ψˆ α + (̟) is symmetric.

# STRATEGIES IN CWT APPLICATIONS TO UV SPECTROSCOPY ANALYSIS OF MULTICOMPONENT MIXTURES

For the past 15 years, the potential application of CWT in chemistry, especially in combination with other mathematical methods, leads us to a conclusion that WT has interestingly became a useful algorithm for UV quantitative analysis of pharmaceuticals. Four different models [i.e., continuous wavelet transform-zero crossing (CWT-ZC), ratio spectra-continuous wavelet transform (RS-CWT), ratio spectra-continuous wavelet transform-zero crossing (RS-CWT-ZC), and double divisor ratio spectra-continuous wavelet transform (DDRS-CWT)] were described in the implementation of CWT to UV spectroscopic data for the resolution of overlapping spectra to quantify drugs in different types of samples. The modeling of CWT—UV spectroscopic approaches are detailed below. Fundamentally, these approached can be successfully applied to the UV spectroscopic analysis of binary and ternary mixtures, provided that the law of additivity of absorbance is obeyed.

# CONTINUOUS WAVELET TRANSFORM-ZERO CROSSING

The application of CWT-ZC approach to UV spectroscopic signals was first proposed by Dinç and Baleanu (2003a).

If a mixture of two analytes (M and N) is considered (see **Figure 2A**) and the absorbance of this binary mixture is measured at λi , we can have the following equation (Charlotte Grinter and Threlfall, 1992):

$$\mathbf{A}\_{\rm mix,\,\lambda i} = \alpha\_{\mathbf{M},\,\lambda i} \mathbf{C}\_{\mathbf{M}} + \beta\_{\mathbf{N},\,\lambda i} \mathbf{C}\_{\mathbf{N}} \tag{21}$$

where Amλi is the absorbance of the binary mixture at wavelength λi , and the coefficients are the absorptivities of M and N, respectively. CM and CN represent the concentrations of M and N, respectively.

If CWT is applied to Equation (21), the following function can be obtained as

$$
\psi\_{\text{(a.b),MIX, \lambda i}} = \psi\_{\text{(a.b), M, \lambda i}} \mathcal{C}\_{\text{M}} + \psi\_{\text{(a.b), N, \lambda i}} \mathcal{C}\_{\text{N}} \tag{22}
$$

If ψ(a.b),N,λiCN = 0, then we obtain the following equation

$$
\Psi\_{\text{(a.b)},\text{MIX},\lambda i} = \Psi\_{\text{(a.b)},\text{M},\lambda i} \,\text{C}\_{\text{M}} \tag{23}
$$

Equation (23) shows that CWT (ψ(a.b), M, λi CM) amplitudes of M in the binary mixture are dependent only on CM regardless of CN (see **Figure 2B**).

# RATIO SPECTRA-CONTINUOUS WAVELET TRANSFORM

Apart from CWT-ZC approach, overlapping spectral bands in a binary mixture could be solved by the application of a combined hybrid approach i.e., RS-CWT (Dinç and Baleanu, 2004a,c).

The absorption spectra of M and N compounds, and their mixture are indicated in **Figure 3A**. By being divided by the standard spectrum (AN,λi = βλi C o N ) of one of the compounds in the binary mixture, Equation (21) becomes

$$\frac{A\_{\text{m},\lambda\dot{i}}}{\beta\_{\lambda\dot{i}}\,\,\mathbf{C}\_{\text{N}}\,^{o}} = \frac{\alpha\_{\lambda\dot{i}}\,\,\mathbf{C}\_{\text{M}}}{\beta\_{\lambda\dot{i}}\,\,\mathbf{C}\_{\text{N}}\,^{o}} + \frac{\beta\_{\lambda\dot{i}}\,\,\mathbf{C}\_{\text{N}}}{\beta\_{\lambda\dot{i}}\,\,\mathbf{C}\_{\text{N}}\,^{o}}\tag{24}$$

**Figure 3B** shows the ratio spectra of analytes and their binary mixture. If CWT is applied to Equation (24), the following equation can be obtained

$$CWT\left[\frac{A\_{m,\lambda i}}{\beta\_{\lambda i} \, ^{\circ}C\_{N}^{\circ}}\right] = CWT\left[\frac{\alpha\_{\lambda i}}{\beta\_{\lambda i}}\right] \frac{^{\circ}C\_{M}}{C\_{N}^{\circ}} + CWT\left[\frac{\beta\_{\lambda i}}{\beta\_{\lambda i}}\right] \frac{^{\circ}C\_{N}}{^{\circ}C\_{N}^{\circ}} \tag{25}$$

If CWT h βλi βλi i CN C o N = 0, then we obtain

$$\text{CWT} \left[ \frac{\text{A}\_{\text{m}, \text{\AA}}}{\beta\_{\lambda i} \text{ C}\_{\text{N}}^{\text{o}}} \right] = \text{CWT} \left[ \frac{\alpha\_{\lambda i}}{\beta\_{\lambda i}} \right] \frac{\text{C}\_{\text{M}}}{\text{C}\_{\text{N}}^{\text{o}}} \tag{26}$$

The ratio-CWT amplitudes of the binary mixture given in Equation (26) depend only on CM and CN o regardless of CN (e.g., see **Figure 3C**).

# RATIO SPECTRA-CONTINUOUS WAVELET TRANSFORM-ZERO CROSSING

In RS-CWT-ZC approach (Dinç et al., 2005a), if a mixture of three analytes (X, Y, and Z) is considered and the absorbance of this ternary mixture is measured at λi , the following mathematical expression (Charlotte Grinter and Threlfall, 1992) would be given

$$\mathbf{A}\_{\text{mix},\lambda i} = \alpha\_{\text{X},\lambda i} \mathbf{C}\_{\text{X}} + \beta\_{\text{Y},\lambda i} \mathbf{C}\_{\text{Y}} + \gamma\_{\text{Z},\lambda i} \mathbf{C}\_{\text{Z}} \tag{27}$$

Where Amix, λi is the absorbance of the ternary mixture at wavelength λi , and coefficients αX, λi , βY, λi , and γZ, λi denote the absorptivities of X, Y, and Z, respectively. CX, CY, and CZ represent the concentrations of X, Y, and Z, respectively.

If Equation (27) is divided by the spectrum of a standard solution (Co X ) of one of the compounds in the ternary mixture, we have the following equation:

$$\frac{\mathbf{A}\_{\text{mix, \lambda i}}}{\alpha \mathbf{x}\_{\text{, \lambda i}} \mathbf{C} \mathbf{x}^{o}} = \frac{\alpha\_{\text{X, \lambda i}} \mathbf{C}\_{\text{X}}}{\alpha \mathbf{x}\_{\text{, \lambda i}} \mathbf{C} \mathbf{x}^{o}} + \frac{\beta\_{\text{Y, \lambda i}} \mathbf{C}\_{\text{Y}}}{\alpha \mathbf{x}\_{\text{, \lambda i}} \mathbf{C}\_{\text{X}}^{o}} + \frac{\gamma\_{\text{Z, \lambda i}} \mathbf{C}\_{\text{Z}}}{\alpha \mathbf{x}\_{\text{, \lambda i}} \mathbf{C}\_{\text{X}}^{o}} \tag{28}$$

If CWT is applied to Equation (28), the following equation can be obtained

$$\text{CWT} \left[ \frac{\text{A}\_{\text{mix},\lambda i}}{\alpha\_{\text{X},\lambda i} \text{C}\_{\text{X}}^{\text{o}}} \right] = \text{CWT} \left[ \frac{\beta\_{\text{Y},\lambda i} \text{C}\_{\text{Y}}}{\alpha\_{\text{X},\lambda i} \text{C}\_{\text{X}}^{\text{o}}} \right] + \text{CWT} \left[ \frac{\nu\_{\text{Z},\lambda i} \text{C}\_{\text{Z}}}{\alpha\_{\text{X},\lambda i} \text{C}\_{\text{X}}^{\text{o}}} \right] \tag{29}$$

Equation (29) indicates that the CWT amplitudes of the ratio spectra of the ternary mixture are dependent only on CZ and CX o regardless of the concentrations of other compounds.

# DOUBLE DIVISOR RATIO SPECTRA-CONTINUOUS WAVELET TRANSFORM

In addition to RS-CWT-ZC approach, the spectral resolution of ternary mixtures could be effectively done by DDRS-CWT approach (Dinç and Baleanu, 2008a) as follows.

When two compounds in the ternary mixture is used as a double divisor, we have

$$\mathcal{A}^{\bullet}\_{\text{mix, \lambda i}} = \mathcal{a}\_{\text{X, \lambda i}} \mathcal{C}^{\bullet}\_{\text{X}} + \beta\_{\text{Y, \lambda i}} \mathcal{C}^{\bullet}\_{\text{Y}} \tag{30}$$

By dividing Equation (27) and (30), we obtain as follows

$$\begin{split} \frac{\mathbf{A}\_{\text{mix, \\_\text{i}\\_i}}}{\alpha\_{\text{X, \\_i\\_i}}\mathbf{C}\_{\text{X}}^{\text{o}} + \beta\_{\text{Y, \\_i\\_i}}\mathbf{C}\_{\text{Y}}^{\text{o}}} &= \frac{\alpha\_{\text{X, \\_i\\_i}}\mathbf{C}\_{\text{X}}}{\alpha\_{\text{X, \\_i\\_i}}\mathbf{C}\_{\text{X}}^{\text{o}} + \beta\_{\text{Y, \\_i\\_i}}\mathbf{C}\_{\text{Y}}^{\text{o}}} \\ + \frac{\beta\_{\text{Y, \\_i\\_i}}\mathbf{C}\_{\text{Y}}^{\text{o}}}{\alpha\_{\text{X, \\_i\\_i}}\mathbf{C}\_{\text{X}}^{\text{o}} + \beta\_{\text{Y, \\_i\\_i}}\mathbf{C}\_{\text{Y}}^{\text{o}}} + \frac{\gamma\_{\text{Z, \\_i\\_i}}\mathbf{C}\_{\text{Z}}^{\text{o}}}{\alpha\_{\text{X, \\_i\\_i}}\mathbf{C}\_{\text{X}}^{\text{o}} + \beta\_{\text{Y, \\_i\\_i}}\mathbf{C}\_{\text{Y}}^{\text{o}}} \end{split} \tag{31}$$

TABLE 2 | Applications of the continuous wavelet transform-zero crossing technique to UV spectroscopic analysis of pharmaceuticals.

Equation (31) can be simplified to

$$\frac{\mathbf{A}\_{\text{mix, \lambda i}}}{\alpha\_{\text{X, \lambda i}} \mathbf{C}\_{\text{X}} \mathbf{o}^{\text{o}} + \beta\_{\text{Y, \lambda i}} \mathbf{C}\_{\text{Y}} \mathbf{o}^{\text{o}}} = \mathbf{k} + \frac{\gamma\_{\text{Z, \lambda i}} \mathbf{C}\_{\text{Z}}}{\alpha\_{\text{X, \lambda i}} \mathbf{C}\_{\text{X}} \mathbf{o}^{\text{o}} + \beta\_{\text{Y, \lambda i}} \mathbf{C}\_{\text{Y}} \mathbf{o}^{\text{o}}} \tag{32}$$

Where k = αX, λiCX+ βY, λiCY αX, λiC o X+ βY, λi C o Y represents a constant for a given concentration range with respect to λi in a certain region or point of wavelength.

A typical case is when CX o and CY o are the same or very close to each other, namely CX o = CY o or ∼= CX o ∼= CY o . Therefore, we obtain

$$
\alpha\_{\rm X, \lambda \dot{\imath}} \mathcal{C}\_{\rm X}^{\bullet} + \beta\_{\rm Y, \lambda \dot{\imath}} \mathcal{C}\_{\rm Y}^{\bullet} = \mathcal{C}\_{\rm X}^{\bullet} \left( \alpha\_{\rm X, \lambda \dot{\imath}} + \beta\_{\rm Y, \lambda \dot{\imath}} \right) \tag{33}
$$

and Equation (32) can be written as

$$\frac{\mathbf{A}\_{\text{mix, \lambda i}}}{\alpha\_{\text{X, \lambda i}} \mathbf{C}\_{\text{X}}^{\text{o}} + \beta\_{\text{Y, \lambda i}} \mathbf{C}\_{\text{Y}}^{\text{o}}} = k + \frac{\mathbf{y}\_{\text{Z, \lambda i}} \mathbf{C}\_{\text{Z}}}{\mathbf{C}\_{\text{X}}^{\text{o}} \left(\alpha\_{\text{X, \lambda i}} + \beta\_{\text{Y, \lambda i}}\right)} \tag{34}$$

After applying CWT to Equation (31), we have

$$\text{CWT}\_{\text{(a,b)}} \left( \frac{\text{A}\_{\text{mix, \lambda i}}}{\text{a} \chi\_{\text{, \lambda i}} + \text{\beta} \chi\_{\text{, \lambda i}}} \right) \frac{1}{\text{C}\_{\text{X}}^{\text{o}}} = \text{CWT}\_{\text{(a,b)}} \left( \frac{\chi\_{\text{Z, \lambda i}} \text{C}\_{\text{Z}}}{\left( \alpha\_{\text{X, \lambda i}} + \beta\_{\text{Y, \lambda i}} \right)} \right) \frac{1}{\text{C}\_{\text{X}}^{\text{o}}} \tag{35}$$

or

$$\text{CWT}\_{(a,b)}\left(\frac{\text{A mix},\text{\AA}}{\text{a}\text{X},\text{\AA}+\text{\beta}\text{Y},\text{\AA}}\right) = \text{CWT}\_{(a,b)}\left(\frac{\text{YZ},\text{\AA}}{\left(\text{a}\_{\text{X},\text{\AA}\text{i}}+\text{\beta}\text{Y}\_{\text{\text},\text{\AA}\text{i}}\right)}\right)\text{CZ}\tag{36}$$

In Equation (36), CZ is to proportional to the coefficients, CWT(a,b) Amix, λi αX, λi+ βY, λi , at λi. If this procedure is separately applied for pure Z and its ternary mixture, the CWT(a,b) coefficients are coincided at some characteristic point or region of wavelength, independent upon both CX and CY.

# WAVELET TRANSFORM-BASED UV SPECTROSCOPIC ANALYSIS OF PHARMACEUTICALS

Typical applications of CWT and FWT algorithms for UV spectroscopic analysis of pharmaceuticals are displayed in **Tables 2**–**5**. It is worth mentioning that WT could be solely applied to raw spectra and ratio spectra (as above-specified) as well as utilized as a hybrid approach (FWT-derivative, FWT-CWT-zero crossing, WT combined with multivariate calibration) for the simultaneous determination of analytes in pharmaceutical binary and ternary mixtures. It was shown that wavelet analysis of UV spectroscopic data was performed by using Wavelet Toolbox and m-file in MATLAB software. The numerous works provided by Dinç and co-workers have clearly highlighted the success of WT-based UV spectroscopic analysis for multicomponent synthetic mixtures, veterinary and pharmaceutical dosage forms as well as different types of test (e.g., assay, in vitro dissolution, stability indicating). Most studies proved it to be suitable for the routine analysis of dosage forms with good precision and accuracy, comparable to HPLC.

# CONCLUSIONS

In the point of view of UV spectroscopic analysis of multicomponent mixtures, CWT-based UV spectroscopic

TABLE 4 | Applications of the ratio spectra-continuous wavelet transform, ratio spectra- continuous wavelet transform-zero crossing approaches to UV spectroscopic analysis of pharmaceuticals.

methods have outperformed both conventional and derivative UV spectroscopy in resolving spectrally binary and ternary mixtures. Nevertheless, wavelet analysis may not also have a sufficient power to resolve overlapping spectra of analytes in samples due to similarity of molecular structures and signal frequencies in some cases. They may not give desirable results for a complex mixture containing more than three compounds and/or a significant difference in ratios of active ingredients. In such a case, the use of WT coupled with chemometric PLS and PCR calibrations is advisable. Undoubtedly, however, wavelets can still be used as

## REFERENCES

a mathematical prism for signal analysis because they can offer many possibilities such as baseline correction, noise removal and resolution of overlapping peaks, when the frequencies of analyzed components are significantly different from each other.

# AUTHOR CONTRIBUTIONS

Contributions of ED are planning and writing of the review paper. Contributions of ZY are literature review, collection, editing, and format arrangement.

of the fractional fourier transform bands for resolving two component mixture. Signal Image Video P. 9, 801–807. doi: 10.1007/s11760-013-0503-9

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Dinç and Yazan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Chemometric Methods for Spectroscopy-Based Pharmaceutical Analysis

Alessandra Biancolillo and Federico Marini\*

*Department of Chemistry, University of Rome La Sapienza, Rome, Italy*

Spectroscopy is widely used to characterize pharmaceutical products or processes, especially due to its desirable characteristics of being rapid, cheap, non-invasive/non-destructive and applicable both off-line and in-/at-/on-line. Spectroscopic techniques produce profiles containing a high amount of information, which can profitably be exploited through the use of multivariate mathematic and statistic (chemometric) techniques. The present paper aims at providing a brief overview of the different chemometric approaches applicable in the context of spectroscopy-based pharmaceutical analysis, discussing both the unsupervised exploration of the collected data and the possibility of building predictive models for both quantitative (calibration) and qualitative (classification) responses.

### Edited by:

*Cosimino Malitesta, University of Salento, Italy*

### Reviewed by:

*Daniel Cozzolino, Central Queensland University, Australia Andreia Michelle Smith-Moritz, University of California, Davis, United States*

\*Correspondence: *Federico Marini<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Analytical Chemistry, a section of the journal Frontiers in Chemistry*

Received: *07 July 2018* Accepted: *05 November 2018* Published: *21 November 2018*

### Citation:

*Biancolillo A and Marini F (2018) Chemometric Methods for Spectroscopy-Based Pharmaceutical Analysis. Front. Chem. 6:576. doi: 10.3389/fchem.2018.00576* Keywords: spectroscopy, chemometrics and statistics, component analysis (PCA), partial least squares (PLS), classification, partial least squares discriminant analysis (PLS-DA), soft independent modeling of class analogies (SIMCA), pharmaceutical quality control

# INTRODUCTION

Quality control on pharmaceutical products is undoubtedly an important and widely debated topic. Hence, in literature, various methods have been proposed to check quality of medicines, either qualitative (e.g., for the identification of an active pharmaceutical ingredient, API; Blanco et al., 2000; Herkert et al., 2001; Alvarenga et al., 2008) or quantitative (quantification of the API; Blanco et al., 2000; Yao et al., 2007; Cruz Sarraguça and Almeida Lopes, 2009); involving either destructive or non-invasive online techniques. Recently, due to the benefits they bring, several non-destructive methodologies based on spectroscopic techniques (mainly Near-Infrared NIR) combined with chemometric tools have been proposed for pharmaceutical quality check (Chen et al., 2018; Rodionova et al., 2018).

Despite the development of analytical methodologies and the commitments of national and supranational entities to regulate pharmaceutical quality control, substandard and counterfeit medicines are still a major problem all over the world.

# Chemometrics as Tool for Fraud/Adulteration Detection

Poor-quality pharmaceuticals can be found on the market for two main reasons: low production standards (mainly leading to substandard medicines) and fraud attempts. Counterfeited drugs may present different frauds/adulterations; for instance, they could contain no active pharmaceutical ingredient (API), a different API from the one declared, or a different (lower) API strength. As mentioned above, several methodologies have been proposed in order to detect substandard/counterfeit pharmaceuticals; among these, a major role is played by those based on the application of spectroscopic techniques in combination with different chemometric methods. The relevance of these methodologies is due to the fact that spectroscopy (in particular, NIR) combined with exploratory data analysis, classification and regression method can lead to effective, high performing, fast, nondestructive, and sometimes, online methods for checking the quality of pharmaceuticals and their compliance to production and/or pharmacopeia standards. Nevertheless, the available chemometric tools applicable to handle spectroscopic (but, of course not only those) data are numerous, and there is plenty of room for their misapplication (Kjeldahl and Bro, 2010). As a consequence, the aim of the present paper is to report and critically discuss some of the chemometric methods typically applied for pharmaceutical analysis, together with an essential description of the figures of merit which allow evaluating the quality of the corresponding models.

# EXPLORATORY DATA ANALYSIS

In the large part of the studies for the characterization of pharmaceutical samples for quality control, verification of compliance and identification/detection of counterfeit, fraud or adulterations, experimental signals (usually in the form of some sorts of fingerprints) are collected on a series of specimens. These constitute the data the chemometric models operate on. These data are usually arranged in the form of a matrix **X**, having as many rows as the number of samples and as many columns as the number of measured variables. Accordingly, assuming that samples are spectroscopically characterized by collecting an absorption (or reflection/transmission) profile (e.g., in the infrared region), each row of the matrix corresponds to the whole spectrum of a particular sample, whereas each column represents the absorbance (or reflectance/transmittance) of all the individuals at a particular wavenumber. This equivalence between the experimental profiles and their matrix representation is graphically reported in **Figure 1**.

Once the data have been collected, exploratory data analysis represents the first step of any chemometric processing, as it allows "to summarize the main characteristics of data in an easy-to-understand form, often with visual graphs, without using a statistical model or having formulated a hypothesis" (Tukey, 1977). Exploratory data analysis provides an overall view of the system under study, allowing to catch possible similarities/dissimilarities among samples, to identify the presence of clusters or, in general, systematic trends, to discover which variables are relevant to describe the system and, on the other hand, which could be in principle discarded, and to detect possible outlying, anomalous or, at least, suspicious samples (if present). As evident also from the definition reported above, in the context of exploratory data analysis a key role is played by the possibility of capturing the main structure of the data in a series of representative plots, through appropriate display techniques. Indeed, considering a general data matrix **X**, of dimensions N×M, one could think of its entries as the coordinates of N points (the samples) into a Mdimensional space whose axes are the variables, which makes this representation unfeasible for the cases when more than three descriptors are collected on each individual. This is why exploratory data analysis often relies on the use of projection (bilinear) techniques to reduce the data dimensionality in a "clever" way. Projection methods look for a low-dimensional representation of the data, whose axes (normally deemed components or latent variables) are as relevant as possible for the specific task. In the case of exploratory data analysis, the most commonly used technique is Principal Components Analysis (PCA) (Pearson, 1901; Wold et al., 1987; Jolliffe, 2002).

# Principal Component Analysis

Principal component analysis (PCA) is a projection method, which looks for directions in the multivariate space progressively providing the best fit of the data distribution, i.e., which best approximate the data in a least squares sense. This explains why PCA is the technique of choice in the majority of cases when exploratory data analysis is the task: indeed, by definition, for any desired number of dimensions (components) F in the final representation, the subspace identified by PCA constitutes the most faithful F-dimensional approximation of the original data. This allows compression of the data dimensionality at the same time reducing to a minimum the loss of information. In particular, starting from a data matrix **X**(N×M) , Principal Component Analysis is based on its bilinear decomposition, which can be mathematically described by Equation (1):

$$X = T\mathbf{P}^T + E \tag{1}$$

The loadings matrix **P**(M×F) identifies the F directions, i.e., the principal components (PC), along which the data should be projected and the results of such projection, i.e., the coordinates of the samples onto this reduced subspace, are collected in the scores matrix **T**(N×F) . In order to achieve data compression, usually F ≪ M so that the PCA representation provides an approximation of the original data whose residuals are collected in the matrix **E**(N×M) .

Since the scores represent a new set of coordinates along highly informative (relevant) directions, they may be used in two- or three-dimensional scatterplots (scores plots). This offers a straightforward visualization of the data, which can highlight possible trends in data, presence of clusters or, in general, of an underlying structure. A schematic representation of how PCA works is displayed in **Figure 2**.

**Figure 2** shows one of the simplest possible examples of feature reduction, since it describes the case where samples described by three measured variables can be approximated by being projected on an appropriately chosen two-dimensional sub-space. However, the concept may be easily generalized to higher-dimensional problems, such as those involving spectroscopic measurements. **Figure 3** shows an example of the application of PCA to mid infrared spectroscopic data. In particular, the possibility of extracting as much information as possible from the IR spectra recorded on 51 tablets containing

either ketoprofen or ibuprofen in the region 2,000–680 cm−1 (661 variables) is represented.

A large portion of the data variability can be summarized by projecting the samples onto the space spanned by the first two principal components, which account for about 90% of the original variance, and therefore can be considered as a good approximation of the experimental matrix. Inspection of the scores plot suggests that the main source of variability is the difference between ibuprofen tablets (blue squares) and ketoprofen ones (red circles), since the two clusters are completely separated along the first principal component. To interpret the observed cluster structure in terms of the measured variables, it is then necessary to inspect the corresponding loadings, which are also displayed in **Figure 3** for PC1. Indeed, for spectral data, the possibility of plotting the loadings for the individual components in a profile-like fashion, rather than producing scatterplot for pairs of latent variables (as exemplified in **Figure 2**) is often preferred, due to its more straightforward interpretability: spectral regions having positive loadings will have higher intensity on samples which have positive scores on the corresponding component, whereas bands associated to negative loadings will present higher intensity on the individuals falling at negative values of the PC. In the example reported in **Figure 3**, one could infer, for instance, that the ketoprofen samples (which fall at positive values of PC1) have a higher absorbance at the wavenumbers where the loadings are positive, whereas ibuprofen samples should present a higher signal in correspondence to the bands showing negative loadings.

Based on what reported above, it is evident how the quality of the compressed representation in the PC space depends on the number of components F chosen to describe the data. However, at the same time, it must be noted that when the aim of calculating a PCA is "only" data display, as in most of the applications in the context of exploratory analysis, the choice of the optimal number of components is not critical: it is normally enough to inspect the data distribution across the first few dimensions and, in many cases, considering the scores plot resulting from the first two or three components could be sufficient. On the other hand, there may be cases when the aim of the exploratory analysis is not limited to just data visualization and, for instance, one is interested in the identification of anomalous or outlying observations, or there could be the need of the imputation of missing elements in the data matrix; additionally, one could also need to obtain a compressed representation of the data to be used for further predictive modeling. In all such cases, the choice of the optimal dimensionality of the PC representation is critical for the specific purposes and, therefore, the number of PCs should be carefully estimated. In this respect, different methods have been proposed in the literature and a survey of the most commonly used can be found in Jolliffe (2002).

Among the applications described above, the possibility of using PCA for the identification/detection of potential outliers deserves a few more words, as it could be of interest for pharmaceutical quality control. Actually, although outliers—or anomalous observations, in general—could be, in principle, investigated by visually inspecting the scores plot along the first components, this approach could be subjective and anyway would not consider some possible data discrepancies. Alternatively, when it is used as a model to

(highlighted in light red in the leftmost panel) spanned by the first two principal components. Inspection of the data set can be carried out by looking at the distribution of the samples onto the informative PC subspace (scores plot) and interpretation can be then carried out by examining the relative contribution of the experimental variable to the definition of the principal components (loadings plot).

build a suitable approximation of the data, PCA provides a powerful toolbox for outlier detection based on the definition of more objective test statistics, which can be easily automatized or, anyway, embedded in control strategies, also on-line. This is accomplished by defining two distance measurements: (i) a squared Mahalanobis distance in the scores space, which follows the T2 statistics (Hotelling, 1931) and accounts for how extreme the measurement is in the principal component subspace, and (ii) a squared orthogonal Euclidean distance (the sum of squares of the residuals after approximating the observation by its projection), which is normally indicated as Q statistics (Jackson and Muldholkar, 1979) and quantifies how well the model fits that particular individual. Outlier detection is then carried out by setting appropriate threshold values for the T2 and Q statistics and verifying whether the samples fall below or above those critical limits. Moreover, once an observation is identified as a potential outlier, inspection of the contribution plot can help in relating the detected anomaly to the behavior of specific measured variables.

## Selected Examples

PCA is customarily used for the quality control of drugs and pharmaceuticals; several examples of the application of this technique to solve diverse issues have been reported in the literature.  For example, in Rodionova et al. (2005) PCA was applied to both bulk NIR spectroscopy and hyperspectral imaging (HSI) in the NIR range to spot counterfeit drugs. In particular, bulk NIR was used to differentiate genuine antispasmodic drugs from forgeries, whereas HSI on the ground uncoated tablets was employed to identify counterfeited antimicrobial drugs. In both cases, the spectroscopic data were subjected to PCA, which allowed to clearly identify clusters in the scores plot, corresponding to the two kinds of tablets, i.e., genuine and counterfeited. In the case of the imaging platform, where the signal is stored as a data hypercube [i.e., a three-way numerical array of dimension number of horizontal pixels Nx, number of vertical pixels Ny and number of wavelengths Nλ, in which each entry corresponds to the spectral intensity measured

observed differences in terms of the spectroscopic signal is made possible by the inspection of the loadings on PC, which are shown in a "spectral-like" fashion in (C).

at a certain wavelength and a specific spatial position (x-y coordinates)], a preliminary unfolding step is needed. Unfolding is the procedure allowing to reorganize a higher-order array into a two-way matrix, which can be then processed with standard chemometric techniques. In the case of hyperspectral data cubes, this is carried out by stacking the spectra corresponding to the different pixels one on top of each other, in a way to obtain a matrix of dimensions (Nx × Ny and Nλ).

Another relevant application of exploratory analysis is related to quality check. For instance, PCA can be applied to investigate formulations not meeting predefined parameters. In Roggo et al. (2005), PCA was used to inquire a suspicious blue spot present on tablets. Samples were analyzed by a multi-spectral (IR) imaging microscope and PCA analysis was performed on the unfolded data-cube, indicating that the localized coloration was not due to contamination, but actually given by wet indigo carmine dye and placebo (expected ingredients of the formulation).

PCA can also be used for routine quality checks at the end of a production process. For example, in Myakalwar et al. (2011) laser-induced breakdown spectroscopy (LIBS) and PCA were combined with the aim of obtaining qualitative information about the composition of different pharmaceuticals.

# REGRESSION

As discussed in the previous section, exploratory analysis is a first and fundamental step in chemometric data processing and, in some cases, it could be the only approach needed to characterize the samples under investigation. However, due to its unsupervised nature, it provides only a (hopefully) unbiased picture of the data distribution but it lacks any possibility of formulating predictions on new observations, which on the other hand may be a fundamental aspect to solve specific issues. In practice, very often quality control and/or authentication of pharmaceutical products rely on some forms of qualitative or quantitative predictions. For instance, the quantification of a specific compound (e.g., an active ingredient or an excipient) contained in a formulation is a routine operation in pharmaceutical laboratories. This goal can be achieved by combining instrumental (e.g. spectroscopic) measurements with chemometric regression approaches (Martens and Naes, 1991; Martens and Geladi, 2004). Indeed, given a response to be predicted y and a vector of measured signals (e.g., a spectrum) **x**, the aim of regression methods is to find the functional relationship that best approximates the response on the basis of the measurements (the predictors). Mathematically, this can be stated as:

$$
\gamma = \hat{\boldsymbol{\gamma}} + \boldsymbol{e} = \boldsymbol{f} \; (\boldsymbol{\mathfrak{x}}) + \boldsymbol{e} \tag{2}
$$

where yˆ is the predicted response (i.e., the response value approximated by the model), f (**x**) indicates a general function of and **x** and e is the residual, i.e., the difference between the actual response and its predicted value. In many applications, the functional relationship between the response and the predictors f (**x**) can be assumed to be linear:

$$\hat{\boldsymbol{y}} = \boldsymbol{f}\left(\mathbf{x}\right) = b\_1\boldsymbol{\chi}\_1 + b\_2\boldsymbol{\chi}\_2 + \dots + b\_M\boldsymbol{\chi}\_M = \mathbf{x}^T\boldsymbol{\mathbf{b}}\tag{3}$$

where x1, x2 . . . xM are the components of the vector of measurements **x** and the transpose indicates that it is normally expressed as a row vector, while the associated linear coefficients b1, b2 . . . bM, which weight the contributions of each of the M Xvariables to y, are called regression coefficients and collected in the vector **b**. Building a regression model means to find the optimal value of the parameters **b**, i.e., the values which lead to the lowest error in the prediction of the responses. As a direct consequence of this consideration, it is obvious how it is mandatory to have a set of samples (the so-called training set) for which both the experimental data **X** and the responses **y** are available, in order to build a predictive model. Indeed, the information on the **y** is actively used to calculate the model parameters. When data from more than a single sample are available, the regression problem in Equations (2, 3) can be reformulated as:

$$\mathbf{y} = \hat{\mathbf{y}} + \mathbf{e} = \mathbf{X}b + \mathbf{e} \tag{4}$$

where the vectors **y**ˆ and **e** collect the predictions and residuals for the different samples, respectively. Accordingly, the most straightforward way of calculating the model parameters in Equation (4) is by the ordinary least-squares approach, i.e., by looking at those values of **b**, which minimize the sum of squares of the residuals **e**:

$$\min\_{\mathbf{b}} \mathbf{e}^T \mathbf{e} = \min\_{\mathbf{b}} \sum\_{i=1}^N e\_i^2 \tag{5}$$

ei being the residual for the ith sample and N being the number of training observations. The corresponding methods is called multiple linear regression (MLR) and, under the conditions of Equation (5), the regression coefficients are calculated as:

$$b = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} \tag{6}$$

Here it is worth to highlight that, if one wishes to use the same experimental matrix **X** to predict more than one response, i.e., if, for each sample, instead of a single scalar yi , there is a dependent vector

$$\mathbf{y}\_i^T = \begin{bmatrix} \mathbf{y}\_{i1}\mathbf{y}\_{i2} \cdots \mathbf{y}\_{iL} \end{bmatrix} \tag{7}$$

L being the number of responses, then each dependent variable should be regressed on the independent block by means of a set of regression coefficients. Assuming that the L responses measured on the training samples are collected in a matrix **Y**, whose columns **y**l are the individual dependent variables,