Hybrid Denosing and Trend Removal Techniques for Modulated Spectroscopy

The application of a Trend Correction Method (TCM) becomes imperative for accurate analysis of spectra. TCM algorithms can be categorized into two distinct groups based on the nature of the information extracted from the original signals. The first group encompasses algorithms that require an understanding of the background, the blurring effect, and noise. These algorithms typically engage with signals by utilizing knowledge about signal components such as background shape, position, and Signal-to-Noise Ratio (SNR). Included in this category are methods like the median filtering method, Signal Removal Method (SRM), and Threshold-Based Classification (TBC) [5, 6].

The second group of TCMs is predicated on knowledge about the frequencies of signal components. It is well understood that noise and background exhibit distinct characteristics: noise generally manifests as a high-frequency phenomenon, whereas the background is characterized by low-frequency components within the signal. Signal processing techniques in this category include the Fourier Transform (FT) [7] and the Wavelet Transform (WT) techniques [8].

The challenge of background (trend) interference in spectroscopic analysis can be addressed through various algorithms. A notable approach combines the Signal Removal Method (SRM) and the Continuous Wavelet Transform (CWT) technique, as developed by Kandjani et al. [9]. In this approach, the initial step involves utilization of CWT to compute the second derivative of the noisy signal, thereby identifying signal peak positions in the experimental spectrum. Subsequently, the SRM is applied to discern the trend, which is then subtracted from the original spectrum. This process effectively corrects the signal by removing peak components of the spectrum and fitting the remaining spectrum.

Another significant contribution by Schulze et al. [10] introduced an automated baseline-removal method for Raman spectra, specifically addressing sloping baselines. This method employs an Automated and model-free Baseline Estimation (ABE) procedure, remarkably requiring minimal human intervention. The efficacy of ABE is monitored using the v2-statistic, which serves as a computational stopping condition based on a zero-order Savitzky–Golay filter. However, this approach has limitations, particularly in correcting baselines with significant curvatures. This often necessitates smaller filter windows, leading to peak erosion, and in some cases, the occurrence of undesired local minima at a maximum window size, resulting in incomplete baseline removal. Additionally, edge artifacts present a challenge.

In the realm of noise reduction, spatial filters like mean and median filters are common. Nonetheless, these filters tend to smooth data, which may be a drawback in preserving signal integrity. The wavelet transform, known for its multi-resolution capabilities, emerges as a superior tool in signal and image processing. Its advantages, such as sparsity and multi-resolution structure, make wavelet denoising more effective than other image denoising techniques. Over the past two decades, as the wavelet transform gained popularity, various denoising techniques in the wavelet domain were introduced, leading to the preference of the wavelet transform domain over time and frequency domains [11].

4.1 Proposed algorithm for removing trend

Initially, to simulate Raman spectral peaks, we utilized a Gaussian function, defined as follows:

$f(x)=a \cdot e^{\left(\frac{-0.5(x-c)^2}{\sigma^2}\right)}$    (2)

where, a represents the intensity, while c and σ denote the mean and variance of the Gaussian peak, respectively. Subsequently, we introduced three distinct types of trends for analysis. Figure 2 shows the original spectrum.

Figure 2. Original spectrum

4.1.1 Linear trend

The first type, namely linear trend, is mathematically expressed as:

Trend $=a \cdot x+b$    (3)

In this equation, a signifies the slope of the linear trend, and b is a constant. The graphical representation of this linear trend is illustrated in Figure 3(a).

4.1.2 Sigmoidal trend

The next type of trend analyzed in our study is the sigmoidal trend, defined by the following equation:

Trend $=\frac{1}{1+\exp (-a(x-c))} I+O$      (4)

where, a determines the gradient at the inflection point, c specifies the location of the inflection point, I acts as an intensity controller, and O represents the offset. This trend is graphically illustrated in Figure 3(b).

(a) Linear trend

(b) Sigmoidal trend

(c) Sinusoidal trend

Figure 3. Spectrum with a) Linear trend, b) Sigmoidal trend, and c) Sinusoidal trend

4.1.3 Sinusoidal trend

The third type is the sinusoidal trend, which mimics the shape of a sinusoidal wave. This form is frequently encountered, and it is mathematically expressed as:

Trend $=x^{1.5} \sin \left(\frac{x}{a}\right) \cdot I+O$ (5)

where, a modulates the frequency, while I and O control the intensity and offset, respectively. Noise in our study is considered as white Gaussian noise, added in accordance with the calculated SNR in decibels. Figures 3(a)-(c) collectively depict the three types of trends assumed in our analysis.

In our methodology, the estimation of the trend involves a two-step process. First, we apply local variance estimation to accurately identify the peak locations. Variance, as a statistical measure, quantifies the degree of spread in a set of data points. A variance of zero indicates uniformity, with all values being identical. Conversely, a small variance value indicates that the data points are closely clustered around the mean or expected value, revealing minimal spread. On the other hand, a high variance value implies a significant dispersion of data around the mean, indicating a wider spread.

Upon determining the peak locations, the next step involves cubic spline interpolation. This technique is employed to effectively reconstruct the data points in the regions where peaks were identified and subsequently removed. By doing so, we ensure a smooth transition and continuity in the data, facilitating a more accurate trend estimation. The application of this two-step approach in trend estimation and elimination is illustrated in Figure 4.

Figure 4. Variance of data segment (x)

The process of estimating the local variance for a signal, denoted as x(n), is a crucial aspect of our methodology. This estimation is articulated through the following mathematical expression [12]:

$\hat{\sigma}_x^2(n)=\frac{1}{(2 K+1)} \sum_{k=n-K}^{n+K}(x(k)-\hat{x}(n))^2$    (6)

where, the term $(2 K+1)$ represents the number of samples within a short segment, which is utilized for variance estimation. Additionally, the local mean of the signal, represented as $\hat{x}(n)$, plays a pivotal role in this calculation. The expression for the local mean is defined as:

$\hat{x}(n)=\frac{1}{(2 K+1)} \sum_{k=n-K}^{n+K} x(k)$     (7)

In the context of Raman spectroscopy, estimating the local variance of a spectrum on a point-by-point basis is pivotal. Large variance values typically indicate the presence of spikes or peaks, which are integral to identify and remove for accurate trend estimation. Following the removal of these spikes, interpolation is employed to reconstruct the peak positions.

Interpolation, in its essence, is the process of making an informed estimation about unknown values within a set of known data points. It can also be described as a model-based method for recovering continuous data from discrete data within a defined range of the abscissa. In the field of digital signal processing, interpolation is often treated as a digital convolution operation, which can be efficiently executed using a digital filtering approach.

Spline interpolation, a form of polynomial interpolation, is particularly advantageous in our analysis. The interpolant in spline interpolation is a piecewise polynomial known as a spline. We opted for spline interpolation over traditional polynomial interpolation due to its lower interpolation error, even when utilizing low-degree polynomials for the spline. Furthermore, spline interpolation effectively circumvents the issue of Runge's phenomenon, a problem typically associated with high-degree polynomial interpolation [13]. The structural formation of a spline is depicted in Figure 5.

Figure 5. Spline shape

To understand the calculations in spline interpolation, let us delve into its formulation. The spline interpolation of a data sequence can be mathematically represented as follows [13]:

$\hat{f}(x)=\sum_{k \in Z} c\left(x_k\right) \beta^n\left(x-x_k\right)$       (8)

where, $\beta(x)$ denotes the interpolation basis function. Here, $x$ and $x_k$ represent continuous and discrete spatial distances, respectively. The terms $c\left(x_k\right)$ are referred to as the interpolation coefficients, which must be estimated before proceeding with the interpolation process. The basis function specific to cubic spline interpolation is given by:

$\beta^3(x)=\left\{\begin{array}{lc}\frac{\frac{2}{3}-|x|^2+\frac{|x|^3}{2}}{\frac{(2-|x|)^3}{6}} & 0 \leq|x|<1 \\ 0 & 1 \leq|x|<2 \\ 0 & 2 \leq|x|\end{array}\right.$    (9)

From Eqs. (8) and (9), we obtain:

$\begin{aligned} & \hat{f}(x)=c\left(x_{k-1}\right) \\ & {\left[(3+s)^3-4(2+s)^3+6(1+s)^3-4 s^3\right] / 6} \\ & +c\left(x_k\right)\left[(2+s)^3-4(1+s)^3+6 s^3\right] / 6 \\ & +c\left(x_{k+1}\right)\left[(1+s)^3-4 s^3\right] / 6+c\left(x_{k+2}\right) s^3 / 6\end{aligned}$      (10)

where, s is a distance represented as follows:

$s=x-x_k$   (11)

During the simulation, we continuously estimate variance point-by-point until a spike is detected. Upon identifying a spike, the estimation process is temporarily halted to facilitate the spike removal. Once the spikes are removed, cubic spline interpolation is applied to the resulting vacant locations. This step allows us to accurately estimate the trend shape and amplitude. Subsequently, the estimated trend is subtracted from the original Raman spectra to obtain the final corrected spectra [14].

4.2 Basics of denoising techniques

Variable pattern noise in signal processing can generally be categorized into two models: additive and multiplicative. In the additive noise model, the noise is inherently additive, a characteristic commonly observed in Raman spectra noise [15]. This paper delves into various noise reduction and removal techniques, including wavelet denoising, moving average filtering, median filtering, and thresholding. These methodologies will be explored in detail in the subsequent sub-sections, each tailored to address the nuances of additive noise in Raman spectra.

4.2.1 Wavelet denoising

The process of noise removal using wavelets begins with the decomposition of the noisy signal into wavelet coefficients, which include both approximation and detail coefficients. After this decomposition, the detail coefficients undergo modification, which is governed by a thresholding function. This requires the selection of an appropriate threshold value. The final step in this process is the reconstruction of the signal by applying the inverse wavelet transform to the modified coefficients [16].

Two types of thresholding functions are commonly employed in this context: hard and soft thresholding [14, 15, 17]. The hard thresholding function retains the input value only if it exceeds the threshold, setting all other values to zero as illustrated in Eq. (12). In contrast, the soft thresholding function is illustrated in Eq. (13).

$f_{\text {hard }}(x)= \begin{cases}x & |x| \geq T H \\ 0 & |x|<T H\end{cases}$   (12)

$f_{\text {soft }}(x)=\left\{\begin{array}{cc}x & |x| \geq T H \\ 2 x-T H & T H / 2 \leq x<T H \\ T H+2 x & -T H<x \leq-T H / 2 \\ 0 & |x|<T H / 2\end{array}\right.$     (13)

where, TH denotes the threshold value, and x represents the coefficients of the high-frequency components of the DWT.  

4.2.2 Moving average filter

The moving average filter serves as an effective smoothing tool for data processing. Its implementation is based on a straightforward characteristic function, which can be mathematically represented as follows [18]:

$H(z)=\left(z+1+z^{-1}\right) / 3$    (14)

4.2.3 Median filter

Median filtering represents a form of nonlinear signal smoothing, primarily aimed at mitigating signal spikes induced by impulsive noise. In the median filtering process, an odd number of signal samples can be collected and organized in a specific order. The median value, which is the middle value after sorting, is then extracted. For a median filter with a length $N=2 K+1$, the output of the filter can be mathematically expressed as follows [19, 20]:

$Y(n)=M E D[X(n-K), \ldots, X(n), \ldots, X(n+K)]$  (15)

where, $X(n)$ and $Y(n)$ represent the $n^{\text {th }}$ sample of the input and output sequences, respectively.

It is important to note that this form of median filtering is non-recursive. This means that the estimate of the median filter output at any given sample time is independent of previous outputs of the median filter. However, there exists another variant of median filtering, known as recursive median filtering. For a recursive median filter with a window length $N=2 K+1$, the output is defined as [17]:

$\begin{aligned} & Y(n)=M E D[Y(n-K), Y(n-K+1) \ldots, Y(n-1), X(n), \ldots, X(n+K)]\end{aligned}$ (16)

This recursive process incorporates a type of feedback mechanism, which is found to be more effective in noise reduction.

4.2.4 Thresholding method

The soft thresholding function is a valuable tool for denoising and can be effectively applied within a transform-domain representation. A critical aspect of this application is the orthogonality of the transform. When the transform is orthogonal, the noise in the transform domain retains the same correlation function as that of the original noise present in the signal domain. Consequently, this means that when the transform is orthogonal, white noise observed in the signal domain is equivalently transformed into white noise in the transform domain. This orthogonality is essential for the effective functioning of the soft thresholding function in noise suppression. The process of noise suppression using the soft thresholding function in the transform domain relies on the maximum function, as detailed in reference [21].

$F_{1 d}=F_1 \times a b s\left(F_1\right) / \max \left(a b s\left(F_1\right)\right)$     (17)

where, F₁ is the noisy signal and F_1d is the corrected signal.

Then, the noise variance is varied from 0 to 1 by a step of 0.1, and the MSE is estimated every time according to the Eq. (18) [22]:

$M S E=\frac{1}{n} \sum_{i=1}^n\left(x-x_c\right)^2$    (18)

where, $x$ is the original signal with length $n$ and $x_c$ is the corrected signal.

4.2.5 FFT filter

The Fourier filtering methods utilize the Fourier transform, a pivotal tool that decomposes a signal into its constituent frequency components. By strategically removing high-frequency components, these methods can effectively achieve denoising. In a continuous-time context, the process can be represented as follows [23]:

$f_{\text {den }}(t)=F^{-1}\left(F\{f(t)\} . X_A(f)\right)$     (19)

In this equation, F symbolizes the Fourier transform, and X_A is the characteristic function of the set A. The parameter X_A represents the cut-off frequency, which is a crucial factor in determining which frequency components are to be filtered out in the denoising process.

This paper presented a study of simultaneous denoising and trend removal, specifically focusing on the Raman spectrum. The choice of the Raman spectrum is due to its characteristics of low Signal-to-Noise Ratio (SNR) and trends with large magnitudes. A key finding from this study is that the wavelet denoising emerges as the most effective denoising technique across all the cases examined. Additionally, it has been observed that applying denoising to a signal prior to trend removal yields better results than conducting the trend removal, initially. Specifically, in the sinusoidal case, trend removal was found to be more straightforward compared to the sigmoidal and linear trend cases. The insights garnered from this study are particularly relevant and applicable in the fields of physics, chemistry, and biology. Given its broad applicability and effectiveness, this research is of considerable value, especially in developing regions such as countries in Africa, where such advanced analytical techniques can greatly enhance scientific and technological capabilities.

•Limitations of this study

Scope of spectrum types: While the study extensively explores the Raman spectrum, the application of the proposed denoising and trend removal techniques may have varying degrees of effectiveness on other types of spectra. Future studies could expand to include a wider range of spectrum types.

Algorithm complexity: The proposed algorithm, especially with cubic spline interpolation, may involve computational complexities that could affect its efficiency in real-time applications. Furthermore, optimization of the algorithm could be explored to enhance its applicability.

Noise variability: The study primarily focuses on scenarios with low SNR and trends with large magnitudes. The performance of the proposed techniques in situations with different noise characteristics, such as higher SNR or different noise types, needs to be thoroughly examined.

Hardware limitations: The effectiveness of the proposed techniques in a real-world settings may be influenced by hardware capabilities.

•Future work

Broadening the spectrum range: Future research could investigate the application of the proposed algorithms to different types of spectroscopic data, such as infrared or ultraviolet spectra, to assess their versatility and adaptability.

Algorithm optimization for real-time analysis: Developing more computationally efficient versions of the proposed algorithms would be valuable, especially for real-time spectral analysis applications.

Exploring automated trend identification: Automating the process of trend identification could enhance the usability of the proposed techniques, making them more accessible to users without extensive technical expertise.

Adaptation to different noise types and levels: Further research could focus on adapting the proposed techniques to handle a wider range of noise types and levels, ensuring broader applicability.

Hardware-software integration: Investigating the integration of the proposed techniques with different types of spectroscopic hardware, particularly in resource-limited settings, would be crucial to enhance their practicality and accessibility.

Cross-disciplinary applications: Exploring the application of these techniques in various scientific fields, such as environmental monitoring or medical diagnostics, could open new directions for research and development.