1. Introduction
Compressed sensing (CS) theory provides the methods to solve an underdetermined linear equation y = Ax within a set of sparse vectors x ∈ ℂ^{n}. The unconstrained linear equation for an m × n matrix A of a maximum rank may have infinitely many solutions for m < n. However, the number of sparse solutions happens to be finite. Vectors in ℂ^{n} having at most l ∈ ℕ nonzero coordinates are referred to as the lsparse vectors.
All lsparse vectors fulfilling the equation y = Ax may be found by considering $(\begin{array}{c}n\\ l\end{array})$ linear subproblems of solving y = Ax with the constraint x ∈ ℂ_{I}, where ℂ_{I} ⊂ ℂ^{n} denotes the subspace of vectors supported by I ⊂ [1,n]. This provides all lsparse solutions, but the overall problem was proven to be NPhard. Remarkably, the sparsest vector x̂ solving y = Ax may be found as the solution of the l_{1}minimization problem:
A number of CS algorithms has been proposed, including: iterative soft thresholding (IST) [2], orthogonal matching pursuit (OMP) [3], compressive sampling matching (COSAMP) [4], stagewise orthogonal matching pursuit (StOMP) [5] and iterative reweighted least squares (IRLS) [6]. Inspired by the greedy CS algorithm, OMP, and motivated by a case of nuclear magnetic resonance (NMR) spectroscopy, we propose an approach that is designed specifically to recover a signal being a combination of exponentially decaying components. The spectrum of such a signal consists of Lorentzian peaks. The algorithm matches a single peak in each of its iterations, and thus, the name Lorentzian peak matching pursuit (LPMP) is proposed. The peaks' centers and heights are found as in the OMP algorithm, whereas the widths of the peaks are determined in a novel way.
2. NMR Spectroscopy and CS
The sources of the free induction decay (FID) signal measured in NMR spectroscopy are magnetic moments of atomic nuclei with nonzero spin, e.g., ^{1}H,^{13}C,^{15}N,^{19}F, that undergo coherent precession when excited by radiofrequency irradiation in the presence of a high static magnetic field [7]. The time dependence of an ideal FID signal f(t) is given by the following linear combination:
The precession frequencies (and thus, the frequencies of the emitted FID signal) are in the order of MHz, but their band is very narrow, typically a few kHz. In such a narrow range, tens or even hundreds of peaks are present. Often, the spectral resolution (peak width compared to the difference in resonance frequencies) is too low to resolve peaks and, thus, makes an analysis difficult. Fortunately, the peak dispersion can be improved by the acquisition of the multidimensional NMR signal [8] where additional (“indirect”) time dimensions are introduced. Peaks in the multidimensional NMR spectrum correspond to groups of nuclei coupled by some physical interaction. Using the above notation, the formula of the twodimensional FID is:
The reason for the application of CS in multidimensional NMR spectroscopy is a timeconsuming sampling of the indirect dimensions of the FID signal (t_{1} in Equation (2)). There are three reasons for the lengthy acquisition. Firstly, each measurement point in the indirect dimension takes a few seconds of “real” time. Secondly, the sampling density has to fulfil the Nyquist theorem [9]. Finally, the sampling has to reach far into the indirect time domain to provide sufficient spectral resolution [10]. The latter, together with the Nyquist condition for the sampling density, means that, often, many thousands of sampling points have to be acquired. This makes an experiment unacceptably long (days in the case of three or four dimensions) or completely impossible in the case of higher dimensionality (5D and more). Thus, there is a need for alternative sampling approaches that allow one to save expensive spectrometer time preserving the spectral information. Most of them use the nonuniform sampling, where the major part of the points is not measured, but reconstructed after the experiment based on certain assumptions [11,12]. The CS, where the assumption is spectral sparsity, has been successfully applied in NMR exploiting the algorithms based on convex l_{1}norm or nonconvex l_{p}norm minimization [13,14]. The history of the application of the “greedy” algorithm, OMP, in NMR is long; in fact, the algorithm was transferred from radio astronomy, where it existed under the name, CLEAN [15]. Direct application of CLEAN in NMR is possible [16], but far from being perfect, because of the aforementioned approximate sparsity of the signal. Thus, various modifications of CLEAN have been proposed [1719]. Importantly, none of these works refer to OMP formalism and do not explain the modifications by the sparsity criteria. We propose an approach that directly exploits the fact that the NMR spectrum consists of a low number of Lorentzian peaks. Thus, LPMP is based on the sparsity concept, which is more relevant in the case of NMR than classical CS methods.
3. OMP Algorithm
Let us fix the notation. For m,n ∈ ℤ and m < n, we write [m, n] = {m, m + 1,…, n}. Let m, n ∈ ℕ. For matrix A ∈ M_{m×n}(ℂ), its range and kernel are denoted ran A and ker A, respectively. The Hermitian conjugate of A is denoted by A*. The number of nonzero coordinates of x ∈ C^{n} is denoted by ‖x‖_{0}, and supp x (the support of x) denotes the set of x nonzero coordinates. Note that ‖x‖_{0} = # supp x (the cardinality of a set I is denoted by # I). We define χ_{n,l} ⊂ ℂ^{n}:
OMP algorithm arose in the context of the sparse approximation problem; see [3]. In order to formulate the problem, let us fix a matrix A ⊂ M_{m×n}(ℂ) satisfying ran A = ℂ^{m}. The matrix A is in this context referred to as the dictionary matrix.
Definition 1
Adopting the above notation, we say that x ∈ χ_{n,l} is an optimal lsparse approximate representation of a vector y ∈ ℂ^{m} if:
We say that y is lrepresentable, if there exists a vector x ∈ χ_{n,l}, such that y = Ax.
For a subset I ⊂ [1, n] of cardinality l, we define a vector subspace ℂ_{I} ⊂ ℂ^{n}:
The restriction of A to ℂ_{I} is denoted by A_{I} : ℂ_{I} → ℂ^{m}.
Definition 2
We say that A ⊂ M_{m×n}(ℂ) is linjective if ker A_{I} = {0} for all I ⊂ [1, n], # I = l.
For an linjective matrix A and I ⊂ [1, n], # I = l, we define:
Since x_{I} is a vector satisfying A_{I}x_{I} = y_{I}, where y_{I} is an orthogonal projection of y onto ran A_{I}, we may see the uniqueness of x_{I} for linjective A. Let us note that x ∈ χ_{n,l} is an optimal lsparse approximate representation of y ∈ ℂ^{m} if:
Although linjectivity of the matrix A does not exclude a vector x′ ∈ χ_{n,l}, x′ ≠ x, such that ‖y — Ax‖_{2} = ‖y — Ax′‖_{2}, we shall usually choose an optimal lsparse approximate and denote it by x_{l}.
The OMP algorithm provides a certain approximation of x_{l} ∈ ℂ^{n}, which we denote by ${x}_{l}^{\text{OMP}}$ Remarkably, ${x}_{l}^{\text{OMP}}$ is found on the basis of ${x}_{l1}^{\text{OMP}}$ Denoting the ${x}_{l}^{\text{OMP}}$ support by ${I}_{l}^{\text{OMP}}$, the index i_{l} ∈ [1, n] is given by:
The OMP vector ${x}_{l}^{\text{OMP}}$ is defined as ${x}_{l}^{\text{OMP}}={x}_{{I}_{l}^{\text{OMP}}}$; see Equation (3) for the definition of ${x}_{{I}_{l}^{\text{OMP}}}$. For the coherence, restricted isometry property (RIP) conditions that imply the equality ${x}_{l}^{\text{OMP}}={x}_{l}$, we refer to [3] and [20], respectively. For the recovery limits of the OMP, see [21]. The OMP algorithm (Algorithm 1) is schematically described below.
Algorithm 1 Orthogonal matching pursuit. 

4. Free Induction Decay Signals in NMR
The free induction decay signal in NMR consists of the combination of exponentially decaying oscillations (see Equations (1) and (2)). The signal can be multidimensional, but for the sake of simplicity, we will limit ourselves to the onedimensional indirect part of the twodimensional signal, corresponding to f(t) from Equation (1). Given the damping factor γ ≥ 0 and the oscillation frequency ω ∈ (x0211D), a single decaying oscillation is given by:
Its Fourier transform:
The FID signal is sampled at the rate given by the Nyquist theorem for the range cut off ν ∈ [–Ω, Ω]. In order to perform the NMR signal processing, the frequency resolution ΔΩ is introduced, where $\mathrm{\Delta}\mathrm{\Omega}=\frac{\mathrm{\Omega}}{2N+1}$ and N ∈ ℕ. For any l ∈ [–N,N], we get the basic complex oscillation e^{2πιωlt}, where ${\omega}_{l}=\frac{l\mathrm{\Omega}}{2N+1}$. The NMR signals can be approximated by the finite combination of the basic complex oscillations:
The approximation being $\frac{2N+1}{\mathrm{\Omega}}$ periodic works only for $0\le t\le \frac{2N+1}{\mathrm{\Omega}}$. Assuming that the series ${\sum}_{l=N}^{N}{a}_{l}{e}^{2\pi l\omega lt}$evaluated at ${t}_{k}=\frac{k}{\mathrm{\Omega}}$ coincides with sampling points of f, i.e., f(t_{k}) for k ∈ [0,2N], we get:
Introducing the time unit $\frac{1}{\mathrm{\Omega}}$ and the frequency unit $\frac{\mathrm{\Omega}}{2N+1}$, we may enumerate the time and frequency by the integers [0, 2N] and [–N, N], respectively. Then, the relation between the peak intensities a and the time samples f is given by a = Ff, where F is the Fourier transform matrix:
l ∈ [–N, N] and k ∈ [0, 2N]. In order to model the Lorentzian peaks within our framework, let us introduce the width peak unit Γ and a damping matrix C ∈ M_{(2}_{N}_{+1)×(2}_{N}_{+1)}(ℂ):
Using the above notation, the FID signal is given by:
5. Lorentzian Peak Matching Pursuit
The spectrum f̂ corresponding to FID signal (Equation (5)) is given by:
We shall assume that j ∈ [0, J], which gives the maximal width JΓ and the width resolution Γ. Let K = {k_{1}, k_{2},…, k_{m}} ⊂ [0, 2N] be the sampling scheme, y ∈ ℂ^{m} the measured signal y_{i} = f(t_{ki}) and A = F_{k} the matrix consisting of rows of the Fourier matrix F enumerated by the elements of K.
The Lorentzian peak matching pursuit algorithm matches a Lorentzian peak ${L}_{l}^{j}$ in each of its iterations. In order to match a peak ${L}_{l}^{j}$ LPMP establishes first the peak's center l. The center l_{1} ∈ [–N, N] of the first peak is determined by:
In order to determine the width of the first peak j_{1} ∈ [0, J], we find bl_{1},_{j} ∈ ℂ, such that:
Then, j_{1} ∈ [0, J] is a number satisfying:
In the second step of the algorithm, we determine the second peak's center l_{2} by:
In order to determine the second peak's width j_{2} ∈ [0, J], we find b_{l}_{2}_{,j}, b_{l}_{1},_{j}_{1} ∈ ℂ minimizing:
We begin the third step of the algorithm by determining the third peak's center l_{3}:
Then, j_{3} and l_{3} are found as in the previous steps.
Note that in the kth step of the algorithm, the peak's width j_{k}Γ is found and not changed afterwards. The b_{lk},_{jk} coordinate may, in turn, vary in the nth step of the algorithm for n > k. The algorithm is performed until the accuracy threshold ‖y – Ax‖ < ε is achieved; for a different stopping criterion, see Section 5.1.
5.1. Stopping Criterion
If a number of Lorentzian peaks in the spectrum cannot be assumed a priori, then the stopping criterion for the LPMP becomes a crucial issue.
Notably, the dimension of a measurement vector y provides a mathematical bound for the number of iterations in LPMP. The implementation of the LPMP algorithm uses an inversion of a matrix, which ceases to be invertible when the number of peaks exceeds the number of measurements. In what follows, we shall introduce the stopping criteria that further restrict the number of iterations. One commonly used criterion is the threshold ε of accuracy used in Algorithm 2.
Let us denote the result of the kth LPMP iteration by f̂_{k}. The k + 1 iteration is performed if:
Alternatively, the noise amplitude σ can be used for a stopping criterion. In this case, LPMP stops on the kth iteration if f̂_{k}(l_{k}) ≤ σ.
Algorithm 2 Lorentzian peak matching pursuit. 

5.2. Mask
NMR experiments are often performed in series, with some or all peaks preserving their positions in at least some of the dimensions. This fact has been exploited to restrict the allowed peak frequencies and to improve the results of NUS (nonuniform sampling) reconstruction, e.g., in hyperdimensional spectroscopy [22], the SIFT method [23] or highdimensional spectra [24]. It is easy to use this information in the LPMP algorithm by introducing the subset Mask ⊂ [−N, N] of admissible peak centers in the frequency domain. In what follows, we give a detailed description of the masked LPMP algorithm (Algorithm 3), taking into an account the stopping criterion described in Section 5.1.
Algorithm 3 Lorentzian peak matching pursuit with a mask. 

6. Results and Discussion
In this section, we present the results of the LPMP performance tests. Extensive simulations on different levels of sampling sparsity provided a quantitative characterization of the method. As an example of an application, we have chosen the challenging NOESY (nuclear Overhauser effect spectroscopy) experiment [25] with a high dynamic range of peak intensities.
The parameters of the synthetic spectrum used in the simulations are given in Table 1. The spectral parameters were chosen to simulate various experimental scenarios, i.e., (I) peak widths differ from peak to peak, which corresponds to the differences in relaxation rates in the case of the NMR experiment; (II) the range of intensities is high; (III) two of the peaks overlap; and (IV) the peaks' centers and widths are offgrid.
An example of the LPMP recovery of an NMR spectrum based on 120 sampling points drawn at random from a 256point full grid is presented in Figure 1. The signal was contaminated with 5% of noise. We performed systemic simulations with a varying sampling level. For each level, 50 different sampling schedules were generated. Denoting the frequency spectrum vector obtained from the full sampling by f̂ and the LPMP recovery by f̂_{LPMP}, we compute the recovery error ‖f̂–f̂_{LPMP}‖, and then, we average it over 50 sampling distributions. The LPMP performance is depicted in Figure 2 (left), where it is also compared with other popular CS algorithms: IST, OMP, COSAMP and StOMP. For all sampling levels, LPMP is revealed to be superior to other CS algorithms.
In order to test the masked version of LPMP algorithm, we introduced the mask of ±5 frequency points around the peaks' centers. The simulation results are presented in Figure 2 (right). As expected, masking has a more pronounced effect on the low sampling levels.
The qualitative aspects of LPMP reconstruction can be evaluated on the data from the 2D NOESY experiment. Each reconstructed 1D slice of the 2D NOESY spectrum contains one dominating peak accompanied by a group of smaller peaks. Such a case of large intensity range is very challenging for CS methods. In Figure 3, we present the results of the LPMP spectrum recovery obtained on the basis of 200, 225 and 250 random measurements of the corresponding FID signal of length 512. The 2D NOESY spectrum was recorded on a 600MHz Varian UNITY Inova spectrometer at 25 °C using a sample of human ubiquitin protein (1 mM concentration).
In order to explain the effectiveness of LPMP algorithm in the NMR context, let us emphasize the key difference between the OMP and LPMP methods. OMP tries to obtain the frequencydomain representation having the minimal number of nonzeros, which cannot be exact, even for a single Lorentzian peak. LPMP improves OMP by replacing every determined spectral nonzero with the Lorentzian peak of the appropriate width. Roughly speaking, this corresponds to the concept of sparsity in the sense of a number of Lorentzian peaks.
7. Conclusions
Our results show that LPMP is superior to the generic CS methods when applied to exponentially decaying signals. Comparison with OMP, StOMP, IST and COSAMP revealed that the assumption of the Lorentzian spectral line shape improves the quality of the reconstruction with LPMP. Moreover, its performance can be further improved by limiting (“masking”) the support of a signal. We believe that this method will find applications in NMR spectroscopy. Other fields of possible application may include different kinds of spectroscopy, wireless communications, sonar and radar. In fact, reconstruction of any nonuniformly sampled exponentiallydamped signal can benefit from the application of LPMP.
Acknowledgments
The authors thank the Polish Ministry of Science and Higher Education for support with Iuventus Plus Grant No. IP2011 023171 and the National Centre of Science for support with Sonata Bis 2 Grant No. DEC2012/07/E/ST4/01386. Krzysztof Kazimierczuk thanks the Foundation for Polish Science for the support with the Team program.
Author Contributions
Paweł Kasprzak developed the core of the algorithm and described it. Krzysztof Kazimierczuk adapted it to NMR requirements and wrote the NMR part of the publication.
Conflicts of Interest
The authors declare no conflict of interest.
References
 Candès, E.J.; Romberg, J.K.; Tao, T. Stable signal recovery from incomplete and inaccurate measurements. Commun. Pure Appl. Math. 2006, 59, 1207–1223. [Google Scholar]
 Papoulis, A. A new algorithm in spectral analysis and bandlimited extrapolation. IEEE T Circuits Syst. 1975, 22, 735–742. [Google Scholar]
 Tropp, J.A. Greed is good: Algorithmic results for sparse approximation. IEEE Trans. Inf. Theory 2004, 50, 2231–2242. [Google Scholar]
 Needell, D.; Tropp, J. CoSaMP: Iterative signal recovery from incomplete and inaccurate samples. Appl. Comput. Harmonic Anal. 2009, 26, 301–321. [Google Scholar]
 Donoho, D.L.; Tsaig, Y.; Drori, I.; Starck, J.L. Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit. IEEE Trans. Inf. Theory 2012, 58, 1094–1121. [Google Scholar]
 Candès, E.J.; Wakin, M.B.; Boyd, S.P. Enhancing Sparsity by Reweighted l_{1} Minimization. J. Fourier Anal. Appl. 2008, 14, 877–905. [Google Scholar]
 Keeler, J. Understanding NMR Spectroscopy, 2nd ed.; Wiley: Chichester, UK, 2010; p. 526. [Google Scholar]
 Jeener, J. AMPERE International Summer School; Basko Polje, Yugoslavia, 1971. [Google Scholar]
 Nyquist, H. Certain topics in telegraph transmission theory. Trans. Am. Inst. Electr. Eng. 1928, 47, 617–644. [Google Scholar]
 Szantay, C. NMR and the uncertainty principle: How to and how not to interpret homogeneous line broadening and pulse nonselectivity. IV. (Un?)certainty. Concept. Magn. Reson. A. 2008, 32A, 373–404. [Google Scholar]
 Kazimierczuk, K.; Stanek, J.; ZawadzkaKazimierczuk, A.; Koźmiński, W. Random sampling in multidimensional NMR spectroscopy. Prog. Nucl. Magn. Reson. Spectrosc. 2010, 57, 420–434. [Google Scholar]
 Coggins, B.E.; Venters, R.A.; Zhou, P. Radial sampling for fast NMR: Concepts and practices over three decades. Prog. Nucl. Magn. Reson. Spectrosc. 2010, 57, 381–419. [Google Scholar]
 Kazimierczuk, K.; Orekhov, V.Y. Accelerated NMR spectroscopy by using compressed sensing. Angew. Chem. Int. Ed. Engl. 2011, 50, 5556–5559. [Google Scholar]
 Holland, D.J.; Bostock, M.J.; Gladden, L.F.; Nietlispach, D. Fast multidimensional NMR spectroscopy using compressed sensing. Angew. Chem. Int. Ed. Engl. 2011, 50, 6548–6551. [Google Scholar]
 Högbom, J.A. Aperture synthesis with a nonregular distribution of interferometer baselines. Astron. Astrophys. Suppl. 174. 15, 417–426.
 Coggins, B.E.; Zhou, P. High resolution 4D spectroscopy with sparse concentric shell sampling and FFTCLEAN. J. Biomol. NMR 2008, 42, 225–239. [Google Scholar]
 Stanek, J.; Kozmiński, W. Iterative algorithm of discrete Fourier transform for processing randomly sampled NMR data sets. J. Biomol. NMR 2010, 47, 65–77. [Google Scholar]
 Coggins, B.E.; WernerAllen, J.W.; Yan, A.; Zhou, P. Rapid Protein Global Fold Determination Using Ultrasparse Sampling, HighDynamic Range Artifact Suppression, and TimeShared NOESY. J. Am. Chem. Soc. 2012, 134, 18619–18630. [Google Scholar]
 Kazimierczuk, K.; Zawadzka, A.; Kozmiński, W.; Zhukov, I. Lineshapes and Artifacts in Multidimensional Fourier Transform of Arbitrary Sampled NMR Data Sets. J. Magn. Reson. 2007, 188, 344–356. [Google Scholar]
 Zhang, T. Sparse recovery with orthogonal matching pursuit under RIP. IEEE Trans. Inf. Theory 2011, 57, 6215–6221. [Google Scholar]
 Wang, J.; Shim, B. On the Recovery Limit of Sparse Signals Using Orthogonal Matching Pursuit. IEEE Trans. Signal Process 2012, 60, 4973–4976. [Google Scholar]
 Jaravine, V.A.; Zhuravleva, A.V.; Permi, P.; Ibraghimov, I.; Orekhov, V.Y. Hyperdimensional NMR spectroscopy with nonlinear sampling. J. Am. Chem. Soc. 2008, 130, 3927–3936. [Google Scholar]
 Matsuki, Y.; Konuma, T.; Fujiwara, T.; Sugase, K. Boosting Protein Dynamics Studies Using Quantitative Nonuniform Sampling NMR Spectroscopy. J. Phys. Chem. B. 2011, 115, 13740–13745. [Google Scholar]
 Kazimierczuk, K.; ZawadzkaKazimierczuk, A.; Koźmiński, W. Nonuniform frequency domain for optimal exploitation of nonuniform sampling. J. Magn. Reson. 2010, 205, 286–292. [Google Scholar]
 Kaiser, R. Use of the Nuclear Overhauser Effect in the Analysis of HighResolution Nuclear Magnetic Resonance Spectra. J. Chem. Phys. 1963, 39, 2435–2442. [Google Scholar]
Center  5.6  22.7  67.8  70.4  194.5  239.25 

width  5.2  5.4  2.6  2.8  4  1.5 
© 2015 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license ( http://creativecommons.org/licenses/by/4.0/).