Denoising NMR time-domain signal by singular-value decomposition accelerated by graphics processing units

We present a post-processing method that decreases the NMR spectrum noise without line shape distortion. As a result the signal-to-noise (S/N) ratio of a spectrum increases. This method is called Cadzow enhancement procedure that is based on the singular-value decomposition of time-domain signal. We also provide software whose execution duration is a few seconds for typical data when it is executed in modern graphic-processing unit. We tested this procedure not only on low sensitive nucleus 29 Si in hybrid materials but also on low gyromagnetic ratio, quadrupole nucleus 87 Sr in reference sample Sr(NO 3 ) 2 . Improving the spectrum S/N ratio facilitates the determination of T/Q ratio of hybrid materials. It is also applicable to simulated spectrum, resulting shorter simulation duration for powder averaging. An estimation of the number of singular values needed for denoising is also provided.


Introduction
The accuracy and efficiency of data post-processing and quantification are paramount for NMR applications.In medical applications, reliable estimates of NMR spectral parameters must be deduced from low signal-to-noise (S/N) ratio spectrum [1].I ns o l i d state NMR, time-domain signal averaging is common practice for improving S/N ratio of spectrum, but this increases NMR machine time.We apply the Cadzow enhancement procedure [2] for denoising NMR time-domain signal, resulting spectrum with better S/N ratio, which permits us to shorten NMR machine time.
This enhancement procedure is based on the singular-value decomposition (SVD) of an Hankel matrix [3] constructed with the NMR time-domain signal.Theoretical investigations focused on signal enhancement of a sum of sinusoids [2,4].This procedure is already involved in numerous NMR applications [5][6][7][8][9][10][11][12].Denoising signal becomes a preliminary step before the application of modelfitting methods to extract spectral parameters.
SVD of a large matrix is a computationally-intense operation, which limited its introduction in NMR processing software.Fortunately, modern graphics processing unit (GPU), whose original purpose is to improve graphics acceleration operations in 3D computer games, becomes programmable for general-purpose computation [13][14][15][16][17]. GPU becomes a co-processor of the main central processing unit for massive parallel computation.
In this article, we first describe the Cadzow enhancement procedure for denoising NMR time-domain signal.Then we test our SVD programs [18,19] not only to insensitive nucleus 29 Si in hybrid materials, but also to low gyromagnetic ratio, quadrupole nucleus 87 Sr (I ¼9/2, 7.02% natural abundance) [20] in reference sample Sr(NO 3 ) 2 studied by Larsen and co-workers [21] and Bowers and co-workers [22].Furthermore, 29 Si MAS spectrum with large S/N ratio facilitates the determination of T/Q ratio of hybrid materials from spectrum decomposition.The Q species describe the connectivity of the silicate network and the T species are indicative of a silicon atom from a silicate network that is bound to a carbon atom.We also denoise time-domain signals simulated with SIMPSON [23], whose line shapes are asymmetric.Thanks to modern GPU, SVD execution duration has been shorten dramatically, from 80 down to 2 s for typical solid state NMR data.

Cadzow enhancement procedure
The enhancement procedure of NMR time-domain signal involves Hankel matrix [3,24], H, in which each descending diagonal from right to left is constant: First, we fill the first row and the last column of H with the complex numbers of the time-domain signal [2,25].Fig. 1A shows a simple 4 Â 3 matrix H where the time-domain digitised data are the series [a, b, c, d, e, f] coloured in red.The remaining matrix elements are filled according to Eq. ( 1).Second, H is decomposed [1,26] as shown in Fig. 1B where the singular-value matrix Σ consists of three positive values Σ 1 , Σ 2 and Σ 3 .Third, we zero the weakest singular value Σ 3 .Fourth, we construct a new 4 Â 3 matrix with this reduced number of singular value (N SVD ) matrix as shown Fig. 1C; however, this new matrix is not of Hankel structure.Finally, we restore this matrix to Hankel one by averaging the elements of each descending diagonal from right to left as shown Fig. 1D.The denoised time-domain signal is provided by the first row and the last column of this matrix, which are coloured in red in Fig. 1D.If the corresponding spectrum is still too noisy, we repeat the procedure to the denoised time-domain signal by keeping the number of singular values smaller than or equal to the previous case.Cadzow [2] used Toeplitz matrix [27] instead of Hankel one in his presentation, in which each descending diagonal from left to right is constant.
From a practical point of view, an m Â n Hankel complex matrix represents a free induction decay (FID) whose size is N¼ mþ n À 1.The number of columns n of H should be much larger than the optimum number of singular values, because the enhanced signal strongly depends on n and N [1].
As we also want to improve the time-domain signals of quadrupole nuclei [28], which usually have asymmetric line shapes, we proceed in another way as shown the flow chart in Fig. 2. The problem is that the optimum number of singular values cannot be guessed.We first apply the line broadening (LB) processing to decrease the noise [4].After the decomposition of H, it is simpler to start with the first singular value, which has the highest value.If the spectrum resulting from the denoised timedomain signal is not satisfactory, we increase the number of singular values until the line shape is recovered.

Experimental
NMR experiments for magic-angle spinning (MAS) 29 Si in hybrid materials were performed with Bruker Avance III 500 at 99.36 MHz and 4-mm diameter rotor.Those for static 87 Sr in Sr(NO 3 ) 2 were performed with Avance III 850 spectrometer at 36.83 MHz and 9.5-mm diameter Varian rotor.

Results and discussion
Fig. 3A shows denoised 29 Si MAS spectra of a layered silicate Magadiite [29] for the number N SVD of singular values varying from 1 to 5. The spectrum is recovered with N SVD Z4, the latter being the number of 29 Si peaks in the spectrum.It consists of a Q 3 species peak located at À 99.2 ppm and three Q 4 species peaks at À 109.5, À 111.1 and À 113.7 ppm.Underestimating the number N SVD of singular values distort the spectrum (Fig. 3A-a and b) but Table 1 Relative integrated intensities of silicon species and T/Q ratio from high-power proton decoupled 29 Si MAS spectra of hybrid materials in Fig. 4, NS being the number of scans.overestimating this number makes spectrum denoising less effective (Fig. 3A-e).Fig. 3B compares line broadening (LB) denoising with SVD denoising.LB denoised spectrum (Fig. 3B-d) whose noise is comparable to that of SVD denoised spectrum (Fig. 3B-a) has broader line widths, resulting a loss of spectrum resolution.Fig. 3B-a denoised with N SVD ¼4 shows that SVD denoising neither broaden the line widths nor improve the spectral resolution.Knowing the number of lines in the spectrum helps us to choose the number N SVD of singular values in the denoising procedure.
Fig. 4A shows TPPM-15 high-power proton decoupled 29 Si MAS spectra of hybrid materials for increasing NS varying from 40 to 240 in steps of 40.Fig. 4B shows SVD denoised spectra of that acquired with 60 h of machine time or NS ¼1200 for increasing number of singular values (N SVD ) varying from 1 to 6 in steps of 1, Fig. 4B-g being the non-SVD denoised spectrum.The spectrum is recovered with N SVD Z5, the latter being the number of silicon species in the spectrum.It consists of lines of T2 and T3 species located at À 72.5 and À80.8 ppm and of Q 2 ,Q 3 and Q4 species located at À 93.5, À 102.7 and À111.9 ppm, respectively.A Q z species signifies a silicon atom bonded to z other silicon atoms via bridging oxygen and 4 -z non-silicon groups X via "non-bridging" oxygen, that is, ðXOÞ 4 À z SiðOSiÞ z .AT z species signifies a silicon atom bonded to z other silicon atoms via bridging oxygen, 3 -z nonsilicon groups R via "non-bridging" oxygen and a non-silicon group R 0 , that is, ðROÞ 3 À z Si ðOSiÞ z R 0 .Fig. 4C shows the denoised spectra of Fig. 4A where the first 6 singular values are selected, that is, N SVD ¼6.It is clearly shown that the S/N ratio of denoised spectrum acquired for 8 h machine time (Fig. 4C-d) is as good as that acquired for 60 h machine time (Fig. 4B-g).These spectra were decomposed as sum of Gaussian line shapes with GSim [30] for determining the relative proportions of the various 29 Si species.
The ratio of the integrated intensities for T to Q species, is reported in Table 1.The T/Q ratios from decompositions/ integrations of non-denoised (Fig. 4A) and SVD denoised (Fig. 4C) spectra are similar.It is $ 0.25.As SVD denoised spectra present clean spectral base lines, the initial decomposition of SVD denoised spectra before fitting (that is, base line correction, line positions, line widths, and line intensities), which is performed by us, is less tedious than those of non-denoised spectra.Error bars on integrated line intensities are not provided by GSim.Fig. 5A shows the superposition of TPPM-15 high-power proton decoupled 29 Si MAS QCPMG spikelet spectrum (Fourier transformation of time-domain echo train) of hybrid materials with that of Fig. 4B-g.Fig. 5B and C shows denoised proton decoupled 29 Si MAS QCPMG spikelet spectra of hybrid materials for increasing values of N SVD , Fig. 5C-f being the non-denoised spectrum.The spectrum is recovered with N SVD Z25, whereas the number of spikelets in the spectrum is 21 in Fig. 5C-f.From a practical point of view, the above examples show that a good guest of N SVD for SVD denoising is the number of resolved peaks in the spectrum.
Fig. 6A shows denoised static spectra of 87 Sr in the reference sample Sr(NO 3 ) 2 acquired with two-pulse, full Hahn echo sequence and 15 h of machine time, N SVD varying from 1 to 8 in steps of 1, Fig. 6A-i being the non-denoised spectrum.The line shape is recovered with N SVD Z6.Fig. 6B shows denoised (N SVD ¼6) 87 Sr spectra acquired with three machine times: 1( Fig. 6B-a), 5 (Fig. 6B-b) and 15 h (Fig. 6B-c), Fig. 6B-d  We have applied SVD denoising to experimental NMR timedomain signals.Now, we extend SVD denoising to SIMPSON simulated FIDs associated with asymmetric line shapes in static experiments.Fig. 7 shows 3 series of static chemical shift anisotropic (CSA) spectra (Fig. 7A-C) and one series of static second-order quadrupole spectra (Fig. 7D).The simulation parameters are reported in Table 2.A st h e number of singular values for denoising without distorting the line shape cannot be guessed, we check the denoised spectra for increasing N SVD number starting from 1.With N SVD ¼3 ,t h ed e n o i s e dl i n es h a p e s look like theoretical spectra but present some imperfections: the lines are not smooth.These figures also show that it is possible to apply N SVD 43 and zero some singular values inside the selected range to smooth the line shape.Applying twice the same series of singular values also smooth the line shape.With N SVD $ 5a n ds o m ez e r o e d singular values inside these 5 values, the line shapes are nearly smooth.We stop the SVD denoising procedure when the smooth line shape can superimpose the non-denoised spectrum.Without SVD denoising, the simulation of similar quality spectra with SIMPSON requires longer computer times for powder averaging.Fig. 7Es h o w s the superposition of the denoised second-order quadrupole line shape (D-g) with its line shape simulated with Dmfit [31].I ti so b v i o u st h a t SVD denoising is not only applicable to MAS spectra of spin I¼1/2 nuclei but also to asymmetric line shapes in static experiments.

Conclusion
We show that the Cadzow enhancement procedure for denoising time-domain signal is applicable not only to FID, but also to two-pulse, full Hahn echo and QCPMG echo train in NMR of spin I¼ 1/2 systems such as 29 Si and of spin I 41/2 systems such as 87 Sr.It is also applicable to simulated spectrum from powder averaging.It does not distort the line shape contrary to LB denoising which broadens the line width.But it does not improve the spectral resolution.Improving the S/N ratio of a spectrum facilitates its component decomposition to determine the T/Q ratio in hybrid materials.The number of singular values necessary for denoising is about the number of peaks in MAS spectrum.An asymmetric line shape is recovered with ca. 5 singular values.Thanks to modern GPU, typical SVD duration lasts a few seconds, which allows us to anticipate that denoising 2D spectrum should last a few minutes.

Fig. 1 .Fig. 2 .
Fig. 1.Cadzow enhancement procedure: (A) 4 Â 3 Hankel matrix H containing the digitised noisy time-domain signal series [a, b, c, d, e, f] in its first row and its last column; (B) decomposition of H where the three singular values are Σ 1 , Σ 2 and Σ 3 ; (C) the time-domain signal series is denoised by zeroing the smallest singular value Σ 3 (N SVD ¼ 2) then a new 4 Â 3 Hankel matrix is recomposed; (D) the matrix elements resulting from (c) are averaged so that the matrix is of Hankel structure; the denoised time-domain signal series is located in its first row and its last column.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. SIMPSON simulated static CSA or second-order quadrupole spectra with parameters reported in Table 2. (A) Static CSA spectra: (a) simulated spectrum; simulated spectra denoised with the first 1 (b), 2 (c), 3 (d), 4 (e) and 5 (f) singular values; (g) the first 5 singular values where Σ 4 ¼ 0; (h) twice with the first 5 singular values where Σ 4 ¼ 0; (i) superposition of (a) and (h).(B) Static CSA spectra: (a) simulated spectrum; simulated spectra denoised with the first 1 (b), 2 (c), 3 (d) and 4 (e) singular values; (f) twice with the first 3 singular values; (g) superposition of (a) and (f).(C) Static CSA spectra: (a) simulated spectrum; simulated spectra denoised with the first 1 (b), 2 (c) and 3 (d) singular values; (e) twice with the first 3 singular values; (f) superposition of (a) and (e).(D) Static second-order quadrupole spectra: (a) simulated spectrum; simulated spectra denoised with the first 1 (b), 2 (c), 3 (d) and 4 (e) singular values; (f) twice with the first 4 singular values; (g) three times with the first 4 singular values; (h) superposition of (a) and (g).(E) Superposition of Fig. 7D-g (blue line shape) with that fitted with Dmfit (pink line shape).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) being the non-denoised spectrum acquired with 15 h of machine time.The denoised line shape of Fig. 6B-b acquired with 5 h machine time is comparable with that of non-denoised line shape of Fig. 6B-d acquired with 15 h machine time.
Crystal_file zcw986 and number of points np ¼ 256 are used.b The nucleus is 43 Ca, I¼ 7/2, whose Larmor frequency is 26.96MHz. a