Denoising applied to spectroscopies – part I: concept and limits

Abstract Some spectroscopies are intrinsically poorly sensitive, such as nuclear magnetic resonance (NMR) and Raman spectroscopy. This drawback can be overcome by using singular value decomposition (SVD) and low-rank approximation to denoise spectra and consequently increase sensitivity. However, SVD limits have not been deeply investigated until now. We applied SVD to NMR and Raman spectra and showed that best results were obtained with a square data set in time domain. Automatic thresholding was applied using Malinowski’s indicators. 6 × 7380 noisy spectra with 41 signal-to-noise ratios were compared to their non-noisy counterparts, highlighting that SVD induces a systematic error for Gaussian peaks but faithfully reproduces shape of Lorentzian peaks, thus allowing quantification. Used carefully, SVD can decrease experimental time by a factor of 2.3 for spectroscopies. This study may help scientists to apply SVD to denoise spectra in a more efficient way, without falling into pitfalls.


A. Introduction
Spectroscopic techniques are of outmost importance in the field of materials analysis.
Among them, Nuclear Magnetic Resonance (NMR) (1) and Raman (2) spectroscopies are very powerful local probes of the chemical structure. Especially, they allow to analyse liquid-state, solid-state, or even gas-state samples. While NMR informs about chemical and magnetic environment of atomic nuclei, Raman spectroscopy provides vibrational and rotational information on chemical entities. Unfortunately, these techniques suffer from a major drawback, namely their intrinsic low sensitivity.
Although focused on NMR and Raman approaches, this work can be easily extended to other spectroscopies.
In the case of NMR, only one nucleus over 10 5 is detected under usual conditions (3). This is due to the low population difference between nuclear spin energy levels, which results from Boltzmann equilibrium. Many factors influence NMR sensitivity: magnetic field strength, sample volume, sample temperature, electronics temperature, radio-frequency coil quality factor and coil filling factor (4). When studying a solid-state sample, situation gets worse due to spectral line broadening, which results either from environment distribution or from relaxation (5). Indeed, NMR relevant anisotropic interactions, as chemical shift anisotropy, dipolar coupling and quadrupolar coupling, are no longer averaged to zero by fast and isotropic molecular motions, leading to spectra spreading over hundreds of ppm or a few megahertz (6).
During the last decades, numerous technical progress has allowed to increase sensitivity of solid-state NMR: Magic Angle Spinning (MAS) up to 110 kHz (7-9), Cross Polarisation (CP) (10), high performance heteronuclear decoupling (11) and NMR magnetic field strength increase (12). Additionally, new very sensitive techniques were 3 developed, for instance micro-coils and microresonators adapted to MAS (13) or Dynamic Nuclear Polarisation (DNP) (14)(15)(16), with which a gain of up to 320 per unit of time was achieved (17).
In the case of Raman spectroscopy, only one photon over 10 6 is detected (18), due to its low scattering cross-section of ≈ 10 -30 cm 2 per molecule, to be compared with ≈ 10 -20 cm 2 for infrared absorption spectroscopy and ≈ 10 -16 cm 2 for fluorescence spectroscopy (19). Furthermore, some samples are fluorescent, thus hiding their Raman spectrum, or can be locally damaged by the laser light during analysis (20). A major advancement in Raman analysis has been achieved with the discovery of Surface-Enhanced Raman Scattering (SERS) (21,22) with an enhancement factor of 10 4 -10 6 (23). However, this effect applies only to transition metals (24) and strongly depends on surface roughness (25). To circumvent these limitations, Raman equipment has been hyphenated with an atomic force microscope or a tunnel effect microscope, leading to Tip-Enhanced Raman Spectroscopy (TERS) (26,27). Another important progress has been achieved through the use of non-linear light sources and picosecond lasers (28,29). Thanks to such improvements, it is now possible to study artworks (30), to map a surface (31) and even to follow cure kinetics of an epoxy resin (32).
In parallel with these instrumental and methodological developments, mathematical and computer tools have become increasingly widespread in the field of data processing (33). In addition to Fourier transform revolution (34,35), other treatments have emerged in NMR spectroscopy and Magnetic Resonance Imaging (MRI). As examples, one may cite Hadamard transform (36), compressed sensing (37), non-uniform sampling (38), or quantitative signal reconstruction from multiple echoes (39). In Raman spectroscopy, pre-processing became of paramount importance to obtain 4 quantitative measurements (40): it is now recommended to suppress fluorescence background (41), to correct cosmic ray spikes (42) and to normalise spectra (43).
Furthermore, a very important mathematical tool family concerns noise reduction (44). Indeed, the above mentioned sensitivity-enhanced spectroscopic equipments are costly and are not accessible to all laboratories, hence the need to denoise. Moreover, even with hardware and methodological improvements, spectra can still be noisy, especially when studying amorphous materials (45,46), for which distribution of bond lengths and bond angles broadens signals and reduce sensitivity. In order to decrease experimental time or in case of unstable samples, signal processing is mandatory to get a reasonable Signal-to-Noise Ratio (SNR). The easiest way to perform noise reduction is smoothing. In NMR, apodisation, especially exponential multiplication, is used prior to Fourier transform (47). In Raman spectroscopy, a polynomial algorithm named Savitzky-Golay is preferred (48). Other options to reduce noise are maximum entropy (49), covariance matrix (50), Wiener's estimation (51), wavelet transform (52), uncoiled random QR denoising (53), and the method initially proposed by Tufts et al. (54) and generalised by Cadzow (55). Singular Value Decomposition (SVD) is an important part of Cadzow's denoising algorithm and its related low-rank approximation (56). An history of SVD can be found in (57). It has been discovered independently by Beltrami in 1873 (58), by Jordan in 1874 (59) and rediscovered by Lanczos in 1958 (60). It is preferred over eingenvalue decomposition, which is less precise, as demonstrated by Laüchli (61).
Recently, Man et al. developed a new SVD application for NMR (81). It was programmed under Java and two versions are currently available: one for processors (82) and the other one for Nvidia graphic cards (83) using CUDA (84). Indeed, graphic cards allow very efficient parallel computations. However, the limits of this approach are not clear: (i) which matrix shape should be preferred? (ii) what is the minimal experimental SNR? (iii) Are denoised spectra quantitative? (iv) is it suitable for other spectroscopies?
Following a previous communication (85), we tried to address these questions.
This work is divided into two parts. In this first contribution (I), we focus on SVD concept and limits. Experimental details will be provided in Section B. Theoretical background on SVD and low-rank approximation concepts is developed in Section C.1.
Hankel and Toeplitz matrices are explored in Section C.2. SNR definitions are given in Section C.3. Section D.1 is devoted to experimental results by applying SVD to solidstate NMR and Raman spectroscopies. The influence of both matrix shape and thresholding are studied in Sections D.2 and D.3, respectively. Time and frequency denoising are compared in Section D.4. The minimum SNR needed to have accurate results is investigated in Section D.5. The impact of SVD on peak shape is considered in Section D.6. Finally, denoising on a real NMR spectrum is analysed in Section D.7.
In a second part (II) (86) we will benchmark SVD using Java, Matlab and Python, on various processors and nvidia graphic cards ranging over 10 and 6 years, 6 respectively. We will try to optimise algorithms, software libraries and hardware capabilities to achieve the fastest possible denoising computation.
This soft chemistry synthetic approach is a suitable route to design hybrid materials that contain both organic (methyltriethoxysilane, MTEOS, T species) and inorganic functions (tetraethylorthosilicate, TEOS, Q species), combining for instance hydrophobicity and high mechanical stability (90). T and Q stand for the number of oxygen on each silicon, namely Tri (3) and Quadri (4), for MTEOS and TEOS, respectively. The letter is associated with a superscript indicating the number of condensed Si-O-Si bridges. It is important to quantify the ratio T / Q and the condensation degree, mainly by 29 Si MAS NMR, in order to properly characterise such hybrid materials. This nucleus suffers from a low natural abundance of 4.7 % and an intermediate resonating frequency at 1 / 5 th of 1 H one, both lowering SNR. It is thus current to average noise over one night or one weekend for a single spectrum. Using denoising is an interesting approach to decrease acquisition time.
Every chemical was used as received with no further purification. The solution was stirred at room temperature for at least one hour at 500 rpm to ensure hydrolysis of the precursors. Controlled condensation occurred during spray 7 drying of the sample, performed using a mini spray dryer B-290 (BUCHI) fitted with an atomiser (nozzle tip diameter = 0.7 mm) and a peristaltic pump. The temperatures at the inlet and outlet of the spray dryer were fixed at 220 °C and within the range of 95-120 °C, respectively. Polydisperse spherical particles of hybrid organic / inorganic amorphous silica (characteristic size: 1-10 µm) were obtained. 29 Si solid-state NMR experiments were performed on a Bruker Avance III spectrometer complex points and cosine multiplication were applied after SVD. This apodisation limits both signal truncation and broadening effects. One may note that SVD was not directly applied to spectra (SPC, frequency domain) because zero-filling increases matrix size and thus computation time.

B.3. Simulation of kinetics studied by Raman spectroscopy under Python
2000 Raman spectra of 2000 points each were calculated in silico using four Gaussian lines at 450, 510, 750 and 900 cm -1 , respectively. Full Widths at Half Maximum 8 (FWHM) ranged from 118 to 212 cm -1 . Such high FWHM are typical of amorphous materials like glasses (45). In order to reflect a kinetic evolution, the amplitude of peaks at 510 and 750 cm -1 were linearly decreased across the series while amplitude of peaks at 450 and 900 cm -1 were linearly increased. This series of spectra may mimic ageing of a material for instance. Homoscedastic white Gaussian noise was added on each spectrum. SVD was applied using Principal Component Analysis (PCA) function from Python Scikit-learn package (92). Computation took only a few seconds under Python Anaconda 3.5. The source code is available in file Figure_I.4a.py of (93).

B.4. Simulation of NMR spectra with known noise under Matlab
NMR complex FID were simulated in silico under Matlab (The MathWorks, Inc., Natick, MA, USA) with a complex exponential at the expected frequency ν and either an exponential decay (Equation 1) or a Gaussian decay (Equation 2), leading to a Lorentzian or a Gaussian peak on spectra after Fourier transform, respectively. While the former is typical of a relaxation-driven shape, the latter highlights a distribution of chemical environments (94) or more complex relaxation phenomena, e.g., strong dipoar coupling. To obtain a Gaussian peak with the same FWHM as a Lorentzian peak, the shape is defined according to Equation 3.
7380 NMR FID were simulated, grouped as follows: 9 • 2 shapes for decay: exponential and Gaussian; • 3 T 2 values of 10, 1.0 and 0.10 ms corresponding to narrow (32 Hz), intermediate (320 Hz) and broad peaks (3200 Hz), respectively, which are typical values obtained by 13 C, 29 Si or 31 P solid-state NMR, for various types of crystalline or amorphous materials; • 41 levels of homoscedastic white Gaussian noise ranging from -20 dB to +20 dB; • 30 random noise patterns at the same noise level.
Additionally, each data set was repeated 6 times with different processing parameters: • 2 with truncation or not, at 5 T 2 , time above which signal is almost no longer existent; • 3 Significance Level (SL, see Section D.3) for SVD automatic thresholding at error level of 5, 7.5 or 10 %.
Each FID was composed of 1024 points for a duration of 49 ms. SVD was applied before zero-filling (if truncation was applied) and Fourier transform. As peak was not at the middle of the spectrum, signal region was defined as the 512 points centred at peak frequency. Noise region corresponded to the other 512 points. Baseline zero-order offset was preliminary corrected by subtracting the mean value of noise region. This step was essential to avoid spectrum aliasing due to Fourier transform, especially for broad peaks. The source codes of SVD automatic thresholding and FID simulations are available in files sfa.m and Figure_I.7_I.8_I.S2.m of (93), respectively.
Computation of the full set of 6 × 7380 = 44280 spectra took 30 minutes with an overclocked Intel Core i5 4670K @ 4.4 GHz processor with Matlab R2016b.

C. Theoretical background
In this section, SVD and low-rank approximation are first developed. Hankel and Toeplitz matrices are then presented. Finally, SNR is defined.

C.1. SVD and low-rank approximation
SVD is a mathematical tool used to decompose a matrix X with m rows and n columns, whatever its size or shape, into the product of three other matrices U, Σ and V T (Equation 4). This is illustrated in Figure 1 by orange hatched rectangles. U and V are unitary square matrices, of size m × m and n × n, respectively. If complex numbers are used, V T , the transpose of matrix V, is replaced by V*, its conjugate transpose. SVD can indifferently be applied on real or complex matrices, the only difference being a double computation time for complex matrices (see part (II) of this work (86)). The central matrix Σ has the same shape as the original matrix X. Nevertheless, it has values only on its main diagonal (green rectangle), sorted by amplitude. These diagonal entries are called singular values and are the non-negative square roots of the eigenvalues of X T X or XX T (95). One can notice that the more elongated is matrix X, the less singular values it has. Using low-rank k, matrix X can be approximated to X k according to Equation 5, where U k , Σ k and V k T are the matrices U, Σ and V T truncated at k values, represented as blue filled rectangles in Figure 1. This process is really useful to compress large data sets (68). In the case of noise-containing data, true signals correspond to low-k values while noise signals are related to high-k values. Thus, selecting the correct k-limit allows to keep all true signals while rejecting noise, including t 1 -noise (73).
Additionally, when a baseline distortion is present, its intense signal correspond to the first singular value which may be removed (96). The following section describes how to apply SVD on one-dimensional (1D) data.

C.2. Hankel and Toeplitz matrices
As stated above, SVD can only be applied to matrices. However, 1D data form only a (complex) row or column but not a matrix. In such a case, a transformation step is required. This can be performed thanks to Hankel matrix, or similarly to Toeplitz matrix. The former is defined by its first row and its last column. All anti-diagonal values are filled identically to the first ones ( Figure 2). Toeplitz matrix is defined by its first column and its first row, all diagonal values being identical (97). Circulant matrices are a special case of Hankel or Toeplitz matrices, where every row of the matrix is a cyclic shift of the row above (98). However, in general case, these are semi-circulant matrices, with one value being replaced from one row to the next one.  Normalisation over (n-1) points is preferred to avoid a bias in standard deviation (99).
This formula is valid only when signal can be measured without any noise, so-called pure signal. However, the only observable parameter on a real noisy signal is the signalplus-noise-to-noise ratio (SNNR) (100) defined by replacing σ signal by σ signal+noise in Equation 6. SNR can then be deducted from SNNR, following Equation 7.
The other possible definition is the analytical chemistry formula (Equation 8), where PSNR is the SNR based on peak amplitude (100 ). While SNR is related to the area of the studied peak, PSNR is related to its 14 height, leading to different results. They can be expressed in decibels (Equation 10), which is more convenient to explore a wide variation range.
Additionally, Currie defined a Critical Level (L c ), a Detection Limit (L d ) and a Quantitative Limit (L q ), at 1.64, 3.29 and 10 σ noise , respectively (103). The corresponding levels, SNR and PSNR measurements are drown in Figure 3.

D. Results and discussion
In this section, SVD is applied to NMR and Raman spectra. We then focus on practical aspects of denoising, namely the impact of the matrix shape and the number of components used for thresholding. Time and frequency denoising are compared. The minimum experimental SNR needed for valid use of SVD and its impact on peak shape are thoroughly investigated. Finally, a limit case is evaluated.

D.1. Denoising of NMR and Raman spectra
Images are already matrices and SVD can directly be applied on them. Twodimensional (2D) spectra can be treated similarly. However one-dimensional (1D) spectra are not directly suitable for SVD. If a series of spectra is available, one just need to stack the successive spectra to obtain a 2D data set. Figure 4a shows such a stack of spectra simulating a reaction kinetics studied by Raman spectroscopy, as described in Section B.3. A similar stack can be obtained when a surface is mapped to analyse species in a sample region (31). The spectrum consisted of four overlapped peaks with varying intensities. Without SVD, the two peaks at 510 and 900 cm -1 were difficult to detect due to the amount of noise (red arrows). However, after SVD, these components were identified as evidenced by green arrows. One should note the very low residual noise. 16  It should be noted that this spectrum already displays a very good SNR. The signal has been enhanced using CP, allowing a non-quantitative spectrum to be acquired in 40 minutes. The second spectrum shows the same sample analysed using CPMG echoes (39), which led to numerous spikelets. The overall shape of the CPMG spectrum 17 is qualitatively similar to the above spectrum. Such an approach was proved to be a mean to increase SNR during acquisition step by discretising broad peaks (104). The original shape with improved SNR can be recovered by summing echoes, but nevertheless this leads to relaxation distortions. The third spectrum of Figure 4b is a vertical zoom to highlight noise level and much less intense spikelets marked by red arrows. By comparison with the same spectrum after SVD processing, such small peaks were highly enhanced and noise has disappeared. This 29 Si CPMG MAS NMR spectrum was used hereafter as a reference to evaluate the performances of the SVD process.

D.2. Matrix shape
It is sometimes argued that efficient denoising can be obtained using an iterative process on a rectangular matrix, with a number of columns higher than the number of signals of interest (105), i.e. roughly the number of peaks. The iteration consists in converting into a matrix, applying SVD, and reverting to a 1D set of denoised data. Due to (anti-)diagonals averaging (see Section C.2), the matrix at the beginning of the second iteration is not exactly the one at the end of the first iteration, which explains that multiple iterations give different results. However in our case, such a procedure led to some residual noise ( Figure S1). Even with m × n = 3901 × 128 points and 10 iterations, noise was only marginally reduced. Nevertheless, this procedure has the advantage of a low computation time per iteration, as Hankel or Toeplitz matrix size is smaller for a rectangular shape than a square shape, when starting with the same number of points in 1D data set (see part (II) of this work (86)).
In a second step, Figure 5a depicts the efficiency of denoising on various matrix shapes, either rectangular ones (m × n = 3997 × 32 points) or square ones (m × n = 2015 × 2014 points). While noise was strongly present for elongated matrices, it decreased when the number of columns increased and finally disappeared for a square matrix. Square matrices were used hereafter. More generally, our results suggest to tend to a square data matrix before applying SVD, as also recommended by Van Huffel et al. (106). For n = 512 and n = 2014, small peaks seem to be missing, highlighted by orange arrows. This feature is explained in next section.

D.3. Thresholding
Another parameter to be varied in order to optimise denoising is the number of singular values k (corresponding to signals) used for low-rank approximation (Figure 5b). In a first attempt, k was set to the number of peaks present, k = 22. However, k = 25 resulted in a better shape for the three more intense spikelets, corresponding to Q 3 and T 3 peaks  and did not display artefacts. It should be noted that singular value thresholding is also available (109) but it needed to adjust too many parameters and to sparsify data.

D.4. Time and frequency domains
In this section, spectra resulting from a same FID with SVD denoising applied either before ( Figure 6  Noisy and denoised spectra are shown in Figure 6b. Interestingly, very different behaviours were observed. Satisfactory denoising was achieved using mult_SPC or Toep_FID. This result explains why SVD was efficient in the case of the afore mentioned Raman spectra mimicking a kinetic process (Figure 4a) and the NMR FID ( Figure 4b). On the contrary, for Toep_SPC, noise -though reduced-was still present.
This clearly denotes that in case of 1D frequency domain, a reverse Fourier transform before converting data into a Toeplitz matrix and applying SVD should be preferred (110). 22

D.5.a. Comparison of SNR dB and PSNR rms dB
An additional question is the sensitivity of SVD, that is to say the minimum SNR needed to get a proper signal detection. The corresponding indicators, SNR dB and PSNR rms dB , were defined in Section C.3. 7380 noisy spectra were simulated with either Lorentzian or Gaussian shape and with a peak width of either 32 Hz, 320 Hz or 3200 Hz (see Section B.4). First, we compared in Figures 7a and 7b    a vertical asymptote at PSNR rms dB = 37 dB. This feature was a consequence of the spectral extension of the wings of Lorentzian peaks that were contributing to amplitude within noise region. On the contrary, Gaussian peak wings are much less intense and in this case, this vertical asymptote was not observed on SNR dB = f(PSNR rms dB ) evolution.

D.5.b. Automatic thresholding
SVD and Malnowski's SL automatic thresholding (see Section D.3) were applied to these simulated NMR FID corresponding to only one peak. When SL error level was set to 5 % and above PSNR rms dB = 20 dB this single peak was detected with k = 1 singular value, whatever peak width and shape (Figures 7c-d and S3a-b, dashed vertical black line). However, this was only an upper limit and many peaks were detected around PSNR rms dB = 17 dB, between the detection limit of 3.3 σ noise and the quantification limit of 10 σ noise , as defined by Currie (103). Thus, to be detected through Malinowski's algorithm, a signal has to be enough different from noise, i.e. between two to three 25 times higher than noise max . Surprisingly, a second singular value was detected with k = 2 for Gaussian shapes above PSNR rms dB = 36 dB, (Figures 7d and S3b). The amplitude of this second component was significant and improved the resulting shape of denoised spectrum. This second singular value will be explained in Section D.6.c. On the full data set of 7380 spectra, only one false detection was observed, as indicated by the black arrow in Figure 7c. Thus, the amount of artefacts was negligible. On the truncated FID ( Figures S2c and S2d), this limit of PSNR rms dB = 20 dB was not so abrupt, due to lack of accuracy on noise measurement. However, an advantage of truncation was a much faster computation for broad peaks, thanks to the smaller matrix used.
When SL error level was set to 7.5 or 10 %, the minimum SNR to get a peak detection was decreasing ( Figures S3c-f). However, the number of false detections also increased noticeably, as indicated by the black arrow. In some rare cases, evidenced with lines being higher than the figure vertical limit (black arrow on Figure S3d), SL was unable to distinguish signal from noise, resulting in a noisy spectrum after SVD.
While an error level of 5 % is really safe, a level of 7.5 % may be necessary to detect tiny peaks. A value of 10 % seems too high to avoid artefacts. Results are summarised in Table 1. SNR is not presented in this table as it is not a relevant parameter, that is too much depending on peak width and shape.

D.5.c. Error measurement
The difference between denoised signal and simulated non-noisy signal (pure signal), 26 was measured using the Root Mean Square Deviation (RMSD), defined on Equation 11, where y i denoised and y i pure are the individual values for denoised and pure SPC, respectively (Figures 7e and 7f). A high RMSD was obtained below PSNR rms dB = 17 dB. As no peak was detected in this range, the obtained value corresponded to RMSD of pure signal compared to zero, which was higher for broad peaks, due to its wide spread range.
Above this PSNR rms dB value, RMSD displayed a steep decrease under 0.1 (dashed horizontal black line). These results emphasised the very good agreement between denoised and pure data. However, RMSD was higher for Gaussian than for Lorentzian peaks (Figure 7f). Above the second threshold of PSNR rms dB = 36 dB, RMSD exhibited a further decrease downto the level obtained for Lorentzian peaks. This confirms the significance of the second singular value. A similar trend was observed on truncated data ( Figures S2e and S2f).

D.6.a. Pure and denoised spectra
The next step that was investigated concerns the possibility to use denoised spectra for the sake of quantification. For each peak width and shape at PSNR rms dB = 20 dB, i.e. at quantification limit, the spectrum with the worst RMSD is presented in Figure 8. These spectra were fitted with a Voigt function with error estimation implemented in Matlab (111). Amplitude, position, shape and width were automatically adjusted. Results are reported in Table 2 and 3 for Lorentzian and Gaussian peaks, respectively. A very large uncertainty occurred on fitting parameters derived from noisy spectra (top traces). For 27 denoised spectra (middle top traces), uncertainty decreased significantly, roughly by a factor of 10, but was still higher than for pure spectra (middle bottom traces). Difference between denoised and pure spectra is presented in bottom trace of Figure 8.

D.6.b. Lorentzian and Gaussian peaks
Surprisingly, despite a RMSD lower than 0.1, the area Percent Error (PE area given by Equation 12) could be as high as 8.5 % and 42.6 % for Lorentzian and Gaussian peaks, respectively. While the former was acceptable at detection limit, the latter evidenced an overestimation. Although PE area decreased after SVD denoising on Lorentzian peaks, it increased for Gaussian peaks. Moreover, difference spectrum on Gaussian peak exhibited a mix of narrow and wide components with opposite amplitudes (Figure 8f).
Such a shape modification was not observed for Lorentzian peaks. Besides, the Gaussian / Lorentzian ratio was around 0.5 instead of 1.0 after denoising on Gaussian peaks (dark grey and light grey rows of Table 3). This result highlighted that SVD induced a change in peak shape from Gaussian peaks to more Lorentzian ones. Above the second threshold of PSNR rms dB = 36 dB, the shape was corrected thanks to the second singular value, giving a pure Gaussian peak after denoising (white rows of Table   3).

D.6.c. Real and extracted errors
The error on the measured area can origin from two sources: first the error introduced by the added noise, known as the real error, and second the error coming from the denoising itself, so-called the extracted error (108). An example of real error is presented in Figure S4a. Two successive measurements with NS = 120 scans gave a different amplitude for Q 2 peak at -91 ppm (red arrow). With a higher noise averaging at NS = 360 (not shown), amplitude ratios were similar to NS = 840. Thus, SNR at NS = 120 was too low and amplitude was tainted by error. A strong apodisation has been used here to artificially improve SNR (see Section D.7.a). Automatic thresholding was unable to correctly discriminate signals from noise and manual thresholding with k = 5 singular values was preferred. Nevertheless, the real error was kept after denoising ( Figure S4b), which demonstrated that manual thresholding is a dangerous tool. Failure of automatic thresholding is thus an indication that SNR has to be improved.
The extracted error was especially present for Gaussian spectra, for which the Gaussian / Lorentzian ratio modification (see Section D.6.b) led to an increase of peak area. Indeed, SVD fits time decays with a sum of exponential (55). When fitting a Gaussian decay with a single exponential component, corresponding to one singular value at low PSNR rms dB , the peak area is correspondingly overestimated by 20 % (112).
Our results were consistent with this value. Gaussian and exponential decays are very different, as Gaussian is flatter around its maximum. Above the second threshold, the Gaussian decay was fitted with two exponential decays, improving peak area value. into account for analysis of our data, PE area was found to be similar between denoised spectra at PSNR rms dB = 20 dB (dark grey rows in Table 2 and 3) (101) and noisy spectra at PSNR rms dB = 30 dB, (light grey rows in Table 2 and 3). In solid-state NMR, another possibility to avoid Gaussian peak error relies on use of CPMG echoes (39) as in Figure   4b. This technique transforms a peak driven by chemical shift distribution (Gaussian shape, inhomogeneous interaction) into multiple narrow peaks driven by relaxation (Lorentzian shapes, homogeneous interaction) (113,114), which are very suitable for SVD, being both sensitive and quantitative.

D.7.a. Pre-processing
A pre-processing step called apodisation can be applied on FID before SVD and Fourier transform. The aim is first to reduce noise, and second to remove truncation artefacts leading to oscillations at peak foot, hence the name. In NMR, one can use for instance either exponential, cosine, or (shifted-)Gaussian decays. Their shape were compared in (115). While exponential is concave, cosine is convex and Gaussian is intermediate. In Figure S4c, we compared the influence of apodisation on initial noise. The resulting denoised SPC with automatic thresholding at an SL error level of 7.5 % are presented in Figure S4d. Without apodisation, SNR was too low to detect Q 2 peak at -91 ppm (red circle) and k = 4 singular values were found. With cosine apodisation, the correct number of peaks was obtained with k = 5 singular values (green ellipsis). With an exponential decay of 20 Hz, corresponding to the intrinsic SPC resolution, a similar result was obtained, but with k = 6 singular values, leading to small baseline distortions.
Surprisingly, with an exponential decay of 50 Hz, a much higher number of k = 94 33 singular values was found, with almost no denoising (orange circle). An explanation was that apodisation changed the amplitude of noise values, especially at the end of the FID. By this way, noise became heteroscedastic, decreasing efficiency of SVD and Malinowski's criterion. When plotting singular values in logarithmic scale ( Figure S4e), the slope moved from a plateau for cosine (purple curve) to a decay for exponential with 50 Hz (blue curve). In such a case, it is harder to discriminate signals from noise, as the slope is similar. Cosine is thus a good compromise before SVD as it decreases noise without changing too much singular values. An alternative would be to combine SVD and Savitzky-Golay smoothing filter (116), which process noise the same way all over the FID.

D.7.b. Denoising
On the sol-gel 50 / 50 MTEOS / TEOS sample, four hours and NS = 240 scans with cosine pre-processing, were needed to have a spectrum with a sufficient SNR to apply SVD and to detect Q 2 peak with automatic thresholding at SL error level of 10 % (not shown). If SL error level was limited to 7.5 %, six hours and NS = 360 scans were necessary (Figure 9 top trace). The corresponding denoised spectrum (middle top trace) was very close to the reference spectrum acquired in fourteen hours and NS = 840 (middle bottom trace), with PSNR rms = 9.7, i.e. at quantification limit. Their difference (bottom trace) was comparable to noise. This was confirmed by peaks integration (Table 4) where very good agreement was obtained between the denoised spectra at NS = 240 or NS = 360 and the reference noisy spectrum at NS = 840. In particular, Q 2 / Q 3 ratio was very consistent. Time gain was thus between 2.3 to 3.5, depending on the SL error level allowed.

E. Conclusion
Singular Value Decomposition is of crucial importance in many mathematical treatments involved in spectrocopies. In this first part (I), SVD with low-rank 35 approximation was successfully applied to denoise NMR and Raman spectra. This approach can easily be generalised to other spectroscopies. We have shown that a better denoising was obtained with square matrices and with SVD applied to time domain signal rather than to the corresponding frequency spectrum. Automatic thresholding was used thanks to Malinowski's Significant Level indicator and a 7.5 % error value was a good compromise between sensitivity and unwanted artefacts. 6×7380 SVD were carried out to compare pure, noisy and denoised spectra with SNR dB ranging over 41 dB. Our results proved that this technique can detect signals as low as twice noise max , i.e. with PSNR max = 2.0 and PSNR rms = 6.6, whatever the peak width. A systematic shape modification has been highlighted for Gaussian peaks with an overestimation of peak area by 20 %. This overestimation for Gaussian peaks is a major result as peak shape is often neglected when denoising, which can give misinterpreted data. A correction step is thus needed if Gaussian / Lorentzian ratio of denoised peak is around 0.5. When used carefully, SVD can lead to similar results between denoised spectra at PSNR rms = 6.6 and noisy spectra at quantification limit (PSNR rms = 10). As PSNR rms is increasing with the square root of time, this difference is equivalent to a considerable gain on acquisition time of 2.3, which is of paramount importance for low sensitivity experiments.
In a second part (II) (86), we will focus on the computation time needed for SVD treatment under Java, Matlab and Python, using both processors and graphic cards.
We will check the influence of algorithms, especially the divide and conquer one, as well as the influence of single precision calculation will be investigated. Software libraries such as MKL (Intel Math Kernel Library) and hardware capabilities such as SSE3 (Streaming SIMD Extensions) (118) will be evaluated. All these optimisations will decrease computation time by a factor of 100. Figure S1: Influence of SVD iterations; Figure S2: SVD applied to 7380 simulated spectra with truncation; Figure S3: Influence of significance level for automatic thresholding; Figure S4: influence of number of scans and preprocessing; simulated data sets and source codes are available online in (93).

G. Acknowledgements
Sylvie Masse and Cédric Lorthioir are thanked for fruitful discussions.