1. Introduction

Coherent anti-Stokes Raman scattering (CARS) microscopy and spectroscopy have long been proposed as a method for label-free chemical analysis and imaging at orders of magnitude higher speeds than traditional Raman micro/spectroscopy and at higher spatial resolution than afforded by infrared absorption [1–3]. CARS uses coherent stimulation of molecular vibrations that enables both more efficient scattering of the probe beam and non-isotropic radiation of the output signal, which dramatically improves collection efficiency [4]. Additionally, the co-generation of a so-called nonresonant background (NRB) acts as a stable homodyne amplifier, greatly amplifying the signal [5,6] — without which CARS would often have limited advantage over spontaneous Raman [7]. Though beneficial, the NRB distorts spectral peak shapes and contributes to the CARS signal having a quadratic polynomial dependence on analyte concentration, potentially spoiling quantitative analysis capabilities. To combat such challenges, optical methods to reduce the NRB were originally proposed; though, they reduced the overall CARS signal [5,8–10]. Later research revealed that computational methods, based on electromagnetic [11] or information theory [12], could remove the distortion from the NRB without a reduction in signal strength [13]. Additionally, these computational methods result in a signal that is linear with analyte concentration, enabling quantitative analysis of biological, chemical, and polymer samples [14–18]. Ideally, these methods require independent measurements of NRB spectra, which is not currently possible, leading to spectroscopic errors and distortions [19]; though, experimental [16] and computational [19] mitigations have been developed. Even with a known NRB, these methods suffer from artifacts when spectroscopic peaks (and “wings”/“tails”) extend beyond the recorded spectroscopic window. These artifacts distort spectral lineshapes and can even generate new features with similar intensities to the analyte signals. In this work, we present a new method of performing the Hilbert transform for the Kramers-Kronig (KK) relation method of extracting Raman spectral features from CARS spectra. The Learned Discrete Hilbert Transform (LeDHT) uses a synthetically generated set of single-peak training data (with known Hilbert transforms) to learn a matrix representation of the Hilbert transform. After learning (“training” in machine learning parlance), arbitrary new input spectra can be transformed via matrix multiplication, which is computationally fast and efficient. Furthermore, even with synthetic training data, the LeDHT matrix is usable on experimental data with superior performance as will be demonstrated with CARS spectra of glycerol.

2. Theory

2.1 CARS spectroscopy and Raman signal extraction using the Kramers-Kronig relation

CARS is a third-order nonlinear optical mechanism in which the difference frequency between a “pump” photon and a “Stokes” photon is capable of stimulating a Raman (vibrational or rotational) transition from which a “probe” photon can inelastically scatter [2,5,6,20–22]. From a macroscopic and scalar perspective, the generated CARS spectral intensity, $I_{\mathrm {CARS}}(\omega )$, may be described as [19]:

Spontaneous Raman spectroscopy records peaked spectra that are proportional to $\Im \{\chi _\mathrm {r}(\omega )\}$ [23], but CARS spectra are proportional to $|\chi ^{(3)}(\omega )|^{2}$; thus, computational or experimental approaches need be taken to ascertain the phase of $\chi ^{(3)}$. The two primary computational approaches (also known as “phase-retrieval” methods) are the maximum entropy method (MEM) [12] and the Kramers-Kronig relation (KK) [11], which have been shown to be functionally equivalent [13]. Mathematically, using the KK relation as derived in [19], the retrieved nonlinear susceptibility (normalized to the NRB) is:

2.2 Hilbert transform

The Hilbert transform is a widely used linear operator implemented to extract analytic representations of signals from real-valued measurements. Simply put, causal systems may be represented by complex-valued functions in the frequency domain with a defined relationship between the real and complex components [24]. The Hilbert transform lets one calculate the imaginary (real) portion of this function when only the real (imaginary) portion are directly measurable. For example, a molecular system with inhomogeneous broadening may present itself as a Gaussian lineshape through a spectroscopy that measures the real (imaginary) part of the system response, and its imaginary (real) part – with analogous information – would be a Dawson lineshape. There is also a similar relationship for when one only measures the magnitude (phase) under certain conditions [25]. The Hilbert transform is used across a myriad of application spaces such as vibrational analysis [26–28], biomedical signal processing (e.g., cochlear implants [29], neural [30,31] and cardiac monitoring [32]), and spectroscopies (e.g., optical [11,19,33,34], terahertz [35], impedance [36], dielectric [37], magnetic resonance [38,39]), where it is often referred to as a Kramers-Kronig “transform” or “relation”. Further, the Hilbert transform is an integral part of the analysis tool known as the Hilbert-Huang transform: a NASA-developed method most frequently applied to geophysical analyses [40,41]. A significant hurdle to the accurate use of the Hilbert transform, however, is that the discrete implementation may deviate significantly from the continuous form.

The Hilbert transform, $\mathcal {H}$, of a function $f(x)$ on the real line is defined as

Using the Fourier transform is the most common approach in-so-much that it is broadly available in such computational software as GNU Octave, MATLAB, SciPy, and SAS. The aforementioned software packages use a fast-Fourier transform (FFT)-based method developed by Marple [48] that uses symmetry properties of the discrete Fourier transform (DFT) to calculate the DHT. Another Fourier transform approach by Henrici [49], which is implemented in LabView (and in previous work of the author [19]), uses the (continuous) Fourier transformed version of Eq. (6) as the starting point for the computation with a FFT. The mathematical forms and the equivalency of these two methods are derived in Supplement 1 subsection 1.A.

2.3 Errors in the DHT

Significant challenges to accurate DHT computation arise in the case of signals that are non-bandlimited, non-compactly supported, and non-periodic — “general inputs” per [50] — such as is common in spectroscopies. The DHT on signals with these characteristics may have dramatic errors near the window edge — “end-”/“edge-effect” errors; though in fact, the error exists across the entire signal. These errors, as will be demonstrated, may cause significant errors in quantitative analysis of peak heights and relative ratios and may even generate new features.

Alternative DHT algorithms beyond Henrici and Marple primarily focus on computational efficiency and/or accuracy away from the window-edge; though, some methods show improvements for window-centered signals [51,52]. To address window-edge errors, padding schema have been proposed [53] as well as analytical [54] and machine learning methods [55,56] to artificially extend the input signal as to push the edge-effect errors further from the original data. These methods, however, do not stop the generation of the edge effect errors, but rather move it further away from the window center to an artificially larger window. This, in some circumstances may be adequate, but as previously mentioned the edge effect errors permeate the entire window.

The aforementioned errors arise from the discretization and truncation of the input data. For example, the use of the DFT/FFT to effectively perform a circular convolution would see a difference in values from each edge of the signal window as a discontinuity. Similarly, a sudden change in derivative between the edges of the window would act as a form of discontinuity that would ultimately affect the correct calculation (related and bounded by the Paley-Wiener theorem [50]). Two properties of the DHT that are derived below (the second being previously unpublished to our knowledge), provide a means to calculate the total error of using the DHT (supporting derivations are found within Section 1 of Supplement 1):

Next, we show that the variance of a signal is the same as that of the DHT of the signal, which is not necessarily true for the continuous Hilbert transform. First, using Parseval’s theorem:

Though not provided in this paper, a detailed analysis of the DHT through analyzing bounds of decay relations in the DFT/FFT (e.g., see Chapter 6.7 in [50]) can provide insights and bounds onto the shape of the errors for known functions.

2.4 LeDHT

As described above, the DHT has properties that adversely affect its accuracy — especially in the common case of the spectrum extending to or beyond the recording window edge. Additionally, the DHT has fundamental properties that are aberrant when compared to the continuous Hilbert transform. Thus the aim of the LeDHT is to calculate the Hilbert transform on discrete data (spectra) with significantly fewer errors than the DHT. For CARS spectra, this means more accurate phase retrieval and Raman signal extraction than conventional methods. At its core, the LeDHT is a square matrix, $\mathbf {H} \in \mathbb {R}^{N\times N}$, which applies the transform via matrix multiplication: $G=F\mathbf {H}$, where $F$ is a matrix of input spectra and $G$ are the calculated Hilbert transforms. Both $F$ and $G$ $\in \mathbb {R}^{M\times N}$, where $M$ are the number of spectra and $N$ is the length of each spectrum. This matrix approach evolved from the fact that linear discrete transforms may be described by matrix formulations [57] (e.g., Hilbert [58,59] and Fourier [50,60] transforms). These matrix forms, however, are equivalent and generate the same errors and distortions. In the LeDHT, though, the matrix is not prescribed but rather learned through linear regression (least-squares) of synthetic training data with known [continuous] Hilbert transforms. This matrix form may be learned once for a given spectroscopic system and reused indefinitely – assuming one’s choice of lineshape(s) and training width ranges and center ranges are still optimal for the sample and optical systems. The three steps to training and applying the LeDHT are described pictographically in Fig. 1 and will be elaborated on below.

 

Fig. 1. Infographic describing the three steps of implementing the LeDHT: generation of synthetic training data, learning an optimal matrix form of the Hilbert transform, and application to new spectra via matrix multiplication. Note: Scatter plot points sub-sampled (5x) for visual clarity.

Download Full Size | PDF

The first step to implementing the LeDHT is the generation of synthetic training spectra. These $M$ spectra will be used to learn the square LeDHT matrix, $\mathbf {H} \in \mathbb {R}^{N\times N}$ ($M \gg N$), that describes the relationship between a spectrum of length $N$ and its Hilbert transform (See step 1 of Fig. 1). The linear nature of the LeDHT offers significant benefits when generating synthetic data to span the available training space (i.e., possible lineshape parameter values). Of the most consequential benefits is that each training spectrum need only a single peak as more complex spectra can be composed of a linear superposition of single-peak spectra. This greatly reduces the number of possible unique spectra for training when compared to nonlinear models such as neural network-based approaches [61,62] that require full multi-peak spectra. As an example, consider the use of Gaussian lineshapes with corresponding scaled ($2\pi ^{-1/2}$) Dawson functions for the Hilbert transforms. These lineshapes have 2 free parameters: center position and a width parameter. Assuming $n \in \mathbb {Z}^{N}$, center position is an integer value on this line, the width is a positive integer value that also enables at least 3 samples per full-width at half-maximum (FWHM), and the maximum FWHM is $N$, the number of unique peaks is $(N-1)(N/2 -1)\,=\,N^{2}/2 - 3N/2 + 1$. For a 1000-point signal, 498501 unique lineshapes exist under these conditions. Assuming double-point precision, the complete training dataset would be under 4 GB. This is dramatically lower than if full, multi-peak spectra were needed, such as with neural network approaches [61,62], where there are $\approx 2.9\times 10^{85}$ unique spectra for a 1000-point signal length with 1 to 15 possible peaks.

Generation of the synthetic training data has several considerations such as lineshape function and any constraints on training parameter space. From earlier experiments, we found that Dawson lineshapes, of which Gaussian lineshapes are the Hilbert transform, generated more optimal results than dispersive and Lorentzian lineshape pairs. These lineshape types have only two training space dimensions to consider: center position (frequency) and width. A potentially more general lineshape for future inquiry would be using the Faddeeva function in which real part is a Voigt function and the imaginary part is its Hilbert transform. From a practical perspective of this work, the three parameters (two related to width and one to spectral position) would have greatly expanded our training space, limited the fraction of which we could sample (limited by available memory of the computer). The second important consideration are any constraints on the training space parameters. In this work, the minimum peak width (standard deviation for a Gaussian) was 4 spectral pixels, and the maximum was 80 spectral pixels. Additionally, any generated peak must be at least half a FWHM from the window edge; thus, as peaks are centered closer to the window edge, the maximum allowable width is narrower. This gives the training space parameter map a trapezoidal quality as seen in Fig. 1. Additionally, we limited spectral position and width values to integer values (unit of spectral pixels in $n$). A description of how the training data was algorithmically generated is provided within the Methods Section.

The second step in implementing the LeDHT is solving for the $\mathbf {H}$ matrix in $G=F\mathbf {H}$, where $F$ are the input training spectra and $G$ are the known Hilbert transforms. Generally speaking, there are numerous methods to solve for an optimal LeDHT matrix. In this work, we are using ordinary least squares (OLS) which minimizes the residual sum-of-squares: $\mathbf {H} = \mathrm {argmin}_\mathbf {H} ||G - F\mathbf {H}||^{2}$. Beneficially, OLS is the most common linear regression method; thus, there are numerous efficient algorithms for its calculation. We have also tested common “regularized” least squares methods such as ridge and LASSO regression with success though with inferior results for the data within this paper. For different lineshapes and different spectroscopic conditions, these other methods could potentially be superior.

Finally, the LeDHT can be applied to new data via matrix multiplication, which is computationally inexpensive: $G_{exp} = F_{exp}\mathbf {H}$, where $F_{exp}$ is a newly measured input spectrum/spectra and $G_{exp}$ is the LeDHT Hilbert-transformed spectrum to be calculated. Ideally, the solved-for LeDHT matrix for a given spectroscopy system does not need to be re-learned between experiments and can be reused indefinitely, again, as long as the training dataset still optimally reflects the expected future input spectral lineshape characteristics.

3. Methods and materials

3.1 Processing and software

All simulated and experimental data processing was performed on a Dell Precision 7730 laptop with a 6-core i7-8850H processor at 2.6 GHz with 64 GB of RAM. All processing of simulated and experimental data was performed using Python 3.8.5, NumPy, SciPy [63], scikit-learn [64], Matplotlib, and in-house developed Python packages. The traditional DHT used in this paper is the Marple version as available in SciPy. The LeDHT method, several other implementations of DHT methods (and padding schema), and Jupyter Notebook examples have been made available as an open-source Python package, Hilbert, at [65].

3.2 Generation of synthetic data

Synthetic data is generated in a stochastic manner in which the free parameters are randomly selected and accepted based on whether it is a unique combination and is within given bounds. For the simulations, $N$ is 401 with $n$ in the range [-200, 200]. The generated training data is composed of Dawson function lineshapes and their Hilbert transforms: Gaussian lineshapes. The minimum acceptable Gaussian standard deviation is 4 $n$-pixels and the maximum is 80. Furthermore, the center location of Gaussian lineshapes is bounded [-200,200] (see Fig. 1 “Training Space” subimage) but that for a given peak, the center position must be at least half a FWHM from the window edge. Under these conditions, 100000 of these unique single-lineshape spectra are generated. For further stability and improved results, we augmented the data by repeating the dataset 3 times, adding white Gaussian noise (standard deviation of $1\times 10^{-12}$) and adding a constant offset sampled from a normal distribution. Thus the total training data was 300000 single-peak spectra.

An additional 100000 single-peak spectra were generated as test set for validation. This validation was repeated 3 times, for a total of 300000 spectra. These spectra were noiseless and did not have a DC offset as to represent the ideal, expected spectral case. Results with test data that did contain noise and DC offset provided similar results to those presented in the Results Section.

For the simulated multi-peak spectra experiments, the same training data was used as for the single-peak experiments, but the test data (i.e., multi-peak spectra) was generated from linear summations of Dawson function lineshapes (and their corresponding Hilbert transform Gaussian lineshapes). Each randomly generated spectrum was composed of 1 to 15 peaks (uniform sampling). Each peak had an amplitude that was uniformly sampled in the range of 1 to 2 (no units), a width that was uniformly sampled between 4 and 26 (integer valued), and a center position that needed to be at least half a FWHM from the window edge.

For working with the experimental BCARS spectra that are 250 to 460 n-pixels long, we generated 150000 Dawson and Gaussian lineshapes limited to widths of 3 to 50 n-pixels. Next we replicated (3x) and augmented the data in a similar manner as previously described, adding white Gaussian noise (standard deviation $10^{-3}$), and a DC offset sampled from a normal distribution. This provided a total of 450000 training spectra.

3.3 CARS spectrometer

In this work, spectra of neat glycerol (BioUltra, Sigma) were collected on a previously developed broadband CARS (BCARS) microscope that is described in detail elsewhere [6]. The glycerol was mounted between a coverslip (VWR) and glass slide (Ward’s) separated by an adhesive well (GraceBio). Each spectrum was collected with a 3.5 ms dwell time and with average power of the probe and supercontinuum of 20.1 mW and 11.2 mW, respectively. Spectra of the coverslip were collected and used for normalization — a process that effectively suppresses a co-generated, coherent nonresonant background [6,19]. Additionally, dark spectra (250) were collected under the same conditions except with the laser sources not temporally overlapped (but still illuminating the sample). The glycerol spectrum shown in this work is an average of 97 individual spectra.

4. Results

4.1 Simulated

As an initial demonstration, we will compare the DHT, the DHT with extrema-value padding (DHT-Pad) – the most common padding scheme in the CARS community – and the LeDHT. Figure 2(a) is a representative Dawson function spectrum that is spectrally centered, and Figs. 2(b)-(d) show the respective outputs of the DHT, DHT-Pad, and LeDHT. The DHT generates the most severe errors, with new downward-facing features at the spectroscopic window’s edges. The magnitude of these errors are such that the mean of the DHT spectrum is zero as described previously in the Theory section. The DHT-Pad, by visual inspection, agrees with the true Hilbert transform with an apparent slight constant offset. Finally, the LeDHT appears to agree with the true Hilbert transform. Using mean squared error (MSE) as the error metric, the LeDHT is approximately 6 orders of magnitude better than the DHT-Pad and 8 orders of magnitude improvement over the DHT.

 

Fig. 2. (a),(e) Example Dawson function input spectra. (b),(f) DHT results showing significant errors. (c),(g) Results using the DHT with extrema-value padding (DHT-Pad). (d),(h) Results using the LeDHT. (i) Mean MSE ($\pm$ 1 standard deviation) calculated for the 300000 stochastically generated test spectra as a function of peak center position. (j) Histogram of optimal method based on minimum MSE for test spectra dataset. Y-axis shows fraction, parenthetical values are number of spectra.

Download Full Size | PDF

Next, we examine the results of using the same input Dawson spectrum but shifted closer to a window edge, as shown in Fig. 2(e). Similar to before, the DHT produces erroneous downward facing features that are significant in amplitude as shown in Fig. 2(f). Additionally, the DHT peak maximum value is approximately 20 % less than the true value. With the spectral peak near the window edge, now the DHT-Pad is demonstrating significant inaccuracies. The peak maximum is significantly above the ideal case ($>20 \%$ at the peak maximum) and there is a clear distortion in the baseline. Figure 2(h) again demonstrates the superior performance of the LeDHT. The MSE of the LeDHT is approximately 550 times lower than the DHT-Pad, and nearly 1200 times lower than for the DHT.

To further examine the peak position-error relationship, Fig. 2(i) shows the MSE (on a log scale) as a function of peak position from the 300000 test data spectra. The LeDHT has the lowest mean MSE of all methods across all peak positions. The next best method is DHT-Pad, which outperforms DHT across the entire window save the edges. And finally, we calculated for each of the 300000 test spectra, which method had the best performance as shown in Fig. 2(j). The LeDHT had the best performance in 299974 of 300000 spectra ($\approx$ 99.99 %), with the DHT having better MSE 26 times ($\approx$ 0.01 %). All 26 results were from narrow peaks centered at the window edges, for which improved training data with some added emphasis on these cases would likely improve LeDHT performance. The DHT-Pad was never the optimal choice, which indicates that the LeDHT always outperforms the DHT-Pad. Section 2.A of Supplement 1 contains the converse scenario in which input spectra are single Gaussian peaks and the resultant Hilbert transforms (Dawson functions). Like the results presented here, the LeDHT is similarly superior. Section 2.C of Supplement 1 also demonstrates using the Gaussian training set above on test data that is Lorentzian lineshapes. Surprisingly, the LeDHT still outperforms the DHT and DHT-Pad on 300000 test spectra.

Next we analyzed the DHT, DHT-Pad, and LeDHT performance when applied to synthetically generated multi-peak spectra. 30000 spectra were generated stochastically, each with 1 to 15 peaks. Figure 3(a) shows an example single spectrum composed of a linear combination of 13 Dawson function lineshapes. Figures 3(b)-(d) shows the corresponding transforms. Looking at the DHT, again one can see downward sloping features at the window edge that are larger than any real peak (though negative) with over 1000 % relative error at the window edges (over 5800 % at the right side). Additionally, one can see that the entire DHT spectrum is offset below ideal outcome leading to a minimum relative error of $\approx$ 29.5 %. Similarly, the DHT-Pad spectrum is entirely shifted above the ideal spectrum with a minimum relative error of $\approx$ 8.5 %. The DHT-Pad also shows amplitude errors that become worse at the window edges: greater than 300 % relative error at the left edge and 550 % at the right edge. In contrast, the LeDHT is significantly less biased with a minimum relative error of $\approx$ 0.04 %. The window edge extrema relative errors were 13 % (left) and 0.15 % (right). Figure 3(e) shows box plots of MSE generated by the three methods. The LeDHT median MSE is $\approx$ 4.5 orders of magnitude less than that of the DHT, and approximately 3 orders of magnitude less than the DHT-Pad. Figure 3(f) further affirms the superiority of LeDHT, identifying that the LeDHT outperformed, based on MSE, the other methods in all 30000 test spectra. Section 2.B of Supplement 1 contains the converse scenario in which input multi-peak spectra are composed of Gaussian lineshapes and the resultant Hilbert transforms, displaying similarly superior results.

 

Fig. 3. (a) Example synthetic spectrum composed of 13 Dawson functions. (b),(c),(d) The true Hilbert transform and the calculated Hilbert transforms using DHT, the DHT applied to extrema-value padded spectra (DHT-Pad), and the LeDHT, respectively. (e) Box plots of MSE for each method. Note the whiskers displays the extreme range of values. (f) Histogram of the optimal method (based on minimum MSE) for 30000 test spectra. Y-axis shows fraction, parenthetical values are number of spectra.

Download Full Size | PDF

As a final note, we compared the computational time of the DHT, DHT-Pad, and LeDHT methods. We do not assume that any of our algorithms are optimal; thus, there is likely much room for improvement. For 100000 single-peak spectra, the DHT and DHT-Pad took approximately 4 seconds and 15 seconds, respectively over 3 trials. The LeDHT, took 10 seconds to train on 300000 spectra, and only 0.4 seconds to 0.5 seconds to calculate 100000 test spectra. And as mentioned previously, the LeDHT can be trained once and reused in future experiments; thus, the LeDHT was approximately an order of magnitude faster than the DHT when already trained, and over 30 times faster than DHT-Pad.

4.2 Experimental CARS spectra

Finally, we applied the LeDHT to experimental CARS spectra of neat glycerol, varying the proximity of the spectrum to the window edge to demonstrate the superior repeatability and stability of the LeDHT method. Though this scenario may seem artificial, for broadband (>3500 cm$^{-1}$) or low-frequency (<500 cm$^{-1}$) applications, peaks extending to the window edge become probable (especially in chemical and biological samples). The Raman phase ($\phi$, Eq. (5)) and Raman:NRB ratio ($\chi _R/\chi _{NR}$, Eq. (5)) were calculated as described in the Theory and Methods section. These quantities were calculated at different signal bandwidths: from $\approx$ 2080 cm$^{-1}$ down to $\approx$ 1020 cm$^{-1}$ in 23 increments. The training data was composed of 450000 Gaussian lineshapes augmented with noise and DC bias as described in the Methods section.

Figure 4(a) shows the input to the Hilbert transforms, $-0.5 \textrm { ln} (I_{\mathrm {CARS}}/I_{\mathrm {NRB}})$, where the NRB was recorded from coverslip glass. Note the negative sign when compared to Eq. (5), due to our experimental data vectors being oriented with decreasing $\omega$ vectors (i.e., the $\omega$ vector is numerically high-valued to low-valued). Figures 4(a) and (g) show the Raman phase and Raman:NRB ratio, respectively, as extracted using the DHT. As expected, the DHT generates large negative features at the window extrema. Figures 4(b) and (h) show the DHT-Pad results that are by inspection significantly improved. One can see, though, that as the bandwidth of the spectrum is changed, the resulting Hilbert transform is modulated in amplitude, especially at the baseline near 2600 cm$^{-1}$ and at the OH-stretch maximum near 3400 cm$^{-1}$. Conversely, the peak at 2949 cm${-1}$ is stable. A common approach in CARS spectroscopy is to set the baseline to 0 at the window edge. Doing so, in this case, at approximately 2600 cm$^{-1}$ would adversely cause the stable amplitude peak at 2949 cm$^{-1}$ to modulate. Finally, Figs. 4(c) and (i) show the LeDHT results that are significantly more stable with window width and only appear to fluctuate mildly within the OH-stretch region. To qualify the stability of the LeDHT algorithm, Figs. 4(d) and (j) present violin plots of the peak phase and Raman:NRB at 2887 cm$^{-1}$, 2949 cm$^{-1}$, and 3400 cm$^{-1}$ across the various window widths. In all six plots, the vertical extent of the LeDHT violin (indicating variation) is significantly smaller than for DHT and DHT-Pad. For example, for the 3400 cm$^{-1}$ peak, the mean $\pm$ standard deviation phase for the DHT is (0.50 $\pm$ 0.12) rad, (0.8 $\pm$ 0.04) rad for the DHT-Pad, and (0.83 $\pm$ 0.01) rad for the LeDHT. Thus in short, the LeDHT standard deviation is about 12 times less than the DHT and 4 times less than the DHT-Pad. For the Raman:NRB calculations, the LeDHT standard deviation is similarly 10 times less than the DHT and 3 times smaller than the DHT-Pad.

 

Fig. 4. (a) Black line associated with right y-axis: negative half natural log of experimental BCARS spectrum from glycerol normalized to coverslip, which is the input to the Hilbert transform for phase calculations. Green lines: DHT-calculated phases when the input is truncated at different spectral widths. (b) Phase calculations using extrema-valued padded spectra into the DHT (DHT-Pad). (c) Phase calculation using the LeDHT. (d)-(f) Violin plots of the phase magnitude for three prominent peaks. (g)-(i) Raman-to-NRB ratio calculations using DHT, DHT-Pad, and LeDHT methods, respectively. (j)-(l) Violin plots of the extracted Raman-to-NRB ratio for prominent peaks.

Download Full Size | PDF

5. Discussion and conclusion

In this work, we have presented a learned matrix approach for calculating the Hilbert transform with superior performance compared to the DHT and a padded DHT schema. These demonstrations were performed on a variety of synthetic signals as well as experimental CARS spectra. And though the LeDHT demonstrated orders-of-magnitude improvement in MSE on the synthetic datasets and heightened stability with experimental data, there are still numerous opportunities for further gains and new application spaces to consider.

There are many extensions and improvements to LeDHT that could further enhance and expand the method. For example, the training data could be augmented to explicitly diminish noise. In this case, it may be advantageous to repeat the training dataset with random additions of noise solely on the input. Our initial attempts at this were successful though areas of low-signal and high-noise were affected in a similar manner to traditional [convolution] filtering methods (e.g., bandpass, Savitzky–Golay) with erroneous lineshape-like features being generated within the noise. Traditionally, we have had more success with separate processes of statistical denoising [19].

Another opportunity for improvement is methods to effectively train with smaller datasets, which in turn enables the generation of larger and larger training datasets with unique features. For examples, we have experimented with generating large training datasets and factorizing with singular value decomposition (SVD), which effectively generates a much smaller training dataset with each entry being more complex than a single lineshape. The 300000 training spectra used for Figs. 2 and 3, for example, could be reduced down to 401 (matching the size of $N$) – a nearly 750x reduction. In the long term, this could enable a training dataset to be originally generated from numerous different lineshapes, potentially enabling a more universal LeDHT matrix that could be shared.

We are also exploring alternative regression/optimization methods for finding the LeDHT matrix. This matrix may be more optimal in terms of MSE, or it may be more optimal in the sense of having a symmetry (e.g., a Procrustes problem) that can be more easily human-interpreted or implemented in hardware. Thus far we have worked on the former problem: implementing Ridge and LASSO regression methods; though linear regression has performed most optimally.

And finally, we are looking to apply LeDHT to other signal types and spectroscopic modalities. Although CARS is a relatively niche technique, methods such as NMR [38,39], dielectric [37], and terahertz [35] spectroscopies have also demonstrated use of the DHT.

In closing, we have developed LeDHT and demonstrated significantly improved MSE over the traditional DHT and with signal extrema padding. Our method is fast, stable, and improves opportunities for quantitative analysis of spectroscopic signals. Additionally, an open-source Python library for this method is freely available to enable community development and foster adoption.