SNR enhancement in brillouin microspectroscopy using spectrum reconstruction.

Brillouin spectroscopy can suffer from low signal-to-noise ratios (SNRs). Such low SNRs can render common data analysis protocols unreliable, especially for SNRs below ∼10. In this work we exploit two denoising algorithms, namely maximum entropy reconstruction (MER) and wavelet analysis (WA), to improve the accuracy and precision in determination of Brillouin shifts and linewidth. Algorithm performance is quantified using Monte-Carlo simulations and benchmarked against the Cramér-Rao lower bound. Superior estimation results are demonstrated even at low SNRs (≥ 1). Denoising is furthermore applied to experimental Brillouin spectra of distilled water at room temperature, allowing the speed of sound in water to be extracted. Experimental and theoretical values were found to be consistent to within ±1% at unity SNR. Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

Introduction

Brillouin spectroscopy is a near century old technique relying on the scattering of incident photons from thermally excited acoustic fluctuations in a medium (i.e. phonons) [1]. By virtue of energy conservation, the scattered photons can possess a different frequency to those incident if a phonon is either created or annihilated as part of the interaction. Measurements of this frequency shift, requiring a spectrometer with sufficiently high spectral resolution and throughput, can elicit useful information about phonon velocities and have therefore been extensively employed to study the elastic properties of materials [2,3]. Moreover, study of the linewidth of the Brillouin spectral peaks enables measurement of phonon lifetimes thereby also enabling viscous characteristics to be quantified [4,5]. More recently, Brillouin spectroscopy has seen a renaissance as it has developed from a spectroscopic technique into a more powerful hyper-spectral imaging modality in which mechanical information is mapped with micron level resolution [6-10]. This is particularly attractive for its applicability in, for example, in vivo diagnostics and cellular imaging (see e.g. [11] for a recent review of biomedical applications).

One drawback of using spontaneous Brillouin scattering, however, is the intrinsically low signal-to-noise ratio (SNR) which in turn leads to relatively long image acquisition times. This issue is further exacerbated in fibre-based systems and remains one of the main challenges in designing a clinical Brillouin endoscope device [12]. The advent of improved instrumentation, for example, tandem Fabry-Pérot interferometers, has aimed to address this challenge. Although a spectral contrast of up to can be achieved [13], the acquisition time is nevertheless on the order of seconds per spectrum [14] and is hence less suitable for in vivo applications. Whilst Brillouin measurements can in principle still be made in spite of weak signal strengths, it has been shown that the resulting accuracy and achievable precision of extracted frequency shifts or linewidths are dependent on the SNR [15,16]. Stimulated Brillouin scattering offers an alternative means to overcome this limitation, and has been used for example in combination with heterodyne detection to efficiently measure large tissue ( mm) samples [17]. Although a promising addition to the Brillouin toolbox allowing large volume measurements to be made in a practical amount of time, the use of spontaneous Brillouin scattering remains the more attractive option for in vivo applications because it requires lower incident power, which is in turn more cost-efficient and helps avoid photo-damage to the sample. Consequently, mitigation of low SNR in spontaneous Brillouin spectroscopy is highly desirable.

Numerous techniques to increase the SNR of a Brillouin system optically can be found in the literature, using for example adaptive optics [18], heterodyning [19] or interferometric filtering [20]. Here, however, we present a complementary software based approach, which can be used either independently or in conjunction with existing experimental techniques. By adopting a purely signal processing approach we therefore enable Brillouin measurements to be made in systems that are typically more challenging to optimise optically e.g. an endoscope. There exist numerous algorithms for signal reconstruction from noisy data. Here, we restrict our attention to the maximum entropy reconstruction (MER) and wavelet analysis (WA) techniques, due to their applicability in a wide scope of experimental situations, particularly in experiments with low SNR and little knowledge of the sample. These methods are introduced in Sections 2.1 and 2.2 respectively. The fundamental precision limit given by the Cramér-Rao lower bound is introduced in Section 3.1, which is then subsequently used to benchmark the performance of the proposed numerical schemes as a function of SNR in Section 3.2. Performance improvements relative to more conventional, pure Lorentzian fitting strategies are also presented. Finally, we apply the reconstruction algorithms to experimental Brillouin spectra of distilled water in Section 4. Note that throughout this work we define the SNR as the ratio of the spectral intensity to the standard deviation of additive noise sources, such as detector noise.

Numerical methods

Inverse problems are commonly encountered in optics and often require the fitting of empirical models to extract useful information from otherwise incomplete and noisy data [21]. Within the context of Brillouin spectroscopy, the frequency shifts (or linewidths) of the inelastic spectral peaks are commonly extracted through least-square fitting of multiple peaks (typically Lorentzians) to noisy intensity data taken with a spectrometer. This is similar to the well-known estimation problem encountered in localisation microscopy, for which the position of individual fluorophores is estimated to extract spatial information with super-resolution [22-24]. Existing estimation techniques can be loosely categorised as fitting or non-fitting, depending on the a priori information that is available to the algorithm [25]. The family of algorithms that fit to a known model, such as least-square fitting, usually produce more reliable results at the cost of requiring that some initial parameters are known beforehand. There are also non-fitting methods, such as centroid [26] and learning-based [27,28] techniques, which require no or minimal information to be known a priori, making them potentially faster and more robust, albeit they are frequently not applicable in all noise regimes [29]. Although these are successful strategies in their own right, the resulting solutions are susceptible to noise and especially in the low SNR regime can suffer considerably from bias and poor precision. Noise reduction algorithms can, however, be applied to the collected data before processing so as to enable improved parameter extraction. Linear or nonlinear smoothing filters, for example, are standardised pre-processing procedures that are applied to Raman and infra-red spectral data [30,31] with reasonably good SNR. Although computationally more demanding, more complex algorithms, such as statistical estimators [32,33] and principal component analysis [34], can often reduce noise whilst preserving spectral features to a greater degree. The effectiveness of the techniques described so far either require higher SNRs or certain prior knowledge and assumptions. In this work we therefore consider reconstruction techniques which require minimal prior knowledge and are able to successfully estimate spectral parameters even at low SNRs, namely maximum entropy reconstruction and wavelet analysis. The former is noted for being effective even at low SNRs, whereas the latter is a ‘tried-and-true’ signal denoising method that is especially effective for data degraded by broadband noise. Each method is presented below and may be used independently or in tandem to construct a versatile denoising scheme that can handle a variety of datasets from Brillouin experiments.

Maximum entropy reconstruction (MER)

Experimental data in spectroscopic applications typically comprises of a series of discrete intensity values, which we here denote by . When spatially or angularly dispersive spectrometers are used, each datum corresponds to the intensity recorded on each of the individual pixels of a segmented detector. For scanning etalon type spectrometers, the data correspond to intensities recorded sequentially in time for each etalon configuration. Observed spectra are however corrupted by noise and thus the determination of the ground truth spectrum is an ill-posed inverse problem. A criterion is therefore required to select a unique spectrum from all possible solutions. The MER method provides one such selection criterion [35]. Specifically, the MER technique selects the spectral distribution (described by , ) which maximises the associated entropy subject to the experimental data, i.e. where is the system response function. For dispersive spectrometers, the entropy can be physically interpreted as the number of bits of information needed to encode on which pixel a single photon was detected [36].

Practically speaking, in the presence of noise Eq. (2) is too restrictive and a solution to the maximisation problem does not exist in general due to the noise degradation of the data. As such it is more appropriate to relax the constraint to allow for an error between the reconstruction and the experimental data. In the presence of additive Gaussian noise for example, the normalised mean square error, is distributed [37,38], where is the noise variance on pixel . Maximising the entropy, , subject to the weaker constraint , where is a constant, can therefore allow a solution to be found. Through variation of , the confidence level of the reconstruction can be varied, with a confidence interval corresponding to [39]. Imposition of the constraint on can be achieved using the method of Lagrange multipliers [40] whereby we construct the Lagrangian function where is a positive Lagrange multiplier. The Lagrangian, , can be maximised with respect to for any fixed value of . Formally the maximum entropy solution corresponds to the value of for which . In practice, however, we have found that similar solutions are found over a broad range of such that to improve algorithm speed we use a fixed value of and it is only changed as an initial parameter as the SNR of the data varies. Maximisation of can then be performed by iteratively updating until the extremum is reached. A simple update strategy is that based on gradient ascent [41] whereby where denotes the iteration number and here is the gradient of with respect to . When is chosen appropriately convergence of to the value which maximises is guaranteed, thus solving the problem: Convergence of this method can, however, be improved by modifying the search direction using the conjugate gradient technique [42]. In this method, instead of using to update , a direction that is conjugate to a set of previous search directions is used. We note, for each iteration , it is possible to find a value of that allows to be locally maximised within the subspace spanned by the previous search directions using the Barzilai-Borwein method [43]. In this work however we instead use a line search approach with two such ‘conjugates’ (i.e. ) based on the Wolfe conditions [44] since this approach is less computationally demanding. Specifically, for each iteration our algorithm updates the value of for varying until the Wolfe conditions are met. The process is iterated over until a termination condition is met, namely that falls below a predefined threshold of . This is a direct measure of the balance between the entropy and constraint during the optimisation, and the true maximum entropy solution should yield a value of 0 [36]. Using this algorithm, which we implemented in Matlab, reconstruction of a single Brillouin spectrum takes ms on a standard office workstation, albeit run times have not been optimised. The speed and stability of the algorithm can be further improved for example by using more advanced search schemes [36], differentiation algorithms [45] or parallelisation [46].

Although the MER algorithm can give good results, it is important to note that its ability to produce meaningful spectra is not without bound. The requirement that the reconstructed spectrum is consistent with the data equates to placing an upper limit on the error i.e. as discussed above. If, however, the unconstrained maximum entropy solution for satisfies this inequality, no meaningful reconstruction is possible since this implies the data is too noisy. This is directly related to the fact that contours of in the vector space defined by f are convex ellipsoids. Moreover, the entropy surfaces are also convex, implying that if a maximum entropy solution exists it lies on the boundary where [47].

Unlike the majority of localisation methods in microscopy, we note that one advantage of the MER is that it makes no assumptions about the data and in theory requires no a priori information. Indeed the algorithm intrinsically assumes we are maximally ignorant about the underlying spectrum [35], whereas localisation algorithms usually function optimally with a well-defined set of assumptions, for example that molecules do not spatially overlap. Nevertheless, it is possible to include a priori information in our reconstruction, if for example the position or width of the likely peaks are approximately known. In turn, such additional information can lead to improved speeds and a higher success rate. Incorporation of such knowledge requires introduction of additional terms into the Lagrangian function. A careful balance must, however, be struck between consistency of the reconstruction with the actual data and the bias towards the a priori information introduced.

Reconstruction by wavelet analysis (WA)

Many conventional noise filtering techniques, such as Fourier filtering or smoothing filters, rely on the fact that the signal of interest and the noise are spectrally diverse. Consequently they can struggle to give good results in the presence of broadband noise. Wavelet analysis (WA) based techniques, however, exploit amplitude filtering in the wavelet domain and thus do not suffer in such scenarios [48]. Even with the best detectors, the performance of Brillouin spectroscopy is typically limited by broadband noise [16]. In particular, obtained spectra suffer from intensity dependent shot-noise, which for a large enough average signal level can be approximated as following a Gaussian distribution with a uniform power spectrum [39].

To illustrate the WA technique consider an additive noise model whereby the observed spectral data can be written in the form where , and are defined as above and denotes a random noise contribution to the th datum (). This experimental spectrum is first decomposed into a chosen basis of wavelets using the discrete wavelet transform (DWT) [49], such that where we have defined , and is the vector formed from the amplitude coefficients of each individual wavelet (known as a level) in the observed data (similarly and contain wavelet coefficients of each level for the observed spectrum and the noise). Under the assumption that the noise is smaller in amplitude than the spectrum of interest, the noise contribution can be reduced by applying a nonlinear thresholding (commonly referred to as shrinking) of the coefficients of the DWT, each according to a threshold . Nonlinear shrinkage of the coefficients is such that smaller amplitude coefficients (i.e. those describing the noise) are reduced to a greater extent than those of the spectrum. Thresholding can be performed using either hard or soft limits, although soft limits are usually preferable in spectroscopy in order to preserve spectral features. By limiting the consideration to only white noise, the thresholding strategy is further simplified as it can be proven statistically that a universal threshold (i.e. the same for each level) can be set [50]. Finally, an inverse DWT is performed on the ‘denoised’ wavelets, yielding a reconstructed spectrum [51]. Although WA is a quick and easy means of extracting spectral features from noisy data, at lower SNRs the argument for universal thresholding is somewhat questionable and the process is then seen to fundamentally rely on some a priori knowledge of the shape and location of anticipated features thus rendering the process inevitably subjective. In these cases, there is also a limitation which comes from the assumption that the signal is strong enough compared to the noise level which allows the separation in amplitude, meaning that useful denoising may not always be possible.

In this work we employ the Matlab Wavelet Analyzer application and use an adaptive thresholding algorithm to denoise spectra based on the SNR which allows reconstruction of a single spectrum in less than 10 ms. The threshold is generally set to be level-independent and follows the relationship given in Ref. [50]: where is the experimentally determined noise level, which is constant in the case of white noise. In principle, even in the presence of intensity-dependent noise i.e. shot noise, this can be set to adaptively vary for different levels. For example, it is possible to incorporate level-dependent thresholding schemes where for each level can be estimated statistically [52].

Assessing algorithm performance

Informational limit of Brillouin systems

To assess the performance of denoising algorithms in Brillouin spectroscopy it is necessary to define suitable performance metrics and to have a benchmark to which they can be compared. Given the objective of denoising algorithms, a natural choice is to consider the average SNR in the reconstructed spectra as compared to the unprocessed data. Ultimately, however, Brillouin spectroscopy aims to quantitatively study the scattering induced frequency shifts. So as to better reflect the aims of Brillouin spectroscopy, we hence choose to consider the accuracy and precision to which the frequency shift can be determined. From a statistical perspective, the accuracy of an estimate of can be parameterised using the bias, , which is defined as the difference between the true value of as compared to the average frequency shift returned from an estimation protocol i.e. . Note here we use the standard statistical notation whereby denotes an estimate of and denotes an ensemble average. The precision can similarly be quantified using the estimator variance . Although extensive research quantifying precision limits in localisation microscopy exists [53,54] it is worth noting that in the context of inelastic spectroscopy at least two peaks are always needed to estimate the relative frequency shift in contrast to the single peak used in localisation microscopy. Consequently the estimation problem is more akin to determining the separation of two molecules [55,56]. Derivation of the relevant precision limits in the context of Brillouin spectroscopy is given explicitly in Ref. [16].

Fundamentally there is no limit as to how small the bias of an estimator can be, however, the limit of precision to which can be determined follows from the well known Cramér-Rao lower bound (CRLB). Specifically, the CRLB bounds the variance, , of any estimate of the frequency shift according to , where is the associated (noise dependent) Fisher information [35]. Assuming a dispersion limited spectrometer with a pixelated detector and restricting to Gaussian white noise of variance (as used in our simulations below), the obtainable precision is bounded according to [16] Here is the linear scale factor between the spatial position on the detector and optical frequency (i.e. ), is the intensity of the Stokes and anti-Stokes peaks (assumed the same) relative to the Rayleigh intensity, is the size of the pixels, is the width of the detector, is the FWHM of the spectrometer response function, is the FWHM of the Brillouin peaks and the average signal to noise ratio on each pixel is , where is the integrated spectral intensity and is the noise variance.

Reconstruction of simulated data

With a suitable benchmark in hand, it is now possible to quantitatively assess the performance of the proposed denoising algorithms, by means of Monte Carlo simulations. The ground truth Brillouin spectrum used assumed the Rayleigh (elastic) peak of frequency THz ( nm) was centred on a detector of total width mm and pixel size , where pixels were used to correspond to an effective bandwidth of GHz. Note that these values were chosen to resemble a realistic experimental set-up. In reality, such coarse pixelation of the spectrum may have an adverse effect on the estimation precision and bias of the relative shift. In particular, pixelation effects become significant when recorded spectral peaks have widths on the order of the pixel size as discussed further in Ref. [16]. To enable a fair comparison between different algorithms we use a fixed pixel size for all simulations. The spectral width of both the Rayleigh peak and Brillouin peaks was assumed to be 1 GHz. The Brillouin peaks were taken to have a relative shift of GHz, whilst the peak intensity of the (anti)-Stokes and the Rayleigh peak was assumed to be and respectively. For simplicity an ideal system response function was also assumed (i.e. where is the Kronecker delta, whereby ). A total of 5000 independent realisations of white Gaussian noise (of varying ) were added to the ‘true’ spectrum. The spectral parameters (most importantly ) of each noisy spectrum were then estimated using Lorentzian fitting before and after the denoising algorithms were applied. The bias and standard deviation of the estimated Brillouin shift were then calculated and are shown in Figs. 1 and 2 respectively as a function of SNR. Note that due to poor performance at lower SNRs we found it necessary to incorporate prior information of the number and positions of spectral peaks into the MER algorithm i.e. the approximate position and width of the Lorentzian function, as discussed further in Ref. [57].

Fig. 1.

Bias of estimates of the Brillouin shift found by Lorentzian fitting of simulated noisy spectral data subject to different denoising algorithms, as calculated using 5000 realisations of simulated noise.

Fig. 2.

Logarithm of the standard deviation of estimates of the Brillouin shift, found by Lorentzian fitting of simulated noisy spectral data subject to different denoising algorithms, as calculated using 5000 realisations of simulated noise.

Inspection of Figs. 1 and 2 demonstrates that both denoising algorithms generally facilitate more accurate and precise estimation of the Brillouin shift. Furthermore, we note that the MER produces superior performance to the WA based method, especially at low SNRs. Although performance of all algorithms are found to degrade at lower SNRs, even for an SNR of unity, the MER allows to be determined with a bias of approximately 1% and an uncertainty of %. It should however be noted that spectrum reconstruction was not possible for all realisations in this case as circa 20% of the generated spectra did not possess ME solutions and were regenerated until the total number of realisations were fulfilled. In contrast, the unconstrained Lorentzian fitting algorithm without denoising was unable to retrieve any meaningful information of the spectra at the same SNR. Both the bias and the uncertainty were well outside of reasonable bounds for SNR , thus they were not plotted for comparison. Fitting performance of unprocessed spectra at an SNR of 10 is still notably degraded by noise yielding a bias (standard deviation) of % (%) in estimates of the Brillouin shift. Similar trends were also found (data not shown) when estimating the linewidths of the obtained Brillouin peaks. In this case, MER also yielded the smallest bias relative to the ground truth width (1 GHz) and the best precision, whilst WA again outperforms Lorentzian fitting. As an indicative example, for an , pure Lorentzian fitting performed poorly giving an average linewidth of GHz, whilst both the MER and WA were noticeably superior, producing linewidths of and GHz respectively.

With respect to the WA method, we note that whilst its utility for denoising is evident from Figs. 1 and 2 for SNRs , at lower SNRs it is unable to function reliably. This degradation in performance is, however, unsurprising because in this high noise regime the strength of the noise is comparable to the signal in contradiction to the assumptions of the WA denoising algorithm. As such, the shape of the reconstructed spectrum becomes very sensitive to the precise manner in which the coefficients threshold is chosen and thus erroneous shifts and distortions of the individual peaks can result. This is particularly detrimental to accurate retrieval of the Brillouin linewidth. Wavelet shrinkage can hence degrade potentially useful information from the signal. Finally we note that although both considered denoising algorithms produce better estimates of the Brillouin shift and linewidth, there remains further scope for improvement since neither achieve the CRLB, although the difference between the theoretical lower bound and the simulation results may be partially accounted for by pixelation effects mentioned above. Statistical methods such as the method of Maximum Likelihood Estimation (MLE) [58] can routinely achieve the CRLB but have yet to be applied to Brillouin spectroscopy. This remains an interesting area for future study.

Denoising of experimental spectra

To demonstrate the applicability of the denoising algorithms to real spectral data, Brillouin spectra of pure water at room temperature were taken using a custom-built Brillouin set-up (see [59] for full details). Variation of the SNR was realised through use of different exposure times. A weak incident power of mW was also used to operate in the low SNR regime. Typical examples of the raw and denoised spectra for an SNR of are shown in Fig. 3 . Critically, the central Rayleigh peak seen in Fig. 3 is saturated which is a feature commonly encountered in Brillouin spectroscopy. This meant that although MER was sometimes possible, even at SNRs as low as 1, the resulting reconstruction did not generally lie within a confidence error due to the large deviation in peak shape that saturation produces. To help mitigate the issues associated with the saturated Rayleigh peak, interferometric suppression of the central Rayleigh peak [60] was used, which in turn improved the contrast of the Brillouin peaks, algorithm convergence rates and reliability of the fitting results. Importantly, Rayleigh suppression can be fully accounted for in MER through the form of the a priori information introduced. If the exact form of the suppressed Rayleigh peak is not known a priori, the Rayleigh peak can alternatively be excluded from the MER entirely, as is the case for the data shown in Fig. 3.

Fig. 3.

Example of a typical unprocessed experimental Brillouin spectrum of distilled water obtained using a 100 ms acquisition time () (blue). Corresponding reconstructed spectra as found using the WA (orange) and MER (yellow) algorithms are also shown. Note that spectra have been vertically shifted to improve visibility.

Upon calibrating the spectrometer and fitting of the (denoised) experimental spectra, the Brillouin frequency shift was found and the corresponding speed of sound determined [61]. Table 1 compares the values obtained using each denoising procedure for different SNRs (or equivalently exposure times). The respective RMS fitting error in percentage is also listed to indicate the relative fitting accuracy. Note that poorer SNRs required stronger interferometric suppresion of the Rayleigh peak. Consequently, algorithm performance for different SNRs is not directly comparable, however, similar relative improvements between algorithms for a given SNR are evident.

Table 1.

Speed of sound in distilled water as obtained from Lorentzian fitting of experimental Brillouin spectra subject to the MER and WA algorithms as compared to no pre-processing.

	SNR ∼ 5(100 ms)	SNR ∼ 2(50 ms)	SNR ∼ 1(30 ms)
MER	1485.6(2.4%)	1490.9(2%)	1493.2(3%)
WA	1481.0(4%)	1450.2(5%)	1508.0(13%)
None	1494.5(34%)	1488.8(27%)	1612.7(37%)

Whilst the speed of sound is sensitive to a range of parameters, it is expected that the speed of sound in distilled water at the laboratory temperature of is around m/s [62]. The MER denoising algorithm was thus once again found to perform best out of the methods considered here, producing experimental values of the speed of sound with low uncertainties and in good agreement with the expected values. Whilst WA also gave good agreement at higher SNRs, a significant bias is seen to occur at low SNRs in agreement with the findings of Fig. 1. Finally, direct Lorentzian fitting was found to produce large errors for all SNRs considered due to the poor quality of the spectrum, such that no meaningful conclusions can be made.

Conclusion

In this work, we have proposed use of two denoising algorithms to improve the accuracy and precision of data analysis in Brillouin spectroscopy. As far as we are aware, this is the first time that such reconstruction schemes have been used to enhance Brillouin spectroscopic data. Specifically, we have discussed both a MER and WA based data processing chain. Through Monte Carlo simulations we have demonstrated that both algorithms enable reduced bias and greater precision when estimating both the frequency shift induced by Brillouin scattering and the associated phonon lifetime. For example, at least a three-fold precision gain was demonstrated when estimating the Brillouin shift when WA was used, whereas up-to an order of magnitude improvement was demonstrated using MER. Indeed, MER was generally found to outperform WA, however, both denoising algorithms enabled superior performance down to SNRs of , a regime in which conventional line fitting fails. Although our simulations assumed white Gaussian noise, as is common place in many experimental scenarios, more complex noise regimes such as coloured and structured noise can also be accounted for through introduction of additional a priori information, use of a pre-whitening processing step or a more elaborate thresholding strategy.

Experimental Brillouin spectra of distilled water were also used to verify the denoising algorithms in a real-world scenario. In particular, improved data analysis was demonstrated with the experimentally determined speed of sound agreeing with theoretical values to within . In conclusion, we have thus shown that a quantifiable SNR enhancement of noisy Brillouin data is easily achievable through numerical methods, thereby unlocking the potential for faster and more reliable measurements. Critically, by adopting a software based approach such improvements are complementary to any technical experimental gains and thus represent an attractive option for Brillouin applications where optical design is impractical, for example in the case of Brillouin endoscopy, or as a means to improve the acquisition speed by enabling operation at lower SNRs. Further development of these denoising schemes can lead to yet further gains, achieving the CRLB, and thus promise to become a useful tool in Brillouin spectroscopy and an embedded part of the data analysis routine in this field.