Enhancing Compression Level for More Efficient Compressed Sensing and Other Lessons from NMR Spectroscopy

Abstract

Modern nuclear magnetic resonance spectroscopy (NMR) is based on two- and higher-dimensional experiments that allow the solving of molecular structures, i.e., determine the relative positions of single atoms very precisely. However, rich chemical information comes at the price of long data acquisition times (up to several days). This problem can be alleviated by compressed sensing (CS)—a method that revolutionized many fields of technology. It is known that CS performs the most efficiently when measured objects feature a high level of compressibility, which in the case of NMR signal means that its frequency domain representation (spectrum) has a low number of significant points. However, many NMR spectroscopists are not aware of the fact that various well-known signal acquisition procedures enhance compressibility and thus should be used prior to CS reconstruction. In this study, we discuss such procedures and show to what extent they are complementary to CS approaches. We believe that the survey will be useful not only for NMR spectroscopists but also to inspire the broader signal processing community.

1. Introduction

Nuclear Magnetic Resonance spectroscopy (NMR) is currently one of the most versatile techniques of chemical and physical analysis. Its range of applications is impressively broad: from analysis of small molecules structures in all states of matter [1], through characterization of complex natural mixtures [2,3], including applications to medical screening (metabolomics) [4] up to the biological studies of structure and dynamics of proteins and ribonucleic acids [5]. The introduction of Fourier transform (FT) in 1966 [6] became a cornerstone of modern NMR spectroscopy, which is based on a measurement of a free induction decay signal (FID) in a time domain. The FID is induced in a receiver coil of an NMR spectrometer by oscillating effective magnetization of nuclear magnetic moments polarized by an external magnetic field and excited by a radio frequency (RF) pulse. Importantly, the precession frequency is dependent not only on the nuclear magnetic moment and the induction of an external magnetic field but also on the electronic surrounding of a nucleus causing shielding or deshielding effect. Thus, the frequencies emitted by the sample are interesting for chemists that can deduct molecular structures from them. Formally, the precession frequency ($\omega$) is dependent on the external magnetic field ($B_0$), shielding tensor ($\sigma$) and magnetic moment being a product of gyromagnetic ratio $\gamma$ and spin vector I.

$$\omega =-\gamma I\left(𝟙-\sigma \right)B_0$$

(1)

The FID signal $s\left(t\right)$ typically takes a form of a sum of oscillatory decaying components:

$$s\left(t\right)=\sum _{k=1}^K{A}_{k}exp\left(i({\omega }_{k}t+{\varphi }_k)\right)$$

(2)

The number of components K is equal to the number of groups of nuclei differing in resonance frequency. Each component has its amplitude ($A_k$) and frequency (${\omega }_k$). The imaginary part of a frequency corresponds to a decay rate of a signal (transverse relaxation). The phase error $\varphi$ stems from various experimental imperfections and is typically either constant or linearly dependent on the frequency ($\varphi ={\varphi }_0+{\varphi }_1\omega$).

The concept of multidimensional data acquisition [7] opened way to measurement of N-dimensional FID signals that are functions of several time variables $t_1,t_2,\dots ,t_N$. Such signals are built of products of N components similar to those in Equation (2):

$$s(t_1,t_2,\dots ,t_N)=\sum _{k=1}^K\prod _{l=1}^N{A}_{k}exp\left(i({\omega }_{k,l}{t}_l+{\varphi }_{k,l})\right),$$

(3)

where the k - index in ${\omega }_{k,l}$ corresponds to the component’s number and l corresponds to the dimension of the signal (there are N dimensions in total, one direct and other indirect).

N-dimensional spectra contain useful information—not only about resonance frequencies ${\omega }_{k,l}$ but also about interactions exploited to trigger excitation transfer between nuclei. This allows resolving the structure of a studied molecule, i.e., to determine which nuclei are connected by single or multiple chemical bonds (transfer via spin-spin couplings), which are close in space (transfer via dipole-dipole cross-relaxation) etc. [1].

However, the acquisition of multidimensional NMR data is very time-consuming, which is cumbersome due to high costs of NMR hardware maintenance, chemical instability of some samples [8] and processes occurring in them [3]. The problem of lengthy acquisition stems from combination of three facts. First, according to the Nyquist-Shannon sampling theorem [9], the sampling rate must be at least equal to the bandwidth of a signal. Secondly, the spectral resolution is proportional to the maximum time sampled [10]. Both requirements must be fulfilled in all spectral dimensions, which means that a number of data points grows exponentially with the dimensionality of a spectrum, reaching many thousands. Finally, every sampled data point in indirect dimensions ($t_1,t_2\dots t_{N-1}$) is acquired as a separate, one-dimensional FID signal ($s\left(t_N\right)$). The excited spin system needs to recover (or at least approach) its equilibrium state before next point is acquired, which usually takes up to a few seconds. Multiplied by several thousand indirect dimension points, this leads to even days-long NMR experiments.

Many methods have been proposed to alleviate the problem of lengthy sampling in multidimensional NMR experiments. Currently, the majority of them are based on the concept of sparse non-uniform sampling, where certain sampling points are removed from the sampling schedule and reconstructed mathematically based on various assumptions about the resulting spectrum. The assumptions may include: maximum entropy of a spectrum [11], presence of empty regions [12], minimum number of FID components [13] or minimum number of meaningful spectral points (“sparsity”) [14,15,16]. The latter assumption is a central pillar of the compressed sensing (CS) method that conquered many branches of technology and science [17], including chemical sciences [18]. The sparsest spectrum is found by minimizing the penalty function which is a sum of two terms: the first measures its accordance with the measured data and the second expressed by the ${\ell }_p$-norm ($0<p\le 1$) of the spectrum corresponds to the spectrum sparseness. The minimum can be found by algorithms like iterative soft thresholding [14,19] (IST) or iteratively re-weighted least squares [20] (IRLS). Approaches related to famous CLEAN method also resemble ${\ell }_p$-norm minimization [21]. The concise presentation of the CS theory can be found below in Section 2.

Other techniques include the projection spectroscopy based on co-sampling of several indirect time dimensions [22], covariance spectroscopy based on non-Fourier analysis of the conventionally sampled data [23], extrapolation of such data using linear prediction [24] or attempts to remove aliasing from sampling below Nyquist rate [25]. Although their effectiveness is also based on the compressibility of the spectrum, the relationship to compressed sensing is loose and thus they are out of scope of this study.

In this work, we will examine the relation between the sampling level and compressibility of a spectrum in the context of various NMR experiments. We survey the data acquisition and signal processing tricks that enhance the compression level and show, for the first time, that it is related to the amount of data required to obtain credible spectral reconstruction with CS methods. The relation, although stems from CS theory, has never been practically verified and demonstrated to the NMR community. On the other hand, readers from outside NMR field may get inspired by the experimental tricks enhancing the compressibility and use their analogues in different contexts. Thus, we find it beneficial to share lessons from NMR spectroscopy with experts of the broadly defined signal processing field.

2. Theory

Compressed sensing theory is based on the concept of sparse representation of a signal and compressibility of a signal. These two features depend on a chosen basis of a signals’ vector space $\mathbf{V}$. The basis ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_n$ is usually referred to as a dictionary. The examples of dictionaries that on the one hand are often used by practitioners and on the other are covered by the CS theory are Fourier basis and wavelet basis [26,27]. As shown by Qu et al. [16], the former is more efficient in the case of NMR spectra. Given a dictionary ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_n$ we say that a signal $\mathbf{s}\in \mathbf{V}$ is k-sparse if it can be written as ${\lambda }_1{\mathbf{v}}_1+\dots +{\lambda }_n{\mathbf{v}}_n$ where at least $n-k$${\lambda }_j$-coefficients are equal to zero. The set of k-sparse signal is denoted by ${\Sigma }_k$. Please note that by simple combinatorics ${\Sigma }_k$ is a union of $\left(\genfrac{}{}{0pt}{}{n}{k}\right)$ vector subspaces of $\mathbf{V}$ but it is not a convex subset of $\mathbf{V}$.

Having the latter in mind let us formulate the central problem of CS theory. For simplicity of formulation, we shall consider the case of $\mathbf{V}={\mathbb{C}}^n$ for which we fix one of the standard dictionaries (e.g., the Fourier dictionary is very useful in the NMR context). For a class of k-sparse signals $\mathbf{s}=(s_1,s_2,\dots ,s_n)\in {\Sigma }_k$ the CS theory determines the minimal number of the coordinates $s_j$ of the signal $\mathbf{s}$ required for the effective determination of the remaining ones. The effectivity requirement excludes the brutal search strategy, which is just checking all possible l-sparse signals with $l\le k$ satisfying the measurement constraint and finding the sparsest one among them. Non-effectivity of this approach is recognized by estimating a time needed for this strategy to be implemented in case of the standard size signals, e.g., $n=512$ and small sparsity, e.g., $k=10$. The time required for solving such a problem by a brute-force would be of the order of hundreds of years (see e.g., an estimation in [28], p. 54). The CS theory provides a useful alternative for a brutal search strategy by replacing a non-convex problem (recall the non-convexity of ${\Sigma }_k$) by its convex version.

Before describing the CS approach in more details, let us note that the sparsity is undoubtedly too strong condition from the practitioner standpoint. Indeed, a signal acquired in a given experiment is contaminated with a measurement noise and the noisy part of a signal excludes its sparsity with respect to the standard dictionaries. Nevertheless, CS is still useful due to the approximate sparsity of a signal, in this work referred to as compressibility.

Qualitatively speaking a signal is compressible if it can be well-approximated by a sparse representation. This feature can be expressed quantitatively by means of the distance of the signal from the subset of the k-sparse signal ${\Sigma }_k$. The measure of this distance can be chosen in many different ways. One possibility is to use the ${\ell }_p$-norm and define

$${\sigma }_k{\left(\mathbf{s}\right)}_{p,{\mathbf{v}}_i}=\underset{\widehat{\mathbf{s}}\in {\Sigma }_k}{inf}{\parallel \mathbf{s}-\widehat{\mathbf{s}}\parallel }_{p,{\mathbf{v}}_i}$$

(4)

where

$$\parallel \sum _i{\lambda }_i{\mathbf{v}}_i{\parallel }_{p,{\mathbf{v}}_i}={\left(\sum _i{\left|{\lambda }_i\right|}^p\right)}^{\frac{1}{p}}.$$

If $\mathbf{V}={\mathbb{C}}^n$ and $({\mathbf{e}}_1,{\mathbf{e}}_2,\dots ,{\mathbf{e}}_n)$ is the canonical basis of ${\mathbb{C}}^n$ then the norm ${\parallel ·\parallel }_{p,{\mathbf{e}}_i}$ will be denoted ${\parallel ·\parallel }_p$:

$$\parallel {({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)}^T{\parallel }_p={\left(\sum _{i=1}^n{\left|{\lambda }_i\right|}^p\right)}^{\frac{1}{p}}.$$

Expanding a signal $\mathbf{s}$ in a dictionary ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_n$

$$\mathbf{s}={\lambda }_1{\mathbf{v}}_1+\dots +{\lambda }_n{\mathbf{v}}_n$$

the best kth approximation is given by keeping k largest components from the coordinates $({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)$ and putting the others to zero.

Compressed sensing provides a methodology for

• highly probable exact recovery of a sparse signal based on limited information about it;
• highly probable approximately exact recovery of a compressible signal based on limited information about it.

High probability refers to the fact that recovery of the signal provided by CS method is (approximately) exact with very high probability, i.e., a wrong recovery is possible but very improbable.

To describe the convexification of CS problem let us consider the function ${\parallel ·\parallel }_{0,{\mathbf{v}}_i}$

$$\mathbf{V}\ni \sum _i{\lambda }_i{\mathbf{v}}_i\mapsto \sum _i{\left|{\lambda }_i\right|}^0\in \mathbf{N}$$

returning the number of non-zero components in the $\left({\mathbf{v}}_i\right)$-expansion. This is often referred to as ${\ell }_0$-norm (actually not being a norm in the mathematical sense). Fix $m\le n$ together with a subset $J\subset \{1,\dots ,n\}$ of cardinality m and for every $j\in J$ fix $s_j\in {\mathbb{C}}^n$. Using ${\ell }_0$-norm one formulates the main CS problem as follows:

$${\mathrm{minimize}\parallel \mathbf{x}\parallel }_{0,{\mathbf{v}}_i}\mathrm{subject}\mathrm{to}{x}_j=s_j\mathrm{for}\mathrm{all}j\in J.$$

(5)

The CS theory is based on the fundamental insight that the ${\ell }_0$-problem as formulated in Equation (5) can be replaced by its ${\ell }_1$-version

$${\mathrm{minimize}\parallel \mathbf{x}\parallel }_{1,{\mathbf{v}}_i}\mathrm{subject}\mathrm{to}{x}_j=s_j\mathrm{for}\mathrm{all}j\in J$$

(6)

at least for a certain class of basis [29]. To be more precise, the solution to the problem Equation (5) coincides with the solution of Equation (6) and for this to hold there must exist sufficiently k-sparse solution $\mathbf{s}$, and the measurement matrix must satisfy so-called uniform uncertainty principle (UUP) also known as restricted isometry property (RIP). The measurement matrix $\mathbf{A}={\left(a_{ji}\right)}_{j\in J,i\in \{1,\dots ,n\}}$ assigned to a basis ${\left({\mathbf{v}}_i\right)}_{i\in \{1,\dots ,n\}}$ of ${\mathbb{C}}^n$ and a fixed measurement schedule $J\subset \{1,\dots ,n\}$, $\left|J\right|=m$ is a matrix $\mathbf{A}$ with m rows and n columns such that for all $j\in J$

$$s_j=\sum _{i=1}^n{a}_{ji}{\lambda }_i.$$

We say that $\mathbf{A}$ satisfies uniform uncertainty principle for k-sparse vector if there is a constant $\delta >0$ (sufficiently small) such that

$${(1-\delta )\parallel \lambda \parallel }_2\le {\parallel \mathbf{A}\lambda \parallel }_2\le (1+\delta ){\parallel \lambda \parallel }_2$$

(7)

for any $\lambda =\left({\lambda }_i\right)\in {\mathbb{C}}^n$ which has at least $n-k$ zero coordinates, or in other words when ${\lambda }_1{\mathbf{v}}_1+\dots +{\lambda }_n{\mathbf{v}}_n\in {\Sigma }_k$. If this condition is satisfied with sufficiently small constant $\delta$ and if the k-sparse solution exists then the solutions of Equation (5) and Equation (6) are equal.

Restricted isometry property for $\mathbf{A}$ is desired by practitioners. For specific dictionaries, this condition can be highly probable in the sense that given the sparseness level k and a random measurement schedule J of size m, RIP holds with high probability for sufficiently large m. An example of a precise criterion in the case of Fourier basis was given in [26] where the authors proved that for every $\rho >0$ there is a constant $D_{\rho }$ such that if

$$m>D_{\rho }klog\left(n\right)$$

(8)

then the random schedule J leads to the measurement matrix $\mathbf{A}$ which satisfies uniform uncertainty principle with probability $1-\mathcal{O}\left(n^{-\rho }\right)$. In particular, we can control the level of RIP-probability of the partial Fourier transform by choosing sufficiently large m, and the good news is that the number m of the measurement of the required signals, is linear in k up to the $log\left(n\right)$ component.

The above discussion considers the idealized case of a sparse vectors’ recovery by CS method. However, the practice, in particular in NMR, immediately leads to non-sparse CS context. The first reason is the measurement error (noise) which forces the strict equality $x_j=s_j$ in Equation (6) to be replaced by an equality up to a certain error $x_j\approx s_j$. Usually ≈ is expressed by ${\ell }_2$-norm of the vector ${(x_j-s_j)}_{j\in J}$ in ${\mathbb{C}}^m$

$$\parallel {(x_j-s_j)}_{j\in J}{\parallel }_2={\left(\sum _{j\in J}{|x_j-s_j|}^2\right)}^{\frac{1}{2}}.$$

(9)

Moreover, in a certain areas of applied CS, for instance in NMR, the strict sparseness assumption must be replaced by the compressibility of signals even in the noiseless case. Indeed, the standard Lorentzian peaks present in an NMR spectrum have infinite support in frequency domain. These two facts justify the replacement of the strictly sparse CS problem Equation (6) by its relaxed approximately sparse noisy version

$${\mathrm{minimize}\parallel \mathbf{x}\parallel }_{1,{\mathbf{v}}_i}\mathrm{subject}\mathrm{to}{\parallel {(x_j-s_j)}_{j\in J}\parallel }_2\le \eta$$

(10)

where $\eta$ reflects the level of the measurement errors. As proved within CS theory [30], the reconstruction algorithms are stable, i.e., for small error level and for the compressible vector $\mathbf{s}$ the solution of Equation (10) is close to $\mathbf{s}$. To be more precise this happens if the measurement matrix has the restricted isometry property—the error of the solution of ${\ell }_1$-CS problem is measured by ${\sigma }_k{\left(\mathbf{x}\right)}_{1,{\mathbf{v}}_i}$, see Equation (4). In other words the k largest $\lambda$’s in the expansion $\mathbf{x}={\lambda }_1{\mathbf{v}}_1+\dots +{\lambda }_n{\mathbf{v}}_n$ will be recover with high accuracy, see [29].

CS in the NMR context is often applied for the recovery of the spectra consisting of peaks with significantly different amplitudes (e.g., an NMR spectrum may contain a dominating peak whose amplitude may be even four orders of magnitude larger than the amplitudes of other peaks, see NOESY spectra). In such a case, many points contributing to the bottom part of the large (Lorentzian!) peak can be “more significant” (have higher values) than smaller peaks and thus be reconstructed in the first place. From the spectroscopist point of view, the hierarchy of importance is opposite—the small, “off-diagonal” peaks in NOESY carry the most important information. Such non-linearity of the reconstruction is also the reason the signal-to-noise ratio (SNR) in NMR spectra reconstructed with CS is not informative [31,32]. Depending on the reconstruction parameters the apparent noise level can vary from zero (high sparsity enforced) to far higher than actual (too low sparsity enforced causing incomplete removal of “NUS artifacts”). This effect can be seen in all Figures below.

To summarize this part, the concept of approximate sparsity and its relationship with the amount of data to be measured is crucial in the context of NMR experiment. Thermal noise, dynamic range of peak intensities and their linewidths, the fact that signal is complex (but only real part is of interest) create a specific framework for the application of CS to reconstruct missing points in the FID signal. Keeping this in mind, we move to the practical considerations—a bunch of Lessons about the effective use of CS in experimental NMR.

3. Results and Discussion

3.1. Lesson 1: Reduce the Number of Peaks

As can be seen from inequality Equation (8), the number of sampling points required for an efficient CS reconstruction (m) is dependent on the number of important spectral points (k). In the language of NMR spectroscopy, k is, roughly speaking, the number of points contributing to peaks. Thus, experimental techniques reducing the number of peaks to the necessary minimum not only improve spectral resolution but also allow reconstruction of the spectrum with lower m and shorten the experimental time.

Many such techniques were proposed. “Pure-shift NMR” (PS-NMR) certainly belongs to the most spectacular improvements over the last decade [33,34]. The idea of PS-NMR is to remove the multiplet structure of NMR spectra by suppressing the effect of J-couplings i.e., interactions between nuclear spins transferred within the molecule via the nearby chemical bonds. The splittings of peaks caused by J-couplings can be informative but lead to the reduction of resolution and requirement of more sampling points for the proper CS reconstruction (since more peaks are present in a spectrum). With the splittings removed, less data points are required for the reconstruction (see Equation (8)). The simulation from Figure 1 shows this effect. The fact that PS-NMR naturally fits to the compressed sensing reconstruction has been discussed extensively by Aguilar and coworkers [35]. Importantly, while the couplings between nuclei of the same kind (i.e., homonuclear) can be removed using selective echo pulse sequence block in both direct and indirect dimensions of an NMR spectrum, the pseudo-random NUS is feasible only in the latter case. In the direct dimension, pure-shift experiment can be performed by sampling “chunks” of an FID signal. However, as shown in several papers [36,37,38,39] such data can also be used as an input for CS algorithms, although the RIP (cf. Equation (7)) is worse and thus more sampling points are required. Interestingly, such “burst sampling” has been reported by some authors to have also certain advantages over other sampling schemes in the indirect dimensions [40].

Besides pure-shift approach, the number of spectral peaks can also be reduced by more selective coherence transfer in correlation experiments. The selectivity is achieved by adjusting the delay time $\Delta$ in the coherence transfer block or additional coherence-selection delay. The coherence of spins coupled with J-constant evolves in an oscillatory manner, typically as ${sin}^n\left(\pi J\Delta \right)$ where n is the number of nuclei J-coupled to the nucleus from which the transfer starts. The classic example is multiplicity selection [41,42] which exploits the fact that coherence transfer between interacting nuclei A and X is dependent on the multiplicity (n) of a spin system $A{X}_n$. For instance, CH groups can be selectively excited in a 2D ¹H-¹³C Heteronuclear Single-Quantum Correlation spectrum (HSQC), making CH₂ and CH₃ peaks invisible. Similarly, one can make the transfer selective by exploiting the differences in J between various pairs of nuclei. In 3D HNCA, the basic experiment used to establish sequential connectivities in spectra of proteins, the excitation is initially transferred from amide hydrogens to amide nitrogens. Then, from each ¹⁵N nucleus the transfer may go two-fold—to $\alpha$ carbon of the same (i) and preceding ($i-1$) amino acid residue. This is caused by the fact that ^HN-C_α coupling constants for both ways are similar (typically 11 Hz for N_i-C_αi and 7 Hz for N_i-C_αi−1) and $\Delta$ can be set to average value in-between. However, the variants of the experiment with an exclusive transfer to one C_α also exist and have been used in combination with non-uniform sampling [43], also due to better compressibility of the spectrum and thus the sampling.

Another approach to achieve the reduction in a number of peaks in protein spectra was proposed by Dötsch and colleagues [44]. They modified a CBCA(CO)NH pulse sequence [45], to acquire a signal for amino acid types selected basing on topology. Only desired amino acid residues give signals in such spectra, which facilitates the analysis and allows efficient low-level non-uniform sampling [46].

Other compressibility-enhancing pulse sequence blocks allow the triggering of an exponential signal decay due to diffusion or relaxation and suppress signals selectively due to differences in the decay rate. The diffusion filter is based on the gradient echo block added to the standard NMR pulse sequence [47]. Used for mixtures of chemical compounds, it suppresses the signal from quickly diffusing smaller molecules (although it exists also in a reverse mode [48]). The somewhat opposite effect is achieved by a $T_2$-filter (Carr–Purcell–Meiboom–Gill sequence block [49]) which suppresses signals from nuclei with short transverse relaxation times (typically belonging to larger molecules). Diffusion-filtering and multiplicity selection, as well as their effect on the required number of sampling points, are shown in practice in Section 3.6.

3.2. Lesson 2: Minimize Dynamic Range

The high dynamic range of amplitudes of signal components does not constitute a significant problem for CS reconstruction in the case of signals with strictly sparse representation. The real FIDs, however, contain noise and are represented by Lorentzian peaks in the Fourier domain, with their half-width being not negligible, but proportional to the signal decay rate. The consequence is an imperfect performance of CS algorithms which are usually based on an iterative deconvolution of a point spread function (PSF) from the spectrum (for the meaning of PSF in NMR context see [50,51,52,53]). For example, one of the most classical CS algorithms, the orthogonal matching pursuit (OMP), does this by seeking for the FT basis function (an “atom”) giving the highest inner product with the FID. Then the approximation that uses only that atom is subtracted from the signal and the process is repeated. Other algorithms, like iterative (soft or hard) thresholding are based on a very similar concept [21]. The noise obviously disturbs the approximation and makes artifact removal less complete. Additionally, as mentioned in Section 2 above, the algorithm will rather tend to improve the bottom points at the sides of large Lorentzian peaks than reconstruct lower components. Thus, whenever possible, it is crucial to avoid high dynamic range of components in the spectra reconstructed with CS.

Very large and very tiny peaks are found together in spectra of mixtures of chemical compounds. It might happen, however, that some of the compounds are not interesting for the spectroscopist and can be suppressed in a spectrum. This is easy to achieve if interesting and non-interesting compounds differ significantly in the molecular size. As mentioned above, the diffusion and relaxation filters can be useful for this purpose. Unfortunately, the small and large peaks may be found even in the pure, single-compound samples. This is the case of spectra based on nuclear Overhauser effect (NOESY and ROESY). The additional difficulty arises from the fact that peaks intensities (especially those of small peaks) are the most informative parameters and thus must be reconstructed with high fidelity. This is troublesome, as series of tiny off-diagonal peaks are accompanied by huge (even ${10}^4\times$ larger diagonal peaks. To deal with this difficulty the diagonal-free NOESY experiments have been proposed [54] and shown to be particularly effective when combined with non-uniform sampling [55,56,57,58].

The Figure 2 shows this effect on a 1D cross-section along the indirect dimension of a simulated NOESY spectrum.

Notably, reducing the dynamic range of a measured object is beneficial in CS applications other than NMR. In a nice example of identification of bacterial species in a mixture by a single Sanger-sequencing reaction [59], Amir and Zuk suggested taking a square root of the data to reduce the differences between “peaks”.

3.3. Lesson 3: Pre-Processing

The reconstruction of missing points in NUS data from NMR experiments is one of the middle steps in the data processing workflow. This, usually, starts from the conventional procedures performed in the direct dimension: filtering, apodization, zero-filling, phasing, Fourier transform and baseline correction [60]. Importantly, before the CS reconstruction of the indirect dimension points is performed, one can apply procedures that will make the frequency representation more compressible. These include: removal of the imaginary part of a spectrum (virtual echo, VE) [61], removal of assumed modulation in the FID (virtual decoupling, VD) [62,63,64] and combination of in-phase and anti-phase (IPAP) complementary sub-signals [65,66].

The first trick, virtual echo, is based on the fact that phase of the signal in the indirect dimensions ($\varphi$ in Equation (3)) is usually known apriori, and thus the signal can be phased before the reconstruction. Thus, imaginary part of the spectrum is not needed for CS procedure and can be removed. This is beneficial, as an imaginary part of the Lorentzian line has dispersive shape and is less sparse than absorptive real part. The effect is achieved by combining FID signal with its conjugate mirror-reflected “copy” in each dimension [61]. Alternatively, a similar effect can be achieved by modifying the minimized term in Equation (10) to calculate only the real part of $\mathbf{x}$ [67]. The Figure 3 demonstrates the concept of VE.

The virtual decoupling is the removal of the manifestation of scalar couplings, i.e., cosine modulation of the FID, by dividing the signal by the assumed $cos\left(\pi Jt\right)$ function [63,64] or by equivalent modification in the algorithm [62]. The operation makes the spectrum sparser but is based on two rigorous assumptions. First, all FID components must share the same modulation (i.e., the same J). Secondly, division by zero must be avoided—either by regularization or by omitting zeros in the sampling schedule [68]. Contrary to virtual echo, the virtual decoupling is beneficial even in the case of fully sampled data, where no reconstruction is required. It leads to resolution and sensitivity enhancement as broad multiplets collapse into narrower and higher singlets. However, the requirement of constant J among all spin systems is rarely fulfilled. Typical examples are limited to adjacent carbons in 3D HNCA [64] (C_α-C_β coupling) or HC-CH TOCSY spectra [63] (coupling between carbon atoms belonging to methyl and a neighboring group).

The NMR spectra are sometimes combined from subsets, as is done in case of IPAP method [65,66]. Two experiments are recorded, first providing doublets in-phase and the other anti-phase. Then, they are added which cancels the doublet components with opposite sign. In standard FT processing, it does not matter whether the addition is performed on FIDs or on spectra. For CS, however, the former is more beneficial, since it makes the spectrum more compressible (reduces the number of peaks) before the reconstruction.

The reduction of the number of components by sample preparation, sophisticated signal excitation or pre-processing is not limited to NMR and can enhance CS reconstruction in other fields. The inspiring examples can be found in CS video processing where common features of neighboring frames can be found by motion-estimation helping to enhance sparsity [69]. The sub-regions of interest in a reconstructed object can be also explicitly defined to reduce the number of significant points, as shown in the field of angiography [70].

3.4. Lesson 4: Match Sampling with the Decay

The vital aspect of NUS in the indirect dimensions of an NMR experiment is that the sampling schedule can be completely arbitrary. For example, its density can be modulated according to the assumed function. Besides J-modulation mentioned above, the useful trick is to avoid large gaps in the sampling schedule [50,71] by approaches known from other fields like jittered sampling [72] or Poisson-disk sampling [73]. Some authors show that gaps should be avoided at the beginning and at the end of a signal [74].

The oldest and most commonly applied modulation of NUS density is relaxation-matched sampling introduced by Barna et al. [75]. The gains on signal-to-noise ratio (SNR) have been analyzed in detail by other authors [31,32,76]. The reason for SNR improvement is quite simple—an FID signal decays exponentially in time, while noise level remains constant. Thus, initial sampling points have higher SNR. The measurement sensitivity aspect of the problem is simple. However, the situation becomes more complicated when analyzed from the point of view of CS theory. It can be shown that although relaxation-matching improves SNR, it worsens restricted isometry property of the measurement matrix [77]. Thus, there is certain balance between the two effects. This fact is also connected to the observation that CS works more efficiently for the interpolation of the data rather than extrapolation [78].

Interestingly, the sensitivity benefit from relaxation-matched CS can be so strong that it can even lead to results better than fully sampled experiment acquired in the same time. The simulation in Figure 4 shows this effect.

In many applications outside NMR field, the sampling density can be adjusted “on-the-fly”, i.e., during the measurement. Such adaptive sampling is well established in image processing where sampling density is a function of local image variance [79,80]. In NMR spectroscopy such examples, although feasible, are still lacking.

3.5. Lesson 5: Non-Stationarity

The parameters of typical NMR signals (frequencies, amplitudes, relaxation rates) do not vary in time when measured for the stable samples. Sometimes, however, at non-stationary conditions, the frequency in some spectral dimensions varies in time [81]. Since sampling of the indirect dimensions can be arbitrary (e.g., pseudo-random), the effective frequency-time dependencies can be complicated. As demonstrated [8,81], the FID frequency varying in time leads to lineshape distortions in the case of “chronological sampling”, i.e., ($t=0,\Delta t,2\Delta t,3\Delta t\dots$) and noise-like artifacts in the case of “shuffled” sampling. Figure 5 shows that NUS of a non-stationary signal lead to spectral quality even better than fully sampled data. This can be explained by the fact that frequency variations within an FID, occurring, e.g., due to chemical reaction in the sample, are reduced due to shorter time needed for data collection. This means that compressed sensing should be the method of choice for samples whose state varies over the time of experiment.

In the scientific literature, one can find several interesting examples of application of NUS/CS in NMR spectroscopy for monitoring of the physical/chemical processes. The monitoring of processes involving biomolecules (e.g., proteins) is particularly interesting. In a paper by Bermel et al. [83], the NUS NMR experiment accompanied by CS reconstruction was successfully applied to monitor protein dynamics in a function of temperature. The application of NUS and CS in their work allowed to precisely track peak positions and intensities during sample heating. The reactions occurring in complex mixtures were also extensively investigated by NUS/CS NMR spectroscopy [3,84]. Usually, such mixtures require at least 2D NMR experiment to resolve important peaks in a spectrum, which makes monitoring troublesome using conventional sampling. One may also apply NUS/CS NMR spectroscopy to track different chemical reactions when a good temporal resolution and the benefits provided by 2D NMR experiments are required [85,86].

Obviously, the fact that NUS/CS experiments are quick compared to full sampling provides benefits in disciplines other than NMR. A good example has been discussed by Vasanawala et al. [87], where authors applied CS in pediatric MR imaging. Since children are rarely able to stay still during the measurement, the undersampled data is often of better quality (resolution) despite the need for reconstruction.

3.6. Practical Example

In this section, we finally move to the experimental example of a 2D HSQC spectrum acquired with NUS and reconstructed using CS (see Figure 6). The 2D HSQC is one of the main workhorses of structural identification, acquired in huge numbers in the most of NMR labs. Thus, optimization of its speed and quality is important.

As described in Section 3.1, NMR offers a variety of different pulse sequence “blocks” that can reduce the number of observed peaks in a spectrum to the necessary minimum and hence increase its compressibility. In this section, we verify the relation between the number of peaks in a spectrum (compressibility) and the reconstruction quality. For that purpose, we acquired 3 variants of the ¹³C HSQC experiment characterized by a different number of peaks in a spectrum. The acquired 2D NMR signals for each HSQC variant were artificially sub-sampled by taking out random points from the full data in the t₁ (¹³C) dimension and reconstructed back to the original size. The reconstructed spectra from the corresponding sub-sampled HSQC experiments are depicted in Figure 6: standard unedited ¹³C HSQC (Figure 6b), ¹³C HSQC with CH-only editing (Figure 6c) and ¹³C HSQC with CPMG filter (Figure 6d). A fully sampled, unedited ¹³C HSQC spectrum is also depicted in Figure 6a) and stands as a quality reference for the reconstructed spectra (Figure 6b–d). The ¹³C HSQC NMR experiments used in this study employ the same core HSQC pulse sequence [88], which allows observing single-quantum ¹H–¹³C correlation signals. The use of appropriate filters (multiplicity-editing—Figure 6c and CPMG—Figure 6d) to the core HSQC sequence (Figure 6a,b) allowed us to reduce the number of components in the signal. The filters were chosen concerning the physicochemical properties of the substances being measured. The sample used for experiments was a mixture of sucrose and heparin dissolved in D₂O. Both compounds are saccharides, but their molecular weights (MW) differ significantly, as heparin is a polysaccharide of MW in the range from 6000 up to 20,000 g/mol, while sucrose is a disaccharide of MW = 342.3 g/mol. We used this fact to suppress signals of fast-relaxing nuclei belonging to large heparin molecules by means of CPMG relaxation-filter (Figure 6d). We also used the fact that structures of sucrose and heparin consist mainly of CH and CH₂ chemical sites, which are the source of ¹H–¹³C single-quantum correlation signals. We employed the multiplicity-editing block to suppress signals that arise from CH₂ chemical sites, thus, only the signals corresponding to CH sites were visible. The unedited ¹³C HSQC spectrum experiment (Figure 6a—reference and Figure 6b—reconstructed spectrum) show all the single-bond ¹H–¹³C correlations regardless of the molecular size and type of chemical site.

The benefits of using editing and filtering “blocks” in 2D NUS NMR experiments can be found through the comparison of stacked spectra in Figure 6. The numerous ¹H–¹³C correlation signals in the unedited ¹³C HSQC spectrum (Figure 6a) were poorly reconstructed using 24 out of 256 t₁ sub-samples (Figure 6b). The effect is visible on heparin signals near 3.55/75.0 ppm and 3.65/80.0 ppm (marked with the black arrows in Figure 6). A reduction of the number of peaks in a spectrum allowed for more reliable reconstruction using the same 24-points sampling level for ¹³C HSQC with CH-only editing (Figure 6c) and ¹³C HSQC with CPMG filter (Figure 6d).

4. Conclusions

Due to high maintenance costs of high-resolution NMR spectroscopy, it is beneficial to apply sparse sampling techniques in multidimensional measurements and save the experimental time. However, since the number of sampling points required for the efficient reconstruction grows with the “complexity” of a spectrum (number of peaks and dynamic range of intensity) it is recommended to minimize it before the CS reconstruction. This can be achieved using dedicated acquisition and processing techniques. It is also noteworthy that in some cases, like strongly decaying or non-stationary signals the sparse sampling followed by the reconstruction leads to results superior to full sampling followed by Fourier transform. In this work we summarized those, often unnoticed, aspects of compressed sensing in NMR.