Robust Algorithms for the Analysis of Fast-Field-Cycling Nuclear Magnetic Resonance Dispersion Curves

Abstract

Fast-Field-Cycling (FFC) Nuclear Magnetic Resonance (NMR) relaxometry is a powerful, non-destructive magnetic resonance technique that enables, among other things, the investigation of slow molecular dynamics at low magnetic field intensities. FFC-NMR relaxometry measurements provide insight into molecular motion across various timescales within a single experiment. This study focuses on a model-free approach, representing the NMRD profile $R_1$ as a linear combination of Lorentzian functions, thereby addressing the challenges of fitting data within an ill-conditioned linear least-squares framework. Tackling this problem, we present a comprehensive review and experimental validation of three regularization approaches to implement the model-free approach to analyzing NMRD profiles. These include (1) MF-UPen, utilizing locally adapted $L_2$ regularization; (2) MF-L1, based on $L_1$ penalties; and (3) a hybrid approach combining locally adapted $L_2$ and global $L_1$ penalties. Each method’s regularization parameters are determined automatically according to the Balancing and Uniform Penalty principles. Our contributions include the implementation and experimental validation of the MF-UPen and MF-MUPen algorithms, and the development of a “dispersion analysis” technique to assess the existence range of the estimated parameters. The objective of this work is to delineate the variance in fit quality and correlation time distribution yielded by each algorithm, thus broadening the set of software tools for the analysis of sample structures in FFC-NMR studies. The findings underline the efficacy and applicability of these algorithms in the analysis of NMRD profiles from samples representing different potential scenarios.

1. Introduction

Fast-Field-Cycling (FFC) Nuclear Magnetic Resonance (NMR) relaxometry is a powerful, non-destructive magnetic resonance technique that allows the exploration of slow molecular dynamics, accessible only at extremely low magnetic field intensities. By varying the strength of the relaxation magnetic fields, FFC enables the assessment of the frequency dispersion of the longitudinal relaxation rate ($R_1$) within a sample, consequently producing NMR Dispersion (NMRD) profiles. Accordingly, FFC-NMR relaxometry measurements can detect molecular motion across a broad spectrum of timescales within a single experiment.

Under the proper constraints [1], spin relaxation theory characterizes relaxation rates as linear combinations of spectral density functions, which are the Fourier transforms of the time correlation functions. These functions capture the frequencies and their intensities present in the correlation function. Several models exist in the literature (see [2] for a survey) where the longitudinal relaxation rates are related to the parameters representing the physical phenomena occurring in the analyzed sample. Thus, all these approaches depend on the specific sample examined. There are several software tools that can be used in these cases, for example, see [3,4,5,6].

An approach that allows greater flexibility is known as the model-free approach [7], where the NMRD profile $R_1$ is expressed as a linear combination of Lorentzian functions.

The model-free approach leads to a data fit problem consisting of an ill-conditioned linear least-squares problem whose stable solution requires appropriate regularization [8,9].

We propose a review of three potential approaches: MF-UPen, which employs locally adapted $L_2$ regularization; an algorithm based on the $L_1$ penalty, referred to as MF-L1; and finally, a method called MF-MUPen, which utilizes both locally adapted $L_2$ and global $L_1$ penalties. In all these approaches, the regularization parameters are computed through automatic procedures founded on the Balancing Principle (BP) [10] and the Uniform Penalty principle [11].

All the algorithms are inspired by two-dimensional time-domain NMR relaxometry techniques. The locally adapted $L_2$ regularization was originally introduced by Bortolotti et al. in [11], where the regularization parameters are determined by applying the Uniform Penalty principle. The global $L_1$ regularization has been applied in [12], which addresses the more complex problem of data exhibiting spurious peaks caused by Quadrupolar Relaxation Enhancement. Finally, the coupling of locally adapted $L_2$ and global $L_1$ penalties was originally introduced by Bortolotti et al. in [13] for inverting two-dimensional NMR relaxation data and has been adapted to NMRD profiles.

The contributions of this work are the following:

• The implementation and experimental testing of the MF-UPen algorithm, featuring a novel rule for automatically computing the threshold parameter ${\beta }_0$.
• The implementation and experimental testing of the MF-MUPen algorithm.
• Development of a “dispersion analysis” procedure, enabling the determination of the existence range for estimated parameters.

Our aim is to illustrate the diversity of results achievable with different algorithms, focusing on fit quality and correlation time distribution. These insights will enrich the suite of software tools available for enhancing sample structure analysis in FFC-NMR studies.

Following this introduction, the mathematical problem and numerical methods are detailed in Section 2 and Section 3, respectively. Section 4 then discusses the results from testing on two sets of NMRD profiles, each representing significant potential scenarios.

2. The Parameter Identification Problem

In the following, we first describe the continuous model for NMRD profiles and then derive its discretization.

Lo Meo et al. [14] developed a heuristic algorithm to model the NMRD profiles $R_1\left(\omega \right)$ based on the model-free approach, where the profiles are defined as a combination of two terms:

$$R_1\left(\omega \right)=R_0+R^{HH}\left(\omega \right)$$

(1)

where $\omega$ is the Larmor angular frequency. The first term $R_0\ge 0$ is an unknown offset, taking into account very fast molecular motions; the second term $R^{HH}\left(\omega \right)$ describes ${ }^{1}H{-}^{1}H$ relaxation as

$$R^{HH}\left(\omega \right)={\int }_0^{+\infty }\left[\frac{\tau }{\left(1+{\left(\omega \tau \right)}^2\right)}+\frac{4\tau }{\left(1+4{\left(\omega \tau \right)}^2\right)}\right]f\left(\tau \right)d\tau$$

(2)

where $\tau$ represents the correlation time, i.e., the average time required by a molecule to rotate one radian or to move for a distance as large as its radius of gyration, and f is the distribution to be recovered. Equation (1) represents a parameter identification problem where the unknown parameters $(f\left(\tau \right),R_0)$ have to be estimated from NMRD profile values.

By discretizing Equation (1), we obtain the following linear system:

$$\mathit{y}=\mathit{K}\mathit{f}+R_0,$$

(3)

where $\mathit{K}\in {\mathbb{R}}^{m\times n}$ is derived from the discretization of the integral for $R^{HH}$ in (2):

$${\mathit{K}}_{i,j}=\frac{{\tau }_j}{\left(1+{\left({\omega }_i{\tau }_j\right)}^2\right)}+\frac{4{\tau }_j}{\left(1+4{\left({\omega }_i{\tau }_j\right)}^2\right)},i=1,\dots ,m;j=1,\dots ,n.$$

(4)

with the vector $\omega \in {\mathbb{R}}^m$ representing the m Larmor angular frequency values ($\omega =2\pi \nu$, with $\nu$ in $MHz$) at which the signal $R_1$ is evaluated. The unknown vector of the correlation times distribution $\mathit{f}\in {\mathbb{R}}^n$ is obtained by sampling $f\left(\tau \right)$ in n logarithmically equispaced values ${\tau }_1,\dots ,{\tau }_n$. The vector $\mathit{y}\in {\mathbb{R}}^m$ represents the observation vector, i.e., $y_i=R_1\left({\omega }_i\right)$, with $i=1,\dots ,m$. Generally, for a typical FFC-NMR experiment, $m\ll n$.

Defining $\mathit{x}\equiv (\mathit{f},R_0)\in {\mathbb{R}}^{n+1}$ and ${\displaystyle {\mathit{K}}_e\equiv \left[\mathit{K}1\right]\in {\mathbb{R}}^{m\times (n+1)}}$, we can reformulate (3) in a more compact form:

$$\mathit{y}={\mathit{K}}_e\mathit{x}.$$

(5)

We underline that the model-free analysis provides a fingerprint of the possible motion regimes regardless of their physical–chemical interpretation. In fact, the integral form of Equation (2) unconstrainedly gives only the number of the possible correlation times describing the dynamics of the overall physical system. Instead, the typical approach used in FFC NMR relaxometry requests an ad-hoc mathematical model containing information about the number and meaning of the correlation times that describe a given system [12].

3. The Algorithms

From a mathematical point of view, the estimation of the parameters $\mathit{f}$ and $R_0$ starting from the experimental observation $\mathit{y}$ represents an ill-conditioned linear inverse problem.

In this section, we describe the proposed algorithms based on locally adapted $L_2$ regularization (MF-UPen), $L_1$-based regularization (MF-L1), and multi-penalty consisting of local-$L_2$ and $L_1$ penalties (MF-MUPen).

3.1. MF-UPen Algorithm

This algorithm, implementing the locally adapted $L_2$ regularization, solves the following constrained minimization problem:

$$\begin{array}{cc}\underset{\mathit{x}\ge 0}{min}\hfill & \left\{{\displaystyle \parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+\sum _{i=1}^n{\lambda }_i{\left(\mathit{L}\mathit{x}\right)}_i^2}\right\}\hfill \end{array}$$

(6)

with $\mathit{L}=[\Delta ,\mathbf{0}]\in {\mathbb{R}}^{n\times (n+1)}$, where $\Delta$ is the discretization of the second derivative operator, according to central finite difference formulas, and $\mathbf{0}$ is the n−components null column vector. Observe that the regularization is imposed only on the parameter $\mathit{f}$ since the sum in (6) ranges for indices i from 1 to n. The regularization parameters ${\lambda }_i$, $i=1,\dots ,n$ are computed according to the following relaxed UPEN principle [11]:

$${\lambda }_i=\frac{{∥\mathit{y}-{\mathit{K}}_e\mathit{x}∥}^2}{n\left({\beta }_0+{\beta }_p{max}_{\mu \in I_i}{\left({\mathit{p}}_{\mu }\right)}^2+{\beta }_c{max}_{\mu \in I_i}{\left({\mathit{c}}_{\mu }\right)}^2\right)},i=1,\dots ,n$$

(7)

where $\mathit{c}=\mathit{Lx}$, $\mathit{p}=[\nabla ,\mathbf{0}]\mathit{x}$, and $I_i$ are the indices subsets related to the neighborhood of the i−th entry, i.e., $I_i=\left\{i-1,i,i+1\right\}$. The $\beta$’s are positive parameters. The parameter ${\beta }_0$ prevents division by zero and is a compliance floor, which should be small enough to prevent under-smoothing and large enough to avoid over-smoothing. The optimum value of $\beta$’s could substantially change with the nature of the measured sample.

The parameters ${\beta }_p,{\beta }_c$ are used to enhance or mitigate the local effects of slope/curvature. A preliminary trial value that often yields satisfactory results is ${\beta }_p={\beta }_c=1$. The parameter ${\beta }_0$, however, is more critical; its value should not exceed the threshold defined by

$$\underset{i}{max}\left\{{\beta }_p\underset{\mu \in I_i}{max}{\left({\mathit{p}}_{\mu }\right)}^2+{\beta }_c\underset{\mu \in I_i}{max}{\left({\mathit{c}}_{\mu }\right)}^2\right\},$$

(8)

while a too-small value, especially in cases where slope $\left(\mathit{p}\right)$ and curvature $\left(\mathit{c}\right)$ approach zero, would lead to an extremely ill-conditioned problem, hence causing computational challenges.

Therefore, we propose an automatic rule for determining ${\beta }_0$ based on the estimate of (8) obtained from a tentative solution $\widehat{\mathit{f}}$ computed by the Truncated Singular Value Decomposition [15] of the matrix $\mathit{K}=U\Sigma V^T$:

$$\widehat{\mathit{f}}=\sum _{{\sigma }_i>{\mathrm{Tol}}_{\sigma }}\frac{{\mathbf{u}}_i^{T}y}{{\sigma }_i}{\mathbf{v}}_i,{\mathrm{Tol}}_{\sigma }={10}^{-6}{\sigma }_1,$$

where ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_i\ge \cdots$ represent the singular values, and the vectors ${\mathbf{u}}_i,{\mathbf{v}}_i$ represent the i-th columns of U and V, respectively [16].

By setting

$$V_i={\beta }_p\underset{\mu \in I_i}{max}{\left({\widehat{\mathit{p}}}_{\mu }\right)}^2+{\beta }_c\underset{\mu \in I_i}{max}{\left({\widehat{\mathit{c}}}_{\mu }\right)}^2,i=1,\dots ,n,$$

where $\widehat{\mathit{c}}=\Delta \widehat{\mathit{f}}$ and $\widehat{\mathit{p}}=\nabla \widehat{\mathit{f}}$, we calculate ${\beta }_0$ as follows:

$${\beta }_0=\rho {\parallel V\parallel }_{\infty },0<\rho <1.$$

(9)

The advantage of this approach is that it substitutes the parameter ${\beta }_0$, which can range in $(0,\infty )$, with the parameter $\rho$, which is confined within the interval $(0,1)$. This substitution ensures that ${\beta }_0$ remains lower than the highest values of V but higher than the lowest ones. This makes determining ${\beta }_0$ more intuitive, particularly when supported by a visual representation of V.

To summarize, MF-UPen is an iterative scheme where, given an initial guess ${\lambda }_i^{\left(0\right)}$, $i=1,\dots ,n$, an approximate solution ${\mathit{x}}^{\left(k\right)}\equiv ({\mathit{f}}^{\left(k\right)},R_0^{\left(k\right)})$ is computed by solving (6) for fixed ${\lambda }_i^{\left(k\right)}$, $i=1,\dots ,n$, and the regularization parameters values are updated according to (7). The minimization problem (6) is solved by the Newton projection method (NP) [11]. MF-UPen is sketched in Algorithm 1; the iterations are stopped when the following condition is satisfied:

$$\sum _{i=1}^n|{\lambda }_i^{(k+1)}-{\lambda }_i^{\left(k\right)}|<\mathrm{Tol}\sum _{i=1}^n\left|{\lambda }_i^{\left(k\right)}\right|$$

(10)

where Tol is a fixed tolerance.

Algorithm 1 MF-UPen

1: Set $k=0$, and choose a starting guess ${\lambda }_i^{\left(0\right)}$, $i=1,\dots ,n$. 2: Compute ${\beta }_0$ according to (9). 3: repeat 4:     $k=k+1$ 5:     NMRD parameters update. By using the Newton Projection method compute $${\mathit{x}}^{\left(k\right)}=\arg \underset{\mathit{x}}{min}{\parallel \mathit{y}-{\mathit{K}}_e\mathit{x}\parallel }_2^2+\sum _{i=1}^n{\lambda }_i^{(k-1)}{\left(\mathit{L}\mathit{x}\right)}_i^2$$ 6:     Regularization parameter update. Set $$\begin{array}{cc}\hfill {\displaystyle {\lambda }_i^{\left(k\right)}}& =\frac{{∥\mathit{y}-{\mathit{K}}_e{\mathit{x}}^{\left(k\right)}∥}^2}{n\left({\beta }_0+{\beta }_p{max}_{\mu \in I_i}{\left({\mathit{p}}_{\mu }^{\left(k\right)}\right)}^2+{\beta }_c{max}_{\mu \in I_i}{\left({\mathit{c}}_{\mu }^{\left(k\right)}\right)}^2\right)},i=1,\dots ,n\hfill \end{array}$$ 7: until converge condition (10) 8: return $(\mathit{f},R_0)={\mathit{x}}^{\left(k\right)}$                                                                      ▹ Result $(\mathit{f},R_0)$

3.2. MF-L1 Algorithm

This algorithm employs an $L_1$-norm-based penalty, which is preferred for inducing sparsity in $\mathit{f}$. This approach is based on the assumption that the $f\left(\tau \right)$ distribution is a sparse function characterized by only a few non-zero terms.

The problem of parameter identification is reformulated as the following optimization problem:

$$\begin{array}{cc}\underset{\mathit{x}\ge 0}{min}\hfill & {\displaystyle \left\{\parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+\alpha {\parallel \mathit{x}\parallel }_1\right\}}\hfill \end{array}$$

(11)

where $\alpha >0$ is the regularization parameter computed according to the Balancing Principle (BP) [10]. Following [12], we rewrite (11) as

$$\begin{array}{cc}\underset{\mathit{x}}{min}\hfill & \left\{{\displaystyle \parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+{\alpha \parallel \mathit{x}\parallel }_1+\eta {\parallel \mathit{x}\parallel }_2^2}\right\}\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & \mathit{x}\ge 0\hfill \end{array}$$

(12)

In this new formulation (12), the last $L_2$-based penalty term, ${\eta \parallel \mathit{x}\parallel }_2^2$, is introduced only to ensure that ${\mathit{K}}_e^T{\mathit{K}}_e+\eta \mathit{I}$ is a definite positive matrix to ensure that (12) is well-posed. It is not a regularization term, and a small positive value for $\eta \approx {10}^{-10}$ is fixed.

The MF-L1 algorithm is an iterative procedure where, starting from an initial guess ${\lambda }^{\left(0\right)}$, at each iteration k, an estimate of the parameters $({\mathit{f}}^{\left(k\right)},R_0^{\left(k\right)})$ is computed by solving the parameter estimation problem (12) for fixed ${\alpha }^{\left(k\right)}$ by the truncated Newton interior-point method [17] (see Algorithm 2, step (4)). Then, a new value ${\alpha }^{(k+1)}$ is determined by using the BP (see [10] for a deep theoretical discussion). The BP selects the regularization parameter $\alpha$ so that the data fidelity and the regularization terms are balanced up to a multiplicative factor $\gamma$, i.e.,

$${\gamma \alpha \parallel \mathit{x}\parallel }_1=\parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+\eta {\parallel \mathit{x}\parallel }_2^2.$$

(13)

We fix $\gamma =1$, and we obtain the following rule for the parameter selection:

$$\alpha =\frac{\parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+\eta {\parallel \mathit{x}\parallel }_2^2}{{\parallel \mathit{x}\parallel }_1}.$$

(14)

The MF-L1 method is summarized in Algorithm 2.

Algorithm 2 MF-L1

1: Set $k=0$, $\eta ={10}^{-10}$, and choose a starting guess ${\alpha }^{\left(0\right)}$. 2: repeat 3:     $k=k+1$ 4:     NMRD parameters update. By using the truncated Newton interior-point method, compute $${\mathit{x}}^{\left(k\right)}=\arg \underset{\mathit{x}\ge 0}{min}\parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+{\alpha }^{(k-1)}{\parallel \mathit{x}\parallel }_1+\eta {\parallel \mathit{x}\parallel }_2^2$$ 5:     Regularization parameter update. Set $${\displaystyle {\alpha }^{\left(k\right)}=\frac{\parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}}^{\left(k\right)}{\parallel }_2^2+\eta {\parallel {\mathit{x}}^{\left(k\right)}\parallel }_2^2}{\parallel {\mathit{x}}^{\left(k\right)}{\parallel }_1}}$$ 6: until ${\displaystyle |{\alpha }^{\left(k\right)}-{\alpha }^{(k-1)}|\le \mathrm{Tol}|{\alpha }^{\left(k\right)}|}$ 7: return $(\mathit{f},R_0)={\mathit{x}}^{\left(k\right)}$                                                                      ▹ Result $(\mathit{f},R_0)$

The Matlab implementation of this algorithm is part of the freeware Matlab tool, ModelFreeFFC, which addresses the more complex issue of data exhibiting spurious peaks caused by Quadrupolar Relaxation Enhancement. The software can be downloaded from https://site.unibo.it/softwaredicam/en/modelfree, accessed on 21 May 2024. For a detailed discussion, please refer to [12].

3.3. MF-MUPen Algorithm

This algorithm implements the multi-penalty approach proposed in [13] for the two-dimensional NMR relaxometry data. MF-MUPen solves the following unconstrained minimization problem:

$$\begin{array}{cc}\underset{\mathit{x}}{min}\hfill & {\displaystyle \left\{\parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+\sum _{i=1}^n{\lambda }_i{\left(\mathit{L}\mathit{x}\right)}_i^2+\alpha {\parallel \mathit{x}\parallel }_1\right\}}\hfill \end{array}$$

(15)

which incorporates both penalty functions from MF-UPen and MF-L1. The regularization parameters are then calculated using (7) for ${\lambda }_i,i=1,n$ and (14) for $\alpha$.

To summarize, MF-MUPen is an iterative scheme where, given an initial guess ${\lambda }_i^{\left(0\right)}$, $i=1,\dots ,n$, a parameter estimate $({\mathit{f}}^{\left(k\right)},R_0^{\left(k\right)})$ is computed by solving (15) for fixed ${\lambda }_i^{\left(k\right)}$, $i=1,\dots ,n+1$. Problem (15) is solved by the FISTA method [18], which is one of the most popular methods for minimizing $L_1$-penalized least squares functions. MF-MUPenis stopped when the following condition is satisfied:

$$\sum _{i=1}^n|{\lambda }_i^{(k+1)}-{\lambda }_i^{\left(k\right)}|+|{\alpha }^{(k+1)}-{\alpha }^{\left(k\right)}|<\mathrm{Tol}\left(\sum _{i=1}^n|{\lambda }_i^{\left(k\right)}|+|{\alpha }^{\left(k\right)}|\right).$$

(16)

MF-MUPen is sketched in Algorithm 3.

Algorithms 1–3 require an initial estimate for the regularization parameters. This can be obtained by computing a rough approximation $\tilde{\mathit{x}}$ to the following non-negatively constrained least squares problem:

$$\underset{\mathit{x}\ge 0}{min}{\parallel \mathit{y}-{\mathit{K}}_e\mathit{x}\parallel }_2^2$$

(17)

and then by using the Balancing and Uniform Penalty principles to obtain the initial guess. More precisely,

$${\displaystyle {\alpha }^{\left(0\right)}=\frac{\parallel \mathit{y}-{\mathit{K}}_e\tilde{\mathit{x}}{\parallel }_2^2+\eta {\parallel \tilde{\mathit{x}}\parallel }_2^2}{\parallel \tilde{\mathit{x}}{\parallel }_1}}$$

in Algorithms 2 and 3, and

$$\begin{array}{cc}\hfill {\displaystyle {\lambda }_i^{\left(0\right)}}& =\frac{{∥\mathit{y}-{\mathit{K}}_e\tilde{\mathit{x}}∥}^2}{n\left({\beta }_0+{\beta }_p{max}_{\mu \in I_i}{\left({\tilde{\mathit{p}}}_{\mu }\right)}^2+{\beta }_c{max}_{\mu \in I_i}{\left({\tilde{\mathit{c}}}_{\mu }\right)}^2\right)},i=1,\dots ,n\hfill \end{array}$$

in Algorithms 1 and 3.

Algorithm 3 MF-MUPen

1: Set $k=0$, and choose a starting guess ${\lambda }_i^{\left(0\right)}$, $i=1,\dots ,n$, ${\alpha }^{\left(0\right)}$. 2: repeat 3:     $k=k+1$ 4:     NMRD parameters update. By using the FISTA method, compute $${\mathit{x}}^{\left(k\right)}=\arg \underset{\mathit{x}}{min}\parallel \mathit{y}-{\mathit{K}}_e{\mathit{x}\parallel }_2^2+\sum _{i=1}^n{\lambda }_i^{(k-1)}{\left(\mathit{L}\mathit{x}\right)}_i^2+{\alpha }^{(k-1)}{\parallel \mathit{x}\parallel }_1$$ 5:     Regularization parameter update. Set $$\begin{array}{c}{\displaystyle {\lambda }_i^{\left(k\right)}}{\displaystyle =\frac{{∥\mathit{y}-{\mathit{K}}_e{\mathit{x}}^{\left(k\right)}∥}^2}{n\left({\beta }_0+{\beta }_p{max}_{\mu \in I_i}{\left({\mathit{p}}_{\mu }^{\left(k\right)}\right)}^2+{\beta }_c{max}_{\mu \in I_i}{\left({\mathit{c}}_{\mu }^{\left(k\right)}\right)}^2\right)},i=1,\dots ,n}\hfill \\ {\alpha }^{\left(k\right)}=\frac{{∥\mathit{y}-{\mathit{K}}_e{\mathit{x}}^{\left(k\right)}∥}^2}{n\parallel {\mathit{x}}^{\left(k\right)}{\parallel }_1}\hfill \end{array}$$ 6: until converge condition (16) 7: return $(\mathit{f},R_0)={\mathit{x}}^{\left(k\right)}$                                                                      ▹ Result $(\mathit{f},R_0)$

4. Results and Discussion

In this section, we report and discuss the results obtained by the proposed algorithms on samples of two different materials, which represent typical case tests well.

In Section 4.1, we present the metrics used to evaluate the results’ quality quantitatively and the experimental setting. In Section 4.2, we report and discuss the results obtained by the algorithms.

Numerical computations were carried out using Matlab R2022b on a laptop equipped with an Apple M1 chip with 16 GB of RAM. Please observe that throughout the section, we refer to the frequencies $\nu$ instead of the angular frequencies $\omega$, i.e., $\nu \equiv \omega /\left(2\pi \right)$.

4.1. Experimental Setting

The fitted NMRD profiles, computed by Algorithms 1–3 are compared to the $R_1$ data by means of the ${\chi }^2$ value, which is defined as follows:

$${\chi }^2=\sum _{i=1}^m\frac{{(e_i-y_i)}^2}{m-1}$$

(18)

where $\mathit{e}$ is the estimated data value, i.e.,

$$\mathit{e}=\mathit{K}\tilde{\mathit{f}}+{\tilde{R}}_0$$

with $(\tilde{\mathit{f}},\widetilde{R_0})$ being the computed parameters. We quantitatively compare the computed correlation time distributions $\mathit{f}$ given by the three algorithms, determining the peak values and the area below $\mathit{f}$ in the neighborhood of such peaks, and defining a value such as SpecificWeight. Let us assume that f has $n_p$ local maxima at the correlation times ${\tau }_{c_{\ell }}$, $\ell =1,\dots ,n_p$. Then, we define the neighborhood of interest ${\mathcal{I}}_{\ell }$ using the Full Width at Half Maximum parameter, i.e.,

$${\mathcal{I}}_{\ell }\equiv [{\tau }_{\ell }^{low},{\tau }_{\ell }^{up}]\mathrm{s}.\mathrm{t}.f\left({\tau }_{\ell }^{low}\right)=f\left({\tau }_{\ell }^{up}\right)=\frac{1}{2}f\left({\tau }_{c_{\ell }}\right)\ell =1,\dots ,n_p.$$

We compute the SpecificWeight value for each peak ${\tau }_{c_{\ell }}$, i.e.,

$${\mathrm{SpecificWeight}}_{\ell }=\sum _{j=1}^{n_{\ell }}{\tau }_{c_j}f\left({\tau }_{c_j}\right),{\tau }_{c_j}\in {\mathcal{I}}_{\ell },$$

(19)

where $n_{\ell }$ is the number of correlation times belonging to ${\mathcal{I}}_{\ell }$, $\ell =1,\dots ,n_p$.

The value of the tolerance parameters used in the stopping criteria of all algorithms is $\mathrm{Tol}={10}^{-2}$. Moreover, a maximum number of 10 iterations k has been set but never been reached. The computational cost is evaluated in terms of execution time.

To test the algorithms’ robustness, we apply them to a set of $s=500$ artificial profiles obtained by adding to $R_1$ uniformly distributed noise within the experimental error intervals. With these tests, referred to as dispersion analysis, we aim to explore the intervals containing the recovered parameter $R_0$ and how the calculated estimates scatter around the average value. Additionally, we want to examine how the position and value of the peaks change in the recovered correlation times distributions.

4.2. Numerical Results from FFC Measures

We present the results obtained by applying all the algorithms to NMRD profiles obtained from two experimental samples, Manganese and Poplar, respectively. Both systems are considered a “gold standard” in relaxometry studies, especially when the involvement of paramagnetic species is necessary. The relaxometric properties of aqueous manganese solutions have been thoroughly investigated [19,20] and, as such, these solutions are routinely utilized to assess the performance and stability of instruments. Additionally, the characteristics of Poplar char have been extensively studied [21], making it an effective model for examining the textural properties and functional mobility of solvents within these porous materials. The NMRD profile for the manganese sample was acquired by the authors, while the data pertaining to Poplar char were taken from [22].

These two samples show how the algorithms’ results can complement each other to improve the overall quality of the information provided. We evaluate the global quality of the examined methods in terms of ${\chi }^2$, offset $R_0$, and computation time.

The $R_1$ data for the Manganese Sample are measured at 26 frequency values $\nu$, ranging from ${10}^{-2}$ to ${10}^1$ MHz. The error intervals for these measurements vary from $\pm 0.4$ s⁻¹ to $\pm 1.1$ s⁻¹. These are illustrated by the black error bars in the left panel of Figure 1. The $R_1$ data for the Poplar sample are measured at 21 frequency values $\nu$, ranging from ${10}^{-2}$ to ${10}^1$ MHz. The error intervals for these measurements vary from $\pm 0.06$ s⁻¹ to $\pm 0.3$ s⁻¹. These are illustrated by the black error bar in the right panel of Figure 1.

Table 1. Manganese sample. Computational results of the proposed methods.

Algorithm	${\mathit{R}}_0$ $\left[{\mathbf{s}}^{-\mathbf{1}}\right]$	${\mathit{\chi }}^2$ $[-]$	Computation Time $\left[\mathbf{s}\right]$
MF-UPen	$6.64\times {10}^0$	$6.23\times {10}^{-1}$	$9.86\times {10}^{-1}$
MF-L1	$1.19\times {10}^1$	$5.55\times {10}^{-1}$	$1.35\times {10}^{-1}$
MF-MUPen	$9.98\times {10}^0$	$4.10\times {10}^{-1}$	$1.24\times {10}^0$

presents the estimated parameter $R_0$, the goodness-of-fit measure ${\chi }^2$, and the computation time in seconds obtained by the three algorithms. We observe that MF-L1 achieves a moderate ${\chi }^2$ value and the shortest computation time. In contrast, MF-UPen shows a slightly higher ${\chi }^2$ and takes longer to compute. Conversely, MF-MUPen achieves the best fit, as evidenced by the lowest ${\chi }^2$. This suggests superior model accuracy, with a modest increase in computation time.

Table 2. Poplar sample. Computational results of the proposed methods.

Algorithm	${\mathit{R}}_0$ $\left[{\mathbf{s}}^{-\mathbf{1}}\right]$	${\mathit{\chi }}^2$ $[-]$	Computation Time $\left[\mathbf{s}\right]$
MF-UPen	$5.40\times {10}^0$	$7.94\times {10}^{-3}$	$7.41\times {10}^{-2}$
MF-L1	$5.41\times {10}^0$	$8.84\times {10}^{-3}$	$6.65\times {10}^{-2}$
MF-MUPen	$5.41\times {10}^0$	$2.19\times {10}^{-2}$	$3.76\times {10}^{-1}$

outlines the computational results obtained for the Poplar sample by the three algorithms. MF-UPen and MF-L1 both report nearly identical values for $R_0$, with minimal ${\chi }^2$ and very short computation times, indicating efficient and effective performance. However, MF-MUPen, while yielding a similar $R_0$ to the other two algorithms, shows a higher ${\chi }^2$ value, suggesting a slightly poorer fit. Additionally, MF-MUPen requires longer computation time.

The outcomes for the Manganese and Poplar samples represent two scenarios, each indicative of the potential variability in sample analysis. This diversity highlights the importance of utilizing multiple methods to fully understand sample characteristics under varying conditions. The peak analysis for both the Manganese and Poplar samples across the three methods is performed by plotting the correlation times amplitudes $\mathit{f}$ computed by each method in Figure 2 and reporting the peak position amplitudes, half-width, and SpecificWeights for each sample in

Table 3. Manganese sample analysis. Position (${\tau }_c$) and amplitude $f\left({\tau }_c\right)$ of the distribution peaks sorted by $f\left({\tau }_c\right)$.

Algorithm	${\mathit{\tau }}_{\mathit{c}}$ $\left[\mathbf{\mu }\mathbf{s}\right]$	$\mathit{f}\left({\mathit{\tau }}_{\mathit{c}}\right)$ $[\mathbf{a}.\mathbf{u}.]$	Half-Width $\left[\mathbf{\mu }\mathbf{s}\right]$	SpecificWeight $[\mathbf{a}.\mathbf{u}.]$
MF-UPen	$8.30\times {10}^{-3}$	$9.78\times {10}^1$	$8.19\times {10}^{-3}$	$6.94\times {10}^0$
	$3.76\times {10}^{-2}$	$2.95\times {10}^1$	$8.47\times {10}^{-3}$	$2.13\times {10}^0$
	$\mathbf{7}.\mathbf{74}\times {\mathbf{10}}^{-\mathbf{1}}$	$9.61\times {10}^0$	$1.18\times {10}^{-1}$	$1.09\times {10}^1$
	$4.86\times {10}^{-1}$	$4.63\times {10}^0$	$5.67\times {10}^{-2}$	$2.25\times {10}^0$
	$3.43\times {10}^{-1}$	$4.85\times {10}^{-2}$	$4.98\times {10}^{-2}$	$2.17\times {10}^{-2}$
MF-L1	$1.05\times {10}^{-2}$	$5.13\times {10}^2$	$1.34\times {10}^{-3}$	$6.30\times {10}^0$
	$4.23\times {10}^{-2}$	$3.63\times {10}^1$	$5.36\times {10}^{-3}$	$1.75\times {10}^0$
	$\mathbf{7}.\mathbf{74}\times {\mathbf{10}}^{-\mathbf{1}}$	$1.16\times {10}^1$	$1.21\times {10}^{-1}$	$1.24\times {10}^1$
	$3.43\times {10}^{-1}$	$2.28\times {10}^0$	$4.00\times {10}^{-2}$	$7.84\times {10}^{-1}$
	$1.63\times {10}^{-3}$	$3.79\times {10}^{-3}$	$7.04\times {10}^{-4}$	$5.17\times {10}^{-5}$
MF-MUPen	$9.33\times {10}^{-3}$	$8.49\times {10}^1$	$1.09\times {10}^{-2}$	$6.98\times {10}^0$
	$4.75\times {10}^{-2}$	$1.09\times {10}^1$	$1.55\times {10}^{-2}$	$1.47\times {10}^0$
	$\mathbf{7}.\mathbf{74}\times {\mathbf{10}}^{-\mathbf{1}}$	$1.15\times {10}^1$	$1.23\times {10}^{-1}$	$1.24\times {10}^1$
	$3.43\times {10}^{-1}$	$1.49\times {10}^0$	$5.24\times {10}^{-2}$	$7.56\times {10}^{-1}$

and

Table 4. Poplar sample analysis. Position (${\tau }_c$) and amplitude $f\left({\tau }_c\right)$ of the distribution peaks sorted by $f\left({\tau }_c\right)$.

Algorithm	${\mathit{\tau }}_{\mathit{c}}$ $\left[\mathbf{\mu }\mathbf{s}\right]$	$\mathit{f}\left({\mathit{\tau }}_{\mathit{c}}\right)$ $[\mathbf{a}.\mathbf{u}.]$	Half-Width $\left[\mathbf{\mu }\mathbf{s}\right]$	SpecificWeight $[\mathbf{a}.\mathbf{u}.]$
MF-UPen	$4.23\times {10}^{-2}$	$1.02\times {10}^0$	$3.92\times {10}^{-2}$	$3.54\times {10}^{-1}$
	$3.05\times {10}^{-1}$	$4.65\times {10}^{-1}$	$1.25\times {10}^{-1}$	$5.00\times {10}^{-1}$
	$1.56\times {10}^0$	$1.75\times {10}^{-1}$	$5.38\times {10}^{-1}$	$8.05\times {10}^{-1}$
	$3.51\times {10}^0$	$5.53\times {10}^{-2}$	$9.32\times {10}^{-1}$	$4.53\times {10}^{-1}$
MF-L1	$4.75\times {10}^{-2}$	$3.77\times {10}^0$	$1.04\times {10}^{-2}$	$3.34\times {10}^{-1}$
	$2.72\times {10}^{-1}$	$1.59\times {10}^0$	$3.50\times {10}^{-2}$	$5.11\times {10}^{-1}$
	$1.56\times {10}^0$	$4.85\times {10}^{-1}$	$1.89\times {10}^{-1}$	$8.20\times {10}^{-1}$
	$3.94\times {10}^0$	$1.15\times {10}^{-1}$	$4.60\times {10}^{-1}$	$4.55\times {10}^{-1}$
MF-MUPen	$4.23\times {10}^{-2}$	$1.87\times {10}^0$	$2.13\times {10}^{-2}$	$3.40\times {10}^{-1}$
	$2.72\times {10}^{-1}$	$7.80\times {10}^{-1}$	$7.49\times {10}^{-2}$	$5.11\times {10}^{-1}$
	$1.56\times {10}^0$	$4.41\times {10}^{-1}$	$2.02\times {10}^{-1}$	$8.19\times {10}^{-1}$
	$3.94\times {10}^0$	$1.14\times {10}^{-1}$	$4.60\times {10}^{-1}$	$4.50\times {10}^{-1}$

.

In Table 3, we observe a perfect agreement among the three methods in locating the peak at the longest correlation time $\tau =7.74\times {10}^{-1}\left[\mathsf{\mu }\mathrm{s}\right]$. Meanwhile, MF-UPen and MF-L1 have quite good agreement at the intermediate times: ${\tau }_c=3.76\times {10}^{-2}\left[\mathsf{\mu }\mathrm{s}\right]$ and ${\tau }_c=4.23\times {10}^{-2}\left[\mathsf{\mu }\mathrm{s}\right]$, respectively. The distribution pattern in Figure 2, left panel, shows similarity features between MF-UPen and MF-MUPen, and reveals a tendency of MF-L1 to resolve multiple components at the shortest times.

In the case of the Poplar sample, as shown in the right panel of Figure 2 and Table 4, we observe a tighter clustering of peaks across the methods, particularly at the highest amplitude peak around ${\tau }_c=4.23\times {10}^{-2}\left[\mathsf{\mu }\mathrm{s}\right]$. This suggests that all three methods are in agreement concerning the main features of the Poplar sample’s distribution.

4.3. Dispersion Analysis

The robustness of the methods is investigated through the dispersion analysis, described in Section 4.1. The boxplots in Figure 3 offer a comparative view of algorithmic performance on the two samples. Each boxplot outlines the algorithms’ interquartile range (IQR) and median of ${\chi }^2$ values.

We observe uniformity in medians for the Manganese sample, with outliers highlighted by red plus symbols, suggesting occasional significant deviations for MF-UPen. The symmetry of the data is evident from the whiskers’ lengths.

Conversely, the Poplar sample exhibits a tighter IQR for each algorithm, denoting less variability. Despite the close median values indicating consistent algorithmic performance, outliers for MF-L1 reveal notable deviations in some cases.

Table 5. Comparison of R0 confidence intervals, mean R0, and median for Manganese and Poplar samples.

Sample	Algorithm	R0 Confidence Interval $\left[{\mathbf{s}}^{-\mathbf{1}}\right]$	R0 Mean $\left[{\mathbf{s}}^{-\mathbf{1}}\right]$	Median $\left[{\mathbf{s}}^{-\mathbf{1}}\right]$
Manganese	MF-UPen	[5.240, 9.253]	$7.25\times {10}^0$	$8.26\times {10}^0$
	MF-L1	[9.251, 12.12]	$1.07\times {10}^1$	$1.12\times {10}^1$
	MF-MUPen	[9.652, 11.38]	$1.05\times {10}^1$	$1.05\times {10}^1$
Poplar	MF-UPen	[5.363, 5.406]	$5.39\times {10}^0$	$5.39\times {10}^0$
	MF-L1	[5.370, 5.413]	$5.39\times {10}^0$	$5.39\times {10}^0$
	MF-MUPen	[5.370, 5.416]	$5.39\times {10}^0$	$5.39\times {10}^0$

compares the $R_0$ confidence intervals [23], mean $R_0$, and medians for both Manganese and Poplar samples across the three algorithms.

The confidence intervals and mean $R_0$ values suggest a wider range of estimates for the Manganese sample, indicating a less uniform agreement among the algorithms. The median values, while closer, still reflect notable variation between the algorithms, suggesting that the model fit depends on the algorithm applied.

Conversely, the Poplar sample demonstrates remarkable consistency, with both confidence intervals and mean $R_0$ values being nearly identical across all three algorithms. The median values also closely align, reinforcing the observation of uniform performance. This indicates that for the Poplar sample, the choice of algorithm does not significantly influence the outcome, and all three algorithms provide equivalent information.

Table 5 represents two distinct scenarios that may emerge when these algorithms are applied to samples with varying characteristics. In the case of the Poplar sample, the outcome from all three algorithms is congruent, implying that the algorithms are robust and interchangeable for this type of sample. Conversely, the Manganese sample demonstrates less consistency across the algorithms, suggesting that additional insights from alternative investigative methods are necessary to supplement the analysis.

Concerning the distribution intensities, we computed the mean distribution obtained by each method and analyzed the peak positions and amplitudes analogously to Table 1 and Table 4 for the single samples.

From

Table 6. Manganese sample. Analysis of mean distribution over 500 realizations. Position (${\tau }_c$) and amplitude $f\left({\tau }_c\right)$ of the distribution peaks sorted by $f\left({\tau }_c\right)$.

Algorithm	${\mathit{\tau }}_{\mathit{c}}$ [μs]	$\mathit{f}\left({\mathit{\tau }}_{\mathit{c}}\right)$ $[\mathbf{a}.\mathbf{u}.]$	Half-Width [μs]	SpecificWeight $[\mathbf{a}.\mathbf{u}.]$
MF-UPen	$9.326\times {10}^{-3}$	$8.640\times {10}^1$	$8.889\times {10}^{-3}$	$6.949\times {10}^0$
	$3.765\times {10}^{-2}$	$1.727\times {10}^1$	$9.075\times {10}^{-3}$	$1.617\times {10}^0$
	$\mathbf{7}.\mathbf{743}\times {\mathbf{10}}^{-\mathbf{1}}$	$9.004\times {10}^0$	$1.385\times {10}^{-1}$	$1.168\times {10}^1$
	$5.462\times {10}^{-1}$	$9.534\times {10}^{-1}$	$1.139\times {10}^{-1}$	$4.464\times {10}^0$
	$5.995\times {10}^{-2}$	$7.407\times {10}^{-1}$	$7.643\times {10}^{-3}$	$1.208\times {10}^{-1}$
MF-L1	$1.048\times {10}^{-2}$	$1.642\times {10}^2$	$4.173\times {10}^{-3}$	$6.216\times {10}^0$
	$6.579\times {10}^{-3}$	$3.024\times {10}^1$	$4.493\times {10}^{-4}$	$5.120\times {10}^{-1}$
	$\mathbf{7}.\mathbf{743}\times {\mathbf{10}}^{-\mathbf{1}}$	$1.061\times {10}^1$	$1.264\times {10}^{-1}$	$1.215\times {10}^1$
	$3.765\times {10}^{-2}$	$8.763\times {10}^0$	$1.897\times {10}^{-2}$	$1.702\times {10}^0$
	$3.854\times {10}^{-1}$	$5.470\times {10}^{-1}$	$7.072\times {10}^{-2}$	$5.642\times {10}^{-1}$
	$4.863\times {10}^{-1}$	$3.761\times {10}^{-1}$	$3.575\times {10}^{-2}$	$4.885\times {10}^{-1}$
	$2.420\times {10}^{-1}$	$1.233\times {10}^{-1}$	$1.700\times {10}^{-2}$	$7.589\times {10}^{-2}$
MF-MUPen	$1.048\times {10}^{-2}$	$8.461\times {10}^1$	$9.725\times {10}^{-3}$	$6.954\times {10}^0$
	$\mathbf{7}.\mathbf{743}\times {\mathbf{10}}^{-\mathbf{1}}$	$9.675\times {10}^0$	$1.466\times {10}^{-1}$	$1.219\times {10}^1$
	$4.229\times {10}^{-2}$	$5.630\times {10}^0$	$2.526\times {10}^{-2}$	$2.163\times {10}^0$
	$4.863\times {10}^{-1}$	$4.064\times {10}^{-1}$	$1.402\times {10}^{-1}$	$9.117\times {10}^{-1}$

and

Table 7. Poplar sample. Analysis of mean distribution over 500 realizations. Position (${\tau }_c$) and amplitude $f\left({\tau }_c\right)$ of the distribution peaks sorted by $f\left({\tau }_c\right)$.

Algorithm	${\mathit{\tau }}_{\mathit{c}}$ [μs]	$\mathit{f}\left({\mathit{\tau }}_{\mathit{c}}\right)$ $[\mathbf{a}.\mathbf{u}.]$	Half-Width [μs]	SpecificWeight $[\mathbf{a}.\mathbf{u}.]$
MF-UPen	$\mathbf{4}.\mathbf{229}\times {\mathbf{10}}^{-\mathbf{2}}$	$7.234\times {10}^{-1}$	$5.094\times {10}^{-2}$	$3.369\times {10}^{-1}$
	$2.719\times {10}^{-1}$	$1.830\times {10}^{-1}$	$8.540\times {10}^{-2}$	$2.547\times {10}^{-1}$
	$1.748\times {10}^0$	$7.503\times {10}^{-2}$	$1.211\times {10}^0$	$1.091\times {10}^0$
MF-L1	$\mathbf{4}.\mathbf{229}\times {\mathbf{10}}^{-\mathbf{2}}$	$1.652\times {10}^0$	$1.794\times {10}^{-2}$	$2.811\times {10}^{-1}$
	$2.719\times {10}^{-1}$	$2.662\times {10}^{-1}$	$1.465\times {10}^{-1}$	$4.595\times {10}^{-1}$
	$8.498\times {10}^{-2}$	$1.744\times {10}^{-1}$	$2.187\times {10}^{-2}$	$6.473\times {10}^{-2}$
	$1.204\times {10}^{-1}$	$1.287\times {10}^{-1}$	$1.081\times {10}^{-2}$	$3.840\times {10}^{-2}$
	$1.963\times {10}^0$	$8.927\times {10}^{-2}$	$1.045\times {10}^0$	$1.007\times {10}^0$
	$1.556\times {10}^0$	$8.810\times {10}^{-2}$	$1.176\times {10}^{-1}$	$3.459\times {10}^{-1}$
	$6.136\times {10}^{-1}$	$2.873\times {10}^{-2}$	$7.136\times {10}^{-2}$	$4.805\times {10}^{-2}$
	$7.925\times {10}^0$	$1.298\times {10}^{-3}$	$5.589\times {10}^{-1}$	$2.304\times {10}^{-2}$
MF-MUPen	$\mathbf{4}.\mathbf{229}\times {\mathbf{10}}^{-\mathbf{2}}$	$1.312\times {10}^0$	$2.518\times {10}^{-2}$	$3.009\times {10}^{-1}$
	$2.420\times {10}^{-1}$	$2.697\times {10}^{-1}$	$1.420\times {10}^{-1}$	$4.568\times {10}^{-1}$
	$1.072\times {10}^{-1}$	$1.047\times {10}^{-1}$	$2.110\times {10}^{-2}$	$5.096\times {10}^{-2}$
	$1.556\times {10}^0$	$7.936\times {10}^{-2}$	$1.204\times {10}^0$	$1.180\times {10}^0$
	$2.205\times {10}^0$	$7.194\times {10}^{-2}$	$1.217\times {10}^{-1}$	$2.996\times {10}^{-1}$
	$6.893\times {10}^{-1}$	$2.774\times {10}^{-2}$	$1.504\times {10}^{-1}$	$8.125\times {10}^{-2}$
	$4.431\times {10}^0$	$7.331\times {10}^{-3}$	$2.851\times {10}^{-1}$	$7.834\times {10}^{-2}$

, we notice that MF-L1 tends to identify a greater number of peaks compared to the other two methods, suggesting a higher sensitivity of the algorithm.

In the case of the Manganese sample, the data reported in Table 6 show that there is a perfect correspondence in peak position at the longest correlation time ${\tau }_c=7.743\times {10}^{-1}$ among the three algorithms, while the peaks at shortest and intermediate times are split into multiple components.

Concerning the Poplar sample (Table 7), we observe that all algorithms show identical peak positioning corresponding to the largest amplitude, which was reached at the shortest correlation time ${\tau }_c=4.229\times {10}^{-2}\left[\mathsf{\mu }\mathrm{s}\right]$. At the longest correlation times, MF-UPen finds a single peak around ${\tau }_c=1.748\times {10}^0\left[\mathsf{\mu }\mathrm{s}\right]$, while MF-L1 and MF-MUPen split the amplitudes into two peaks at ${\tau }_c=1.556\times {10}^0,1.963\times {10}^0\left[\mathsf{\mu }\mathrm{s}\right]$ and ${\tau }_c=1.556\times {10}^0,2.205\times {10}^0\left[\mathsf{\mu }\mathrm{s}\right]$, respectively.

However, despite the differences in the number of peaks identified, Figure 4 shows that all three algorithms exhibit a fundamental robustness in the localization of the positions of the highest peaks.

From Table 5, we observe that MF-MUPen has the smallest confidence intervals in both samples and Figure 3 shows that the number of outliers in MF-MUPen is smaller than in the other methods. This, together with Figure 4, suggests a higher robustness of MF-MUPen compared to the other methods.

5. Conclusions

This study reviews three algorithms designed for the analysis of NMRD profiles from experimental data, showing their efficacy in two different datasets representative of possible application scenarios. The comparative analysis of the three algorithms, while revealing a general consistency in identifying the primary peak positions, shows that MF-MUPen exhibits higher robustness in the presence of data noise.

In conclusion, the proposed algorithms have the potential to significantly improve the quality of NMRD profile analysis and enrich the toolkit available to researchers for investigating the molecular dynamics of various samples.