Abstract

Of late, we have put forward a new branch called high-order derivative signal processing. This investigative strategy is universally relevant for all spectroscopies, the progress of which ultimately depends on resolution improvement and noise suppression. The high-order, non-parametric derivative fast Padé transform (dFPT) simultaneously solves these two problems of utmost importance. The present work goes one critical step further than our previous two studies on this particular topic, by setting up the goal of validating the non-parametric dFPT by its parametric counterpart. This is done by comparing the full lineshapes of derivative envelopes from the non-parametric dFPT with the corresponding derivative component spectra from the parametric dFPT. The non-parametric dFPT, as a shape estimator, never solves the quantification problem (or equivalently, the spectral analysis problem via e.g. an eigenvalue problem, rooting characteristic/secular equations, etc.). The parametric dFPT first solves the quantification problem from which the lineshapes of components and envelopes are plotted. Thus, if the derivative component spectra from the parametric dFPT could be fully reconstructed by the derivative envelopes from the non-parametric dFPT, the goal of achieving quantification would be done by derivative lineshape processing alone. This would amount to providing stand-alone quantification without actually solving the quantification problem (and, of course, without fitting, either). The present study accomplishes this goal, with an important application to data encountered in magnetic resonance spectroscopy for clinical diagnostics of breast cancer.

1 Introduction

We pursue a novel pathway in analyzing spectra by non-parametric high-order derivative signal processing. Hereafter, when dealing with spectra, derivative operator \( {\text{D}}_{\nu }^{m} \) of order m is taken with respect to the sweep frequency v as \( {\text{D}}_{\nu }^{m} = \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \). This is a general methodology for all spectroscopies, including those based on magnetic resonance (MR), which is the focus of the present work within the non-parametric derivative fast Padé transform (dFPT). In our two most recent studies [1, 2], the dFPT was shown to be capable of simultaneously improving resolution and suppressing noise. By contrast, the derivative fast Fourier transform (dFFT) amplifies noise with increased derivative order m [1]. As demonstrated in Refs. [1, 2], using the phase-insensitive magnitude mode, the derivatives of envelopes in the non-parametric dFPT systematically reduce the peak widths and enhance the peak heights. Crucially, these key quantities have been shown [2] to be uniquely related to the corresponding absorptive peak parameters of the non-derivative (m = 0) fast Padé transform (FPT) [2].

To emphasize, in both Refs. [1, 2], the explicit numerical computations were performed by exclusively employing the non-parametric dFPT. The present investigation extends these previous examinations to the parametric dFPT. The goal is to compare the envelopes (from the non-parametric dFPT) and the components (from the parametric dFPT), so as to quantitatively validate the former. Of course, viewed on its own, the non-derivative parametric FPT does not need to resort to the dFPT. This occurs because in the non-derivative parametric FPT, the peak parameters are already available by solving the quantification problem. Nevertheless, if the component lineshapes are computed by the parametric dFPT, the obtained results, as the “gold standard”, would retroactively confirm the entire spectral lineshape profiles reconstructed as the envelopes by the non-parametric dFPT. This is what we have set to accomplish in the present study. This undertaking is necessary in order to validate the full quantitative reliability of the non-parametric dFPT which is, at the onset, exclusively a shape estimator.

As to the application, this entire context is herein focused on magnetic resonance spectroscopy (MRS) for breast cancer diagnostics. In this particular problem of utmost clinical relevance, spectral envelopes abound with closely-overlapping resonances. Therein lie a number of diagnostically informative metabolites, among which of key importance is phosphocholine (PC), as a recognized breast cancer biomarker [3]. The PC peak is entirely invisible in the non-derivative envelope (computed by any processor including the non-derivative FPT), due to the presence of a much taller resonance, that of phosphoethanolamine (PE), which is not a cancer biomarker.

The proof of reliable quantification of PC and reconstruction of its component lineshape with the correct concentration by the non-parametric dFPT would be provided by full agreement with the parametric dFPT. This would be the most stringent benchmarking of derivative magnetic resonance spectroscopy (dMRS) for diagnostics within clinical oncology. Importantly, PC is also a recognized biomarker of other cancers [4, 5]. Of course, the overall power of derivative signal processing is not limited to this pair of strongly coupled resonances. Rather, it also applies to any number of tightly overlapping peaks across the entire spectrum under study.

Although our applications of the dFPT are currently within biomedicine, the design of derivative signal processing has much wider repercussions. The reason is that the present concept of processing generic data by using high-order derivative estimations could be directly of use in all spectroscopies (atomic, molecular, nuclear, sub-nuclear) with several specialized branches, e.g. photoelectron spectroscopy (PES) [6], infrared spectroscopy (IRS) [7, 8] or Raman spectroscopy [9]. This could also be said for other areas independent of spectroscopy and signal processing, i.e. whenever the unknown components or constituents are sought from the given chemical mixture or any other compound material, as encountered in many fields (powder diffraction phenomena [10], etc.). For example, from the viewpoint of signal processing alone, the time signals encoded by ion-cyclotron resonance mass spectrometry (ICRMS) [11, 12, 13] or time of flight mass spectrometry (TOFMS) [14], or IRS [7, 8] are all amenable to exactly the same type of estimation done presently within MRS, which outside medicine is called nuclear magnetic resonance (NMR) spectroscopy [15, 16].

2 Theory

We will separately give the outlines of the main ingredients of the two versions of the derivative fast Padé transform, dFPT, with the non-parametric and parametric aspects. Prior to this stratified presentation, we shall summarize the customary, non-derivative fast Padé transform, FPT, by giving its principal features.

where \( \omega \) is the angular (cyclic) frequency (given in cycles per second) and \( \omega = 2\pi \nu \) with \( \nu \) being the linear sweep frequency in hertz (Hz). The spectrum \( G_{N} \left( {z^{ - 1} } \right) \) is the MacLaurin polynomial, or the so-called \( z \)-transform [17]. This is a truncated version of the MacLaurin series:

2.2 The output spectra

2.2.1 The non-derivative, non-parametric FPT

When \( N \) is finite \( \left( {N < \infty } \right), \) as is always the case in practice (in any measurement or in any computer simulation), the spectrum (2.1), despite being exact, as the finite-rank Green function, is nevertheless ill-adapted to analyzing time signals from all spectroscopies. The reason is in the fact that possible singularities in the analyzed function are not captured by the polynomial \( G_{N} \) from (2.1). On the other hand, the poles as a type of singularity are omnipresent in general spectra. One of the rescues from this obstacle is in the realization that the ratio of two polynomials \( P/Q \) (i.e. a rational polynomial) is the mathematical model of choice to adequately approximate all functions with polar singularities that are, by definition, the roots of the denominator polynomial, \( Q = 0. \) Such an observation naturally leads to the introduction of the following two-fold form of the non-derivative, non-parametric FPT via:

2.2.2 The non-derivative, parametric FPT

Overall, the non-parametric \( {\text{FPT}}^{\left( \pm \right)} \) can provide only the envelope lineshapes, or equivalently, the profiles of total shape spectra. This could suffice if such spectra were to contain solely well-isolated, non-overlapping peaks. However, such an idealization is hardly ever encountered in practice when dealing with spectra generated from encoded time signals. Quite the contrary, real-life spectra abound with overlapped resonances for which no customary (non-derivative) non-parametric shape estimator could complete the main task of all spectroscopies: quantification of the physical components of envelopes. One way out of this obstacle, while still working within the realm of the conventional (non-derivative) estimations, is to resort to parametric signal processing. Here, the word “parametric” implies a suitable and physically justified parametrization of a complicated phenomenon under study by a set of quantities that properly describe the major aspects of the problem. In the case of spectra, these quantities are the peak parameters (peak positions, widths, heights and phases). The object of quantitative spectroscopies is to faithfully reconstruct these four real-valued parameters for each physical resonance in the spectrum. The parametric \( {\text{FPT}}^{\left( \pm \right)} \) accomplishes this task by solving the quantification problem consisting of two basic steps: rooting the denominator polynomials \( Q_{K}^{ \pm } \) and finding the Cauchy residues of the rational polynomials \( P_{K}^{ \pm } /Q_{K}^{ \pm } \). The zeros \( z_{k}^{ \pm } \) of \( Q_{K}^{ \pm } \left( {z^{ \pm 1} } \right) \) give the complex fundamental frequencies \( \omega_{k}^{ \pm } = \mp \left( {i/\tau } \right){ \ln }\left( {z_{k}^{ \pm } } \right) \):

In practice, we avoid polynomial rooting. Instead, this non-linear operation is replaced by the linear eigenvalue problem of the corresponding Hessenberg (companion) matrix [17, 18]. This matrix is extremely sparse with its elements being the polynomial expansion coefficients on the 1st row, unity on the main diagonal and zero elsewhere. Such an extraordinary sparseness of the Hessenberg matrix permits work with huge dimensions. The eigenvalues of the Hessenberg matrix are exactly equal to the roots of the associated polynomial.

Only the simple roots \( \left\{ {z_{k}^{ \pm } } \right\} \) are present in (2.12). The amplitudes \( d_{k}^{ \pm } \) take the more general forms for degenerate spectra that contain one or more confluent fundamental frequencies. In such spectra, more than one fundamental amplitude correspond to the same fundamental frequency. This is modeled by including the multiple (or repeated) roots of the polynomials \( Q_{K}^{ \pm } \). In such a case, the amplitudes become:

where \( m_{k} \) is the multiplicity of the kth root (for brevity, we did not put the superscript on \( m_{k} \)). The essential difference between a non-degenerate and degenerate spectrum is that the former and the latter are complex Lorentzian and non-Lorentzian lineshapes, respectively. The fact that the \( {\text{FPT}}^{\left( \pm \right)} \) can handle non-Lorentzian lineshapes makes the Padé-based processing superior to the other existing parametric estimators that are limited to Lorentzians alone.

After completing the spectral analysis (i.e. after solving the outlined quantification problem), the parametric \( {\text{FPT}}^{\left( \pm \right)} \) can construct the spectra in different forms (canonical, partial fractions, etc.) In practice, the most frequently used is the representation given by the Heaviside partial fraction decomposition. For non-degenerate spectra, this representation reads as:

This is one of the ways in which the FPT achieves its cross-validation by verifying the overall concordance between its two variants, the \( {\text{FPT}}^{\left( + \right)} \) and \( {\text{FPT}}^{\left( - \right)} \).

For noise suppression, Padé-based signal processing uses the so-called signal–noise separation (SNS) concept. The SNS procedure exploits the pole-zero coincidence, which is associated with the emergence of Froissart doublets in spectra. Random noise is inherently unstable. This means that even the slightest perturbation (e.g. varying levels of truncation of the total signal length \( N \)) can significantly change the distribution of the noisy part of the encoded FIDs. This is replicated in the confluence or near-confluence of poles and zeros in the given Padé spectrum. Since a polynomial quotient is a meromorphic function, its zeros and poles are given by the roots of the numerator and denominator polynomials, respectively. The functions whose only singularities are poles are called meromorphic functions. Thus, to implement the SNS concept, we also need to find the roots of the numerator polynomials. If the poles and zeros coincide exactly for certain resonances, their amplitudes (and, consequently, peak heights) will all precisely be equal to zero. Such parts of the polynomial quotients cancel each other. This pole-zero cancellation is most obvious in the canonical representation of the Padé spectrum [19]. It is also obvious in the Heaviside partial fractions (2.15) where \( d_{k}^{ \pm } = 0 \) whenever \( z_{k}^{ \pm } \) satisfies both \( Q_{K}^{ \pm } (z_{k}^{ \pm } ) = 0 \) and \( P_{K}^{ \pm } (z_{k}^{ \pm } ) = 0, \) which occurs for pole-zero coincidences. In the other case, when there are approximate pole-zero confluences, the associated resonances have negligibly small amplitudes. Any perturbation would alter the parameters (position, width, height, phase) of such spectral structures. Being weak and markedly unstable, these resonances share the like characteristics of noise. They are binned as spurious (unphysical, ghost, extraneous) and are viewed as noise resonances (or Froissart doublets because of their association with pole-zero pairs). On the other hand, the resonances with no pole-zero confluences exhibit robust stability in face of perturbations. These stable resonances are categorized as genuine (physical) and are, therefore, associated with the true content of the spectra and the corresponding signals. The genuine and spurious resonances are then retained and discarded, respectively, from the final output list (linelist) in the parametric \( {\text{FPT}}^{\left( \pm \right)} \). This is the essence of the SNS concept, by which the \( {\text{FPT}}^{\left( \pm \right)} \) succeed in automatically identifying and reducing noise. Such an improvement in signal–noise ratio (SNR) by the Padé-based estimation is to be credited not only to the SNS concept, but also to the non-linear mathematical form of the polynomial quotients. Such quotients reduce noise both in the non-parametric and parametric \( {\text{FPT}}^{\left( \pm \right)} \). This is sharply opposed to the FFT, which as a linear processor transfers the entire unaltered noise from the time to the frequency domain:

Strictly speaking, this expression is the discrete Fourier transform (DFT), which we formally call the FFT in (2.19). As is well-known, the FFT and the associated DFT have a very significant computational difference. This difference is in some \( N{ \log }_{2} N \) (fast) and \( N^{2} \left( {\text{slow}} \right) \) multiplicative complexity for the FFT and DFT, respectively, when \( N \) is set to be a composite integer of the form \( N = 2^{m} \left( {m = 0,1,2, \ldots } \right) \) [17, 18].

2.3 Physical interpretation of the output linelist of peak parameters

Overall, the FPT (both parametric and non-parametric) is firmly grounded in the best theory of approximations by means of rational polynomials. According to the generalized Weierstrass theorem, any function (rational, irrational, with or without singularities, continuous or discontinuous) can be approximated to within any prescribed level of accuracy, by suitable rational polynomials, provided that the generated numerator and denominator polynomials are of sufficiently high degree. Even singular functions with branch points and branch cuts can be successfully described by a sequence of chained Froissart doublets [19] from rational polynomials. What is more; the parametric version of the FPT goes beyond the envelope of the analyzed spectrum by peering into its inner structure and uncovering the hidden components. The physical information which is unfolded by the parametric FPT is entirely contained in the output linelist with the peak-by-peak signatures associated with the input time signal \( c_{n} \), which is modeled as:

This is a geometric progression showing that the factor \( c_{n,k} \), as the kth component of the nth time signal point \( c_{n} \), is built from the nth power of the complex harmonic function \( z_{k}^{{}} \) whose strength (complex intensity) is equal to \( d_{k} . \) Each attribute of this peak signature tells the physics story (on the level of MR) about the resonances. Thus, the real part, \( {\text{Re}}\left( {\nu_{k} } \right), \) of the kth fundamental frequency \( \nu_{k} \) informs about the chemical shift (i.e. the frequency at which the given molecule resonates with the externally-applied radiofrequency (RF) pulse, in concert with the gradient and static magnetic fields). This information, in turn, tells us about the extent of shielding of nuclei by the surrounding electronic cloud within the parent molecules. The more electrons in the atomic and molecular environment of MR sensitive nuclei, the more shielding, the more local weakening of the external static magnetic field strength \( B_{0} \). This is the essence of the definition of the resonant frequency, which is not the Larmor frequency \( \nu_{\text{L}} \). Rather, it is a miniscule departure from \( \nu_{\text{L}} \) caused by the mentioned electronic shielding, which locally alters \( B_{0} \). Similarly, the reciprocal of the imaginary part, \( {\text{Im}}\left( {\nu_{k} } \right), \) of the kth fundamental frequency speaks about the time rate constant \( T_{2k}^{*} \) of decay of the resonant unstable state of the nucleus. This nucleus, after undergoing the external perturbation, is first excited to a transient state, and then recovered or relaxed through de-excitation to an energetically lower state.

The linewidth itself \( 2\pi {\text{Im}}\left( {\nu_{k} } \right) \) is interpretationally connected to the spin–spin relaxation time \( T_{2k}^{*} . \) Namely, transient signals that decay slowly (long \( T_{2k}^{*} \)) generate sharp, narrow peaks with smaller \( 2\pi {\text{Im}}\left( {\nu_{k} } \right) \). Conversely, signals with faster decays (short \( T_{2k}^{*} \)) correspond to broad peaks with larger \( 2\pi {\text{Im}}\left( {\nu_{k} } \right) \) in the frequency spectrum. Further, the magnitude \( \left| {d_{k} } \right| \), as the third parameter of the kth resonance, carries key information. This is the doubled peak area for a purely absorptive Lorentzian lineshape resulting from the use of \( c_{n} \) from (2.20). This occurs because the peak area \( a_{k} \) of a non-derivative Lorentzian absorption is proportional to the product of the peak height and peak width. On the other hand, the peak height is proportional to the ratio of \( \left| {d_{k} } \right| \) and the peak width. Therefore, the peak width cancels out from the peak area which, in turn, becomes proportional to \( \left| {d_{k} } \right| \) alone in the setting of the absorptive Lorentzians, so that \( a_{k} = \left| {d_{k} } \right|/2 \). Returning to (2.21), we see that the kth signal component \( c_{n,k} \) is reduced to \( d_{k} \) at the time \( t = 0 \) (i.e. at \( n\tau = 0 \)):

Thus, in the context of (2.21), the halved magnitude \( \left| {d_{k} } \right|/2 \) of the \( k{\text{th}} \) component \( \left| c_{0,k} \right| \) of \( c_{n} \) is the discussed kth peak area \( a_{k} \). This is the reason for the traditionally held view that the NMR amplitude of the signal in the time domain (at the onset of data acquisition, \( t = 0 \)) can be referred to as the doubled peak area of the absorptive Lorentzian lineshape in the spectrum from the frequency domain. The peak area itself is proportional to the metabolite concentration.

Finally, the fourth peak parameter is the phase \( \varphi_{k} \) of the complex amplitude \( d_{k} \) via \( \varphi_{k} = \arg (d_{k} ) \). Any departure of the angle \( \varphi_{k} \) from zero indicates the degree of deviation of the kth lineshape from the bell-shaped symmetric absorptive Lorentzian profile. As per theory for (hypothetically) ideally quadrature-encoded FIDs in NMR spectroscopy, all the physical phases \( \left\{ {\varphi_{k} } \right\} \) in the input time signal should be equal to zero, \( \varphi_{k} = 0 \left( {1 \le k \le K} \right). \) In any realistic measurement, however, various experimental imperfections and uncertainties produce non-zero phases, \( \varphi_{k} \ne 0 \) for most, if not all of the resonances. Among these uncertainties are time delays between the end of the excitation RF pulse and the beginning of the acquisition of the time signal data points. A part of the delay is also due to the recovery of the electronics (the dead time of the detector) after the shock caused by the RF pulse. For these reasons, non-zero values of phase \( \varphi_{k} \) are generally conceived as phase errors, rather than some important quantities that should be taken into account for reconstructing e.g. the concentrations of molecules. The easiest way to correct for this error is to use the phased spectrum as the sum of the real and complex component spectra multiplied by \( \cos (\varphi_{k} ) \) and \( \sin (\varphi_{k} ) \), respectively, instead of the real part of the component spectrum alone [2]. Namely, this former sum as the “phased spectrum” is in the purely absorptive mode, whereas the latter spectrum is a mixture of the absorption and dispersion for a complex response function. Rigorously, this procedure is possible only if every \( \varphi_{k} \ne 0 \) has already been reconstructed, as is indeed the case in the parametric FPT. It has been shown in Ref. [2], that the mentioned linear combination of the real and imaginary part of the component spectrum is equivalent to simply redefining all the reconstructed phases \( \varphi_{k} \) as zero \( \varphi_{k} = 0 \left( {1 \le k \le K} \right). \) This amounts to using the magnitude \( \left| {d_{k} } \right| \) as the real quantity in lieu of the complex amplitude \( d_{k} \). Earlier [20, 21], this resulting spectrum was provisionally called the “ersatz spectrum” instead of the “phased spectrum”. Although both terminologies designate the same purely absorptive spectrum, the name “phased spectrum” is nevertheless more transparent as it transcribes the alteration of the seed spectrum by exactly correcting all the phase errors of each individual component resonance.

Needless to say, such an ideal phase correction is impossible in the current MRS literature which uses merely the zeroth and the 1st order-phase correction, \( \phi_{0} \) and \( \phi_{1} \), respectively. Namely, instead of our exact phasing by the reconstructed angles \( \varphi_{k} \left( {1 \le k \le K} \right) \), the approximate phasing, as performed in the MRS community, amounts to multiplying the entire complex time signal by \( { \exp }\left( {i\phi } \right) \) where \( \phi = \phi_{1} \delta + \phi_{0} . \) Here, the compound phase correction \( \phi \) is arbitrarily assumed to depend linearly on the chemical shift \( \delta \). Clearly, \( \phi \) cannot correct for all the phase errors of many individual components. In e.g. analytical chemistry based on NMR spectroscopy, several hundred resonances are not infrequently encountered in spectra reconstructed from time signals with some \( 10^{4} - 10^{5} \) measured FID data points. The same can also be said for ICRMS [13]. As a result, the zeroth and the 1st-order phase correction \( { \exp }\left( {i\phi } \right) \) leaves many phase errors uncorrected or wrongly modified across the spectrum. Such a spectrum is, therefore, never positive definite throughout the Nyquist range. By contrast, our outlined procedure always generates the strictly positive-definite, “phased spectrum”. Moreover, since such a phased spectrum is in a purely absorption mode, its interpretation is straightforward for the peak parameters and especially for the resulting peak areas and metabolite concentrations.

2.4 Derivative fast Padé transform, dFPT

Recently [1, 2], we found yet another Padé-based way to improve resolution and suppress noise, i.e. to simultaneously overcome the two major obstacles to quantitative analysis of spectra. This alternative methodology called derivative signal processing belongs to non-parametric, i.e. shape estimation from the onset. Yet, at the end of the computations, all the physical resonances are adequately parametrized, although the quantification problem itself is not actually solved. This remarkable occurrence makes the division between the parametric and non-parametric processors elusive and eventually obsolete.

In order to explain this success of the derivative Padé-based non-parametric estimation, we resorted in Ref. [2], to parametric processing. This was necessary since if a signal processor, initially conceived for shape estimation alone, is to yield the correct peak parameters, then in the end, these findings must be validated by parametric analysis which gives the explicit solution of the quantification problem within the same processor. In our first paper [1] on this topic, we presented solely the lineshapes of the envelopes in the non-parametric dFPT. In the second paper [2], we explained how the peak parameters can actually be determined from the magnitude modes of the envelope profiles in the non-parametric dFPT. Moreover, in that same study [2], we derived the explicit expressions showing how the peak parameters (positions, widths, heights) of the general mth derivative envelope for the magnitude mode in the non-parametric dFPT are connected to the absorptive peak positions, widths and heights from the non-derivative (\( m = 0 \)) parametric FPT. In addition to the magnitude modes in the non-parametric dFPT from Ref. [1], we also presented in Ref. [2] the lineshapes from the real (”absorptions”) and imaginary (”dispersions”) parts of the associated complex envelopes. We have shown in Ref. [2] that the central lobes of both these real (\( m \) even) and imaginary (\( m \) odd) derivative envelopes from the non-parametric dFPT can correctly reconstruct the exact input data for the peak heights. In other words, these “absorptions” and “dispersions” carry the same information. This is in accordance with the Kramers-Kronig general relations that permit obtaining the spectrum in an absorption from a dispersion mode and vice versa. Although the peak heights and positions were seen as being correct in the non-parametric dFPT from Refs. [1, 2], it nevertheless remained to be verified whether the peak widths are also adequate. From a fuller visualization context of the overall validation of the findings from Refs. [1, 2], it would, therefore, be advisable to determine whether the entire non-parametrically reconstructed derivative envelope lineshapes could be appropriately scrutinized and eventually confirmed as fully correct. The best way to perform such a stringent test is to compare the results from the non-parametric and the parametric dFPT, as done in Sect. 3 of the present study. To this end, we have performed the computation of the component spectra in the parametric dFPT. The goal is to see how the component spectra in the parametric dFPT agree with the envelopes from the non-parametric dFPT. Moreover, in the parametric dFPT, we also generate the total shape spectra (or envelopes) as the sum of their derivative components.

Both the non-parametric and the parametric dFPT possess the analytical expressions for the general mth order derivatives of the complex spectra. These closed formulae are obtained by applying the derivative operator \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \) to the envelopes in the non-parametric and parametric dFPT, as well as to the components of the parametric dFPT via:

Because of a more general significance going well beyond the topic of derivative signal processing, the mathematical derivation of the analytical expressions corresponding to (2.24)–(2.26) is deferred to a separate paper to be published soon. Presently, as in Refs. [1, 2], we will simply employ these analytical expressions. In the computations, for the purpose of cross-validation, we shall use all the variants of the fast Padé transform, i.e. the \( {\text{FPT}}^{\left( \pm \right)} \) and \( {\text{dFPT}}^{\left( \pm \right)} \) in both the parametric and non-parametric formulations. All the results from the \( {\text{FPT}}^{\left( + \right)} \) and \( {\text{FPT}}^{\left( - \right)} \) ought to coincide and so must those from the \( {\text{dFPT}}^{\left( + \right)} \) and \( {\text{dFPT}}^{\left( - \right)} \). After verifying that these conditions are fulfilled, it suffices to present only the “+” or “−” versions of the FPT and dFPT. This will be done with the \( {\text{FPT}}^{\left( - \right)} \) and \( {\text{dFPT}}^{\left( - \right)} \) in Sect. 3, so as to cohere with Refs. [1, 2].

As mentioned, progress in all spectroscopy is hampered by the presence of overlapping resonances. Particularly in the MRS literature, this problem is viewed as practically unsolvable [22]. To cope with the spectrally crowded envelopes with tightly-overlapped resonances, that are often indiscernibly glued to each other, research practitioners have resorted to “short-cuts”. They encode FIDs at long echo times (TE) to eliminate the faster decaying resonances. The resulting envelopes at e.g. TE = 136 or 272 ms (or the likes) are sparser than those at a short TE, e.g. 22 ms. The intention of this practice is to generate the Fourier total shape spectra with as many as possible semi-isolated resonances that could presumably be more amenable to fitting by Lorentzians, Gaussians (or their linear combination approximating Voigtians) in attempts to extract the peak parameters and, ultimately, the metabolite concentrations. While the aim of obtaining isolated resonances when trying to avoid overlapping peaks is an understandable motivation, the mentioned encoding protocols to achieve such a goal are, in fact, ill-designed. The first reason is that even at longer TEs, there will always be overlapping resonances, often located at chemical shifts that house some of the diagnostically-informative metabolites. The second reason is that even the resonances that are subjectively visualized as being isolated, single peaks at long TEs, could well be comprised of their constituent components. In such cases, the estimated metabolite concentrations would be incorrect for the resonances that are assumed to be structureless, i.e. with no sub-peaks. Thus, altogether this approach of using long TEs, while encoding FIDs, and subsequently fitting the presumed single peaks is misleading, and, as such, fruitless.

As noted, it is understandable to strive to deal with single, isolated structureless peaks, preferably in the positive-definite absorption mode to facilitate interpretation and extraction of metabolite concentrations. We can take that justifiable motivation at its “face value” and ask the question of paramount importance as to whether there could be a mathematically rigorous way to convert a congested envelope (from FIDs encoded at short, middle or long TEs) by lineshape estimation alone into its components, comprised exclusively of single isolated resonances. The answer is in the affirmative for both noise-free and noisy synthesized times signals, provided that the non-parametric dFPT is employed [1, 2]. With this prospect at our disposal, the aim of quantitatively tackling only single, isolated peaks comes to its full fruition. This is where the non-parametric dFPT comes to rescue the long-lasting trouble in the MRS community regarding the earlier stated insurmountable problem of faithfully disentangling all the overlapping peaks in the total shape spectra without recourse to the explicit quantification problem. Extension of this powerful methodology to FIDs encoded at short, middle and long TEs is currently underway, and the results shall be reported shortly.

The mentioned mathematical reason behind the overlapping resonances is a considerable spread of the long tails of resonances, especially those that are strong (tall peaks) [2]. The peak overlap is further exacerbated by the constructive and destructive interference of complex-valued spectral lineshapes of neighboring resonances for \( \varphi_{k} \ne 0 \). A way to avoid this obstacle is to reduce these tails, preferably to an extent where they would be mostly “cut-off”, so as to scale the dominant content around the resonant frequency. This amounts to a simultaneous flattening of the tails and narrowing of the peak widths. Both effects would synergistically improve SNR. Firstly, lowering the tails would reduce the baseline, which, by definition, is the averaged noise in the spectrum. Secondly, the physical signal would be enhanced, since narrowing of the resonance widths would concomitantly result in increasing the peak heights. To put this plausible concept into practice, it remains only to find a mathematical tool which accomplishes the outlined simultaneous width narrowing and tail flattening. Given the form (2.15) of the complex Lorentzian spectrum, it is clear that a transformed spectrum with the higher powers of the component denominators, i.e. using \( 1/\left( {z^{ \pm 1} - z_{k}^{ \pm } } \right)^{m} \) with \( m \) > 1 instead of \( 1/\left( {z^{ \pm 1} - z_{k}^{ \pm } } \right) \) would yield narrower, taller peaks with reduced tails. These higher powers \( 1/\left( {z^{ \pm 1} - z_{k}^{ \pm } } \right)^{m} \) of the component denominators \( 1/\left( {z^{ \pm 1} - z_{k}^{ \pm } } \right) \) in (2.15) are generated precisely by the derivatives \( \left( {{\text{d}}/{\text{d}}z^{ \pm 1} } \right)^{m} \) of the sufficiently high-derivative order \( m \) (for a derivative spectrum, a connecting formula is used between \( \left( {{\text{d}}/{\text{d}}{\nu} } \right)^{m} \) and \( \left( {{\text{d}}/{\text{d}}z^{ \pm 1} } \right)^{m} \)). Such an observation made in Ref. [2] led to the dFPT. The explicit computations in both Refs. [1, 2] fully confirmed this expectation.

When plotting the derivative spectra, the peak heights continue to increase with augmented differentiation order \( m \). Thus, for the purpose of monitoring stabilization and convergence of derivative spectra, a convenient normalization is necessary. This normalization can be made to the peak height of a reference peak (e.g. tetramethylsilane (TMS) in organic solvent or trimethylsilyl-tetradeutero-propionic acid (TSP), as used for the in vitro study of Ref. [23]) or the “pivotal peak” at one of the far ends of the given spectrum (e.g. lactate (Lac) in the case of the applications from Refs. [1, 2], and herein as well). Such a procedure, besides being useful in monitoring stabilization/convergence, is also helpful for the main practical purpose, since the relative concentration of the given metabolite can be read off directly from the normalized peak height and the corresponding FWHM. Passing from such a relative to an absolute concentration of the metabolite is made through multiplication of the former quantity by the number of the MR active nuclei and by the concentration of the chosen reference metabolite [2].

2.4.1 Stopping criteria for completing the computation of the derivative envelopes

There remains a question of certain practical significance: where do we stop with the estimation by the non-parametric dFPT alone? Earlier [1, 2], as well as presently, we used an exhaustive range for \( m \left( {1 \le m \le 50} \right). \) The question is then which of these widely varying derivative orders should be deemed sufficient to formulate a stopping criterion for ending the computations when only the non-parametric dFPT is employed? First of all, such computations in the non-parametric dFPT are so fast, that covering literally hundreds of values of \( m \) would be of no concern whatsoever regarding the used central processing unit (cpu) time on customary personal computers (the usual lap-tops of researchers). Nevertheless, the stopping criterion to finish the computations can be given as that value of \( m \) at which the normalized envelopes from the non-parametric dFPT have converged to a stable distribution throughout the spectral region of interest (SRI). Such a steady distribution of physical resonances can provide the reliable peak parameters of the magnitudes of the derivative spectral profiles, to be related to the absorptive peak positions, widths and heights. We know that the spectral structures that survive are physical (genuine), whereas the unphysical ones (spurious) are washed out by the Padé-implemented higher-order derivative operator \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} . \)

2.4.2 Peak parameters for the non-parametric dFPT

The key stage of estimation at which the non-parametric dFPT for the magnitudes of envelopes becomes quantitative is at the point where this processor can deliver the most correct numerical values of the peak positions, widths and heights as in the physical absorptive non-derivative \( \left( {m = 0} \right) \) component spectra. This is achieved by computing the normalized magnitude spectrum (stabilized with increased derivative order \( m \)) at a very dense grid of the sweep frequencies \( \nu \). Unlike Fourier processing, in the Padé-based estimation, such a grid is entirely independent of the total duration (or acquisition) time \( T \) of the input FID. This computation simultaneously yields two resonance parameters. One of these parameters is the mth derivative peak height \( h_{k,m} \) as \( h_{k,m} = {\text{max }}\left\{ {\left| {{\text{D}}_{\nu }^{m} \left( {P_{K} /Q_{K} } \right)} \right|} \right\}_{{\nu \in {\text{SRI}}_{k} }} \). Here, \( {\text{SRI}}_{k} \) is the spectral region of interest for the sweep frequencies \( \nu \) around the detected kth peak of the mth derivative profile. The other parameter is the corresponding location of the kth peak height, \( \zeta_{k,m} = {\text{Re}}\left( {\nu_{k,m} } \right) \), which is the value of the sweep frequency \( \nu \) at which the magnitude profile attains its maximum. The associated complex frequency \( \nu_{k,m} \) is \( \nu_{k,m} = \zeta_{k,m} + i\lambda_{k,m} \), where \( \lambda_{k,m} \) is the halved FWHM of the corresponding peak of the magnitude profile, \( \lambda_{k,m} = (1/2)\left\{ {\text{FWHM}} \right\}_{k,m} \). The latter profile breadth \( \left\{ {\text{FWHM}} \right\}_{k,m} \) is provided by the non-parametric dFPT as the distance between the two crossing points of the mth derivative envelope and the straight line parallel to the chemical shift abscissa at the level of the halved maximum, \( (1/2)h_{k,m} \). The earlier found peak location \( \zeta_{k,m} \) can be verified to coincide with the middle point of the mentioned two crossings. It is in this way that the same program for the dFPT, which non-parametrically computes the lineshape of the magnitude spectrum \( \left| {{\text{D}}_{\nu }^{m} \left( {P_{K} /Q_{K} } \right)} \right|, \) also automatically generates the peak positions \( \left\{ {\zeta_{k,m} } \right\} \), half-widths \( \left\{ {\lambda_{k,m} } \right\} \) and heights \( \left\{ {h_{k,m} } \right\} \left( {1 \le k \le K} \right) \) of the kth resonance in the mth derivative envelope. This is the case because for sufficiently high derivative order \( m \), the profiles of the envelope \( \left| {{\text{D}}_{\nu }^{m} \left( {P_{K} /Q_{K} } \right)} \right| \) and its kth component \( \left| {{\text{D}}_{\nu }^{m} \left( {P_{K} /Q_{K} } \right)_{k} } \right| \) coincide locally around the kth peak, as will be illustrated in Sect. 3. Once the derivative set \( \left\{ {\zeta_{k,m} , \lambda_{k,m} , h_{k,m} } \right\} \) is produced by the explained procedure, the corresponding non-derivative \( \left( {m = 0} \right) \) absorptive peak positions \( \zeta_{k } \), widths \( \lambda_{k} \) and heights \( h_{k} \) are immediately deduced from the explicit relations from Ref. [2] summarized as:

2.4.3 Advantage of the dFPT over the dFFT regarding noise

The most favorable aspects of derivative signal processing are best appreciated when viewed in the dFPT. This becomes most evident when the input time signals are noise-corrupted, which mimics the situations with encoded FIDs. It has been demonstrated in Ref. [1] that the dFPT suppresses noise, which is, however, amplified by the dFFT. In other words, the same operator \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \) has two opposite effects by acting as a noise filter in the dFPT and as a noise amplifier in the dFFT. Namely, since \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \exp \left( { - 2\pi i\nu n\tau } \right) = \left( { - 2\pi in\tau } \right)^{m} \exp \left( { - 2\pi i\nu n\tau } \right), \) it follows that the dFFT, generated using (2.19), does not process the time signal \( c_{n} \) directly, but rather its windowed (weighted) counterpart \( ( - 2\pi in\tau )^{m} c_{n} \) according to:

Here, the window function \( w = \left( { - 2\pi in\tau } \right)^{m} \) puts a higher weight on the tail parts of the time signals, that in encoded FIDs are comprised mainly of noise. In other words, the dFFT emphasizes the noisy part of the encoded FID and de-emphasizes the earlier encoded data points that house the most intense and, thus, most important oscillations. In sharp contrast, the dFPT does not multiply the input time signal \( c_{n} \) by any weight function. Rather, we first process the original, raw time signal \( c_{n} \) by the non-derivative \( (m = 0)\,{\text{FPT}} \) and it is to the resulting analytical expression \( P_{K}/Q_{K} \) that we apply the \( {\text{D}}_{\nu }^{m} \) operator (2.23) to arrive at the dFPT as in (2.24). This latter step further reduces noise as per the features of derivative envelopes by reference to the discussed width narrowing and the tail flattening.

3 Results

Here, we shall report on the results from the present computations using both the non-parametric and the parametric versions of the dFPT. Of course, we could have chosen any relevant application, but for continuity of the systematic scrutiny, we opt to complement our two previous studies [1, 2] and, thus, we shall also focus here on MRS for breast cancer diagnostics. To this end, we will employ the same synthesized time signals as in Refs. [1, 2], all in accordance with the corresponding encoded FIDs reported in Ref. [23]. Both noiseless and noisy simulated time signals will be subjected to non-derivative and derivative estimations.

Here, the real part of \( \nu_{k} , \,{\text{i}}.{\text{e}}. \)\( {\text{Re}}\left( {\nu_{k} } \right), \) is the chemical shift, whereas the imaginary part, \( {\text{Im}}\left( {\nu_{k} } \right), \) is proportional to the reciprocal of the spin–spin relaxation time \( T_{2k}^{*} {\text{via Im}}\left( {\nu_{k} } \right) = 1/(2\pi T_{2k}^{*} ) \) for the kth resonance. Both \( {\text{Re}}\left( {\nu_{k} } \right) \) and \( {\text{Im}}\left( {\nu_{k} } \right) \) are expressed in dimensionless units called parts per million (ppm). The generally complex amplitude \( d_{k} \), as the magnitude \( \left| {d_{k} } \right| \) multiplied by the phase factor, \( d_{k} = \left| {d_{k} } \right|\exp \left( {i\varphi_{k} } \right), \) is presently set to be real \( \left( {d_{k} = \left| {d_{k} } \right|} \right) \) by choosing \( \varphi_{k} = 0\left( {1 \le k \le 9} \right). \) The magnitudes of amplitudes are in arbitrary units (au). Also, in the sequence (3.1) are the peak heights of metabolites alongside their abbreviations, whose full names are given in the list of abbreviations. The kth metabolite concentration is proportional to the peak area ak, which is equal to \( \left| {d_{k} } \right|/2 \) for an absorptive Lorentzian non-derivative \( \left( {m = 0} \right) \) lineshape as per (2.22). Since all the peak widths are chosen to be the same in (3.1), it is seen therein that lactate, Lac, is the most abundant metabolite which resonates at 1.332 ppm. The listed chemical shifts cover the band of the fundamental frequencies from 1.332 to 3.281 ppm. In the spectral range [1.332, 3.281] ppm, the narrow interval [3.219, 3.222] ppm is of our primary interest because it contains the overlapping PE and PC resonances, separated by a mere 0.001 ppm. The peak heights of PE and PC, being vastly different, 176.715 versus 23.562 au, respectively, represent the biggest challenge to quantitative visualization by estimations when only envelopes are used. The specific values of the parameters in (3.1) are derived from the corresponding estimation with encoded time signals from Ref. [23].

These latter encoded FIDs [23] were of long total length with 65536 data points that the Fourier frequency grid necessitate and yet the ensuing resolution is low in the FFT. In Ref. [23], the bandwidth, BW, was 6000 Hz, resulting in the sampling time \( \tau = \left( {1/6000} \right) \) s. The static magnetic field strength was \( B_{0} \approx 14.1 \) T, which corresponds to the Larmor frequency \( \nu_{\text{L}} = 600 \) MHz. Since we have in mind the FPT and dFPT, the signal length can be much shorter and, thus, we initially set \( N = 2048. \) The entire computation was repeated with \( N = 4096 \) to test the convergence and stability of all the reconstructions in the non-parametric and parametric FPT and dFPT.

Regarding the encoded FIDs and the related spectra [23], the chemical shifts are measured relative to the location (taken to be at the zero frequency) of TSP (not present in the tissue). Relative to this zero frequency of TSP, we set the water resonance to be positioned at \( \nu_{{{\text{H}}_{2} {\text{O}}}} \) (ppm) = 4.68 ppm. All the other chemical shifts \( \nu \) (ppm) are obtained from their counterparts \( \nu \) (Hz) by using the relation:

where, of course, no translational scaling by \( \nu_{{{\text{H}}_{2} {\text{O}}}} \) (ppm) appears in the resonance widths from (3.4).

As to the synthesized noisy time signals, these are generated following Ref. [1]. Therein, a noise model is used according to the following prescription. Given that the main uncertainties in encoded FIDs come from the unavoidable presence of random noise, the synthesized noiseless time signals should be corrupted with a similar type of stochastic perturbations. This is modeled in Ref. [1] by adding complex random noise \( \left\{ {r_{n} } \right\} \) to the noiseless complex FID data points \( \left\{ {c_{n} } \right\} \). Both data sets \( \left\{ {r_{n} } \right\} \) and \( \left\{ {c_{n} } \right\}\,\left( {0 \le n \le N - 1} \right) \) must be of the same length N. The added disturbance \( \left\{ {r_{n} } \right\} \) is random Gaussian zero-mean white noise with a prescribed standard deviation \( \sigma \) (the square root of the variance). In Ref. [1], the simulated noisy time signal \( \left\{ {c_{n} } \right\} + \left\{ {r_{n} } \right\} \) has been processed with \( \sigma = 0.0289 \) RMS, where the root mean square (RMS) error refers to the noiseless part \( \left\{ {c_{n} } \right\}, \) i.e. \( {\text{RMS}} = \mathop \sum \nolimits_{n = 0}^{N - 1} \left| {c_{n} } \right|^{2} /N \). This same standard deviation \( \sigma \) will also be used in the present work.

3.2 Output data

In Sect. 3.2.1, we shall first recapitulate what is already partially known from Refs. [1, 2] on the non-parametric dFPT. To this, we shall add the previously non-displayed lineshapes (on the same graphs) of the magnitudes, absorptions and dispersions, so as to make in evidence their combined patterns and the interference effects. Section 3.2.2 is devoted to a detailed comparison of the envelope lineshapes reconstructed by the non-parametric dFPT with the corresponding components from the parametric dFPT. In the present work, all the computations in the FPT and dFPT have been carried out for both noiseless and noisy time signals. The obtained results are indistinguishable from each other (as was also seen in Ref. [1]), such that the presently-reported spectra refer simultaneously to the Padé-based reconstructions for the noiseless and noisy time signals.

3.2.1 Processing by the dFPT with comparisons to the dFFT

3.2.1.1 The problem of resolution and noise in spectra

The frequency resolution \( \Delta \nu_{\text{FFT}} \) in the spectral envelopes from the FFT is limited by the total signal length \( N \), or equivalently, by the total acquisition time \( T. \) This is implied by the relations \( \Delta \nu_{\text{FFT}} \equiv \nu_{k + 1}^{\text{FFT}} - \nu_{k}^{\text{FFT}} = N/\tau = 1/T \), where \( \nu_{k}^{\text{FFT}} \) is the Fourier grid, \( \nu_{k}^{\text{FFT}} = k/T \left( {0 \le k \le N - 1} \right). \) To approximately determine the spectral poles and the corresponding residues (amplitudes), most researchers in the MRS community first plot the FFT spectrum and then fit the peaks seen therein by some model functions (Lorentzians, Gaussians, …). By contrast, the parametric FPT exactly determines the poles and residues of the spectrum prior to ever plotting the spectrum. This is done by exactly solving the secular equation (the characteristic polynomial equation) in (2.11) which is equivalent to diagonalizing the associated sparse Hankel matrix treated in our Padé-based reconstructions. In fact, for a given set \( \left\{ {a_{j} } \right\} \) of a polynomial \( A_{M} \left( u \right), \) e.g. the Matlab command “roots \( \left( a \right) \)” returns the roots of the equation \( A_{M} \left( u \right) = 0 \) as the eigenvalues of the eigenproblem of the corresponding Hessenberg matrix. The roots of the secular equations are the spectral poles. According to (2.12), the corresponding amplitudes are obtained analytically as the Cauchy residues of the response function given by the rational polynomial in the FPT spectrum. Once these complex fundamental frequencies \( \left\{ {\nu_{k} } \right\} \) (poles) and the complex amplitudes \( \left\{ {d_{k} } \right\} \) (whose absolute values are proportional to the concentrations of the detected metabolites) are reconstructed, the spectrum can be plotted in any desired mode (absorption, dispersion, magnitude, power). Such a Padé-based strategy improves the resolution beyond \( \Delta \nu_{\text{FFT}} \) and, moreover, this is achieved with much shorter signal length than that required by the FFT.

The mentioned characteristics (increased resolution, removed noise) are illustrated for both non-derivative and derivative lineshapes in the FPT \( (m = 0) \) and dFPT \( (m > 0) \) in Figs. 1 and 2, respectively. Thus, Fig. 1 starts with the noisy FID on panel (a). Therein, only the real part of the complex time signal is shown (the imaginary part is of a similar waveform, except for the standard phase mismatch by \( \pi /2 \) radians). The displayed FID used for the Padé-based processing is of the length \( N = 2048. \) This is some formidable 32 times shorter than the corresponding 65536 FID data points required by Fourier analysis from Ref. [23]. Shorter FIDs are advantageous, as they imply faster encoding and this, in turn, leads to more expedient protocols for the patient undergoing the MR scanning. The un-normalized spectral envelopes in the absorption mode from the non-parametric FPT are depicted on panels (b) and (c) in Fig. 1 for all the nine resonances. The peak of lactate, Lac, from panel (b) appears as the dominant resonance from 1.30 to 3.29 ppm. Note that the full names for the metabolites, abbreviated in the figures, are in the list of abbreviations. The input peak heights \( \left\{ {h_{k} } \right\}\,(1 \le k \le K, K = 9) \) from (3.1) in panels (b) and (c) fully match the maximae of the absorptive lineshapes of envelopes for all the seven isolated, single resonances. The exception being the two overlapped resonances (PC, PE) that are seen as a single, amalgamated peak in which PC is completely invisible. Its presence in Fig. 1c is indirectly suggested by the discrepancy between the input peak height for PE and the profile maximum at the chemical shift 3.221 ppm of phosphoethanolamine. The un-normalized absorptive spectra from the parametric FPT are shown in panels (d) and (e) of Fig. 1. Here, the component spectra on panel (d), of course, delineate all the individual resonances, including the two separate peaks at 3.220 and 3.221 ppm for PC and PE, respectively. On panel (e), the total shape spectrum from panel (c) is superimposed on top of the component lineshapes from panel (d). The components and the envelope coincide for the lineshapes above the baseline for the five isolated resonances. The only difference is at the PC–PE location. Therein, the envelope is almost entirely determined by the PE component, which dominates the hidden PC peak. The envelope from panels (c) or (d) reconstructed by the non-parametric FPT is verified to be identical to the envelope (the sum of the component spectra) from the parametric FPT. The same holds true for panel (b).

Non-derivative FPT with a noisy time signal or FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \) for exact reconstructions of the un-normalized non-parametric envelopes and parametric components. The real part of the FID is on (a), spectra are on (b–e). The envelopes are on (b), (c) and (e), while the components are on (d) and (e). The envelope on (c) cannot resolve the recognized cancer biomarker phosphocholine, PC (3.220 ppm) masked by phosphoethanolamine, PE (3.221 ppm). These delineated overlapped peaks are seen in the components on (d) and (e) (Color figure online)

Non-parametric estimations of envelopes for a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \) by 2 processors: comparison of the dFPT (b–d) with the dFFT (e). The un-normalized non-derivative \( \left( {m = 0} \right) \) magnitude profile in the FPT is on (a), while (b–e) are for magnitude modes, normalized to the tallest peak (lactate, Lac) from Fig. 1b. A striking resolution improvement and noise suppression in the dFPT is seen on (b) and (d) with the PC–PE double peak resolved for high derivatives of envelopes (e.g. \( m = 50 \)). The dFFT amplifies noise even at lower derivative orders (e.g. \( m = 6 \)) for the normalized magnitude envelope on (e) (Color figure online)

Figure 2 is for the normalized total shape spectra, reconstructed by the dFFT \( \left( {m = 6} \right) \), as well as by the non-parametric FPT (\( m = 0 \)) and dFPT \( \left( {m = 6, m = 50} \right). \) The magnitude mode is shown alone. On panel (a), the non-derivative (\( m = 0 \)) envelope in the FPT appears as a broad single peak due mainly to PE. This lineshape does not have any noticeable curvature change around the PC location. This perspective is dramatically altered for the better in the derivative envelopes from the dFPT \( \left( {m = 6, m = 50} \right) \) on panel (b). Therein, at \( m = 6 \), a visible shoulder emerges at the PC position alongside the clear PE peak. This shoulder eventually becomes a completely delineated PC resonance, whose maximum entirely coincides with the associated input peak height, just as in the case of PE on the same panel (b). A more expanded presentation of the comparison between the lineshapes for the two derivative orders \( m = 6 \) and \( m = 50 \) in the dFPT is given on panels (c) and (d), respectively. Relative to Fig. 1c for the absorptive envelope in the FPT (\( m = 0 \)), the peak breadths are noticeably narrower for the magnitude modes of the six (out of seven) well-delineated resonances (Cho, PE, GPC, β-Glc, Tau, m-Ins) in the dFPT \( \left( {m = 6} \right) \) on Fig. 2c. As to the PC shoulder, it is seen to be leaning towards the lower part of the PE lineshape. All the input peak heights for the mentioned six profiles exactly match the maximae of the envelopes at the resonance chemical shifts. The top of the PC shoulder is seen to overestimate the input peak height at 3.220 ppm for \( m = 6 \) on panel (c). However, on panel (d) with \( m = 50 \) for the dFPT, the PC and PE peaks are fully separated all the way down to the baseline level. In Fig. 2d, all seven resonances in the displayed window [3.205, 3.290] ppm are extremely thin, especially when compared to the non-derivative (\( m = 0 \)) absorption from Fig. 1c. The results in the FPT and dFPT on panels (a)–(d) are checked to be identical for the noise-free \( \left( {\sigma = 0} \right) \) and noise-corrupted \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \) input time signals. This implies that both the FPT and the dFPT effectively eliminate the entire noise. Particularly, the dFPT is remarkably successful not only in complete noise removal, but also in a total flattening of the long tail ends of all the resonances. This is best appreciated by comparing the Lorentzian tails near the resonance frequencies in Fig. 1c with Fig. 2d. In Fig. 1c, with the FPT (\( m = 0 \)), the lower parts of all the single peaks are quite wide. Moreover, in the same Fig. 1c, they mix together for e.g. m-Ins and Tau through an interference effect which lifts the baseline above the level of the abscissa for the chemical shifts. All such interference effects have completely disappeared in Fig. 2d for the dFPT \( \left( {m = 50} \right) \), where the Lorentzian tails of every single (Cho, GPC, β-Glc, Tau, m-Ins) and initially overlapped (PC, PE) resonances are entirely immersed into the baseline, which is itself embedded into the chemical shift axis. This means that the dFPT for high derivative order \( m \) hugely improves the SNR compared to the non-derivative FPT \( \left( {m = 0} \right) \). In other words, the differentiation transform \( {\text{D}}_{\nu }^{m} \) from (2.23) and (2.24) acts as a powerful noise filter in the dFPT.

This is in sharp contrast to the dFFT from (2.27), where the same derivative transform \( {\text{D}}_{\nu }^{m} \) amplifies noise, as illustrated for \( m = 6 \) on panel (e) of Fig. 2. Such a breakdown of the dFFT for noisy FIDs is checked to worsen for \( m > 6 \) (not shown to avoid clutter with the uninformative Fourier processing). This failure of the dFFT is rooted in the fact that the operator \( {\text{D}}_{\nu }^{m} \) yields the time-dependent factor \( \left( {n\tau } \right)^{m} \) for \( m > 0 \) which weights heavily the noisy portion with larger values of \( n\tau \) in the synthesized \( \left\{ {c_{n} + r_{n} } \right\} \) or encoded FIDs. With this feature, as stated earlier, the derivative operator \( {\text{D}}_{\nu }^{m} (m > 0) \) in the dFFT emphasizes the unphysical, noisy tail of the FID and de-emphasizes the physical, more intense time signal oscillations towards the beginning of the data acquisition.

Such a drawback which flagrantly invalidates the dFFT is absent from the dFPT, which first uses the non-derivative FPT \( \left( {m = 0} \right) \) to process the original time signal, e.g. \( \left\{ {c_{n} + r_{n} } \right\} \), unaltered by the operator \( {\text{D}}_{\nu }^{m} . \) Subsequently, the derivative spectrum in the dFPT \( (m > 0) \) is generated by subjecting the mentioned seed spectrum from the FPT \( \left( {m = 0} \right) \) to the operator \( {\text{D}}_{\nu }^{m} . \) The second step does not introduce any new noise. This occurs because the results of the application of \( {\text{D}}_{\nu }^{m} \) to the closed formulae \( P_{K}^{ \pm } /Q_{K}^{ \pm } \) for the seed spectra in the FPT are also the analytical expressions for any derivative order \( m. \) In other words, since the analytical expressions for \( {\text{D}}_{\nu }^{m} P_{K}^{ \pm } \left( {z^{ \pm 1} } \right)/Q_{K}^{ \pm } \left( {z^{ \pm 1} } \right) \) are the working formulae in the Padé derivatives, the dFPT never actually employs any numerical differentiation, which is known to be prone to computational errors. Moreover, the derivative operator \( {\text{D}}_{\nu }^{m} \) itself removes noise when applied to \( P_{K}^{ \pm } /Q_{K}^{ \pm } . \) The reason for this is that the analytical expressions for \( {\text{D}}_{\nu }^{m} P_{K}^{ \pm } /Q_{K}^{ \pm } \) are also the ratios of two polynomials. It is this non-linearity feature of the polynomial quotients in the dFPT which gives the noise filtering property to the derivative operator \( {\text{D}}_{\nu }^{m} \). Thus, any trace of noise, which is eventually left after passing the input data \( \left\{ {c_{n} + r_{n} } \right\} \) through the zeroth order (\( m = 0 \)) filters \( P_{K}^{ \pm } /Q_{K}^{ \pm } \), will be further suppressed by the subsequent mth order \( (m > 0) \) filters \( {\text{D}}_{\nu }^{m} P_{K}^{ \pm } /Q_{K}^{ \pm } . \) After two such consecutive rational polynomial filters, noise has no chance of surviving to any noticeable extent, especially for large values of \( m . \) This is what is evidenced on panel (d) of Fig. 2. Similarly, the unprecedented noise removal capability of the dFPT is also encountered for other more noise-corrupted FIDs with higher standard deviation \( \sigma \) than \( \sigma = 0.0289 {\text{RMS}} \) used in the present illustrations.

Our aim in this subsection is to superimpose the different modes of spectral lineshapes on the same subplots of multi-trace figures. For systematics, we first do so by simultaneously displaying on Fig. 3 the magnitudes and the real parts of the complex envelopes form the non-parametric dFPT. This is followed by Fig. 4, where each subplot contains derivative lineshapes for magnitudes and the imaginary parts of the complex envelopes. Finally, Fig. 5 shows all three derivative envelope lineshapes (magnitudes, real and imaginary parts). Figures 3, 4 and 5 are for normalized envelopes associated with the varying model orders \( m = 8\left( 8 \right)48. \) Note that for \( \varphi_{k} = 0, \) the real part of the complex non-derivative Lorentzian spectrum (\( m = 0 \)) is a bell-shaped single absorptive symmetric profile. The imaginary part of the same complex spectrum (\( m = 0 \)) is the single dispersive asymmetric lineshape for \( \varphi_{k} = 0 \). These notions cease to apply to derivative spectra of the general order \( m > 0. \) For example, taking the real part of a complex derivative envelope for \( m \) even and odd results in the parity change, since maximae become minimae and vice versa. Moreover, peaks and dips appear for a single kth resonance in either the real or imaginary derivative lineshapes for any order \( m. \) Nevertheless, for an easier and shorter nomenclature, despite the differences among non-derivative (\( m = 0 \)) and derivative (\( m > 0 \)) envelopes, we can still preserve the adjectives “absorptive” and “dispersive” for the real and imaginary derivative lineshapes, as long as we are aware of the meaning behind the quotation marks for these more familiar notions.

Non-parametric dFPT for a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \): comparison of 2 modes, the real parts of normalized envelopes versus their magnitudes. Variation of the derivative order \( m \) from 8 to 48 in the step of 8, i.e. \( m = 8\left( 8 \right)48 \). Magnitude modes exhibit a single peak per resonance. The real parts show multiple lobes, the strongest of which are at the resonant frequencies and with the peak heights coinciding with the associated input data. No positive-definite lobe of the envelope real parts is taller than the magnitude lineshapes (i.e. the lineshapes of the real parts of complex envelopes are all confined within the magnitude profiles) (Color figure online)

Non-parametric dFPT for a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \): comparison of 2 modes, the imaginary parts of normalized envelopes versus their magnitudes. Variation of the derivative order \( m \) from 8 to 48 in the step of 8, i.e. \( m = 8\left( 8 \right)48 \). The imaginary parts show multiple lobes, the upward-oriented intensities of which are always weaker than the maximum of the single peak per resonance in the magnitude mode (Color figure online)

Non-parametric dFPT for a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \): comparison of 3 modes, the real and imaginary parts of normalized envelopes versus their magnitudes. Variation of the derivative order \( m \) from 8 to 48 in the step of 8, i.e. \( m = 8\left( 8 \right)48 \). For the same derivative order m, the interference of the real and imaginary parts of the given envelope reproduces exactly the lineshape profile of the associated magnitude mode (Color figure online)

Naturally, no magnitude mode profile, being strictly positive-definite, can ever cross the abscissae of chemical shifts. This results in the emergence of the single peak per resonance in the spectral magnitude mode of lineshapes, as is clear from Figs. 3, 4 and 5. By contrast, the “absorptions” (Fig. 3) and “dispersions” (Fig. 4) have multiple lobes. These correspond to zeros (per the single kth resonance) of the numerator polynomial of the real or imaginary parts of the complex derivative envelopes for \( m > 0 \) in the dFPT. In the derivative “absorptions” (Fig. 3), the central lobe located at the kth fundamental chemical shift, \( {\text{Re}}\left( {\nu_{k} } \right), \) has the peak heights which coincide with those of the magnitude mode and the corresponding input data for \( h_{k} \) from (3.1). The “absorptive” side lobes (Fig. 3) are placed symmetrically and equidistantly around the central lobe. The heights of the side lobes are systematically and gradually lower than the central lobe with the increased separation from the central chemical shift \( {\text{Re}}\left( {\nu_{k} } \right) \) of the kth resonance.

Both maximae (peaks) and minimae (dips) emerge per each kth “absorptive” lineshape in Fig. 3. With increased derivative order \( m \), Fig. 3 shows that all the “absorption” lobes are narrowed. This follows the same pattern exhibited by the magnitude modes (except that the latter are the strict, single maximum per kth resonance). The maximae of every such “absorption” lobe is strictly bounded by the lineshape of the magnitude mode. In other words, the kth magnitude of the mth-order derivative envelope is itself the envelope for all the maximae (peaks) of the “absorptive” lineshape lobes (Fig. 3).

The situation with the “dispersions” (Fig. 4) differs from their “absorption” counterparts in several aspects. The “dispersive” maximae lobes are also all confined within the magnitude modes. However, their heights are considerably smaller than any point on the envelope for the magnitude mode. Moreover, the highest maximum (peak) of the “dispersion” lobes is not placed at the location \( {\text{Re}}\left( {\nu_{k} } \right) \) of the kth resonance. These features occur because we chose to present the “dispersions” in Fig. 4 for the same derivative order \( m \) (even, \( m = 2m^{\prime}, m^{\prime} = 4,8, \ldots \)) as for the “absorptions” from Fig. 3. Had we chosen the odd value of \( m\,(m = 2m^{\prime} - 1) \) for the “dispersions” in Fig. 4, their central lobes would have the peak heights matching those of the magnitude (and of the input data \( h_{k} \)), as shown in Ref. [2]. The reason for which we have presently chosen the same derivative order \( m \) for both “absorptions” and “dispersions” is to analyze their interference pattern (Fig. 5). Such an analysis would only make sense when, of course, both “absorptive” and “dispersive” lineshapes refer to the same values of the derivative order \( m. \)

The intricacies of “absorptive” and “dispersive” derivative lineshapes are best appreciated when viewed in concert on the same graph. Namely, the interference pattern of derivative “absorptions” and “dispersions” is most clearly delineated when they are both juxtaposed to the magnitude modes on the same plots, as done in Fig. 5. Therein, the phase mismatching (by \( \pi /2 \)) of “absorptions” and “dispersions” is fully reminiscent of the similar behaviors of the ordinary sine and cosine trigonometric functions. In particular, it becomes evident why the sole “absorption” recovers the entire peak height of the magnitude mode (and, in turn, of the input data \( h_{k} \)). This is due to the stated phase mismatch of \( \pi /2 \) between “absorption” and “dispersion” lobes: the former peak is at the resonant chemical shift \( {\text{Re}}\left( {\nu_{k} } \right), \) where the latter attains its zero value. More generally, the constructive-destructive interference of “absorption” and “dispersion” is at play on all the panels of Fig. 5. As a net result of this interplay, the multiple lobes of “absorptions” and “dispersions” are combined in such a way that, in the end, only a single maximum survives, and this is the peak of the magnitude mode.

3.2.2 Non-parametric dFPT versus parametric dFPT

This subsection is devoted to the critical test of the non-parametric dFPT by comparison with the parametric dFPT. Initially, the parametric non-derivative (\( m = 0 \)) FPT has already undergone the quantification step by explicitly solving the spectral analysis problem. Subsequently, the obtained fundamental frequencies and amplitudes \( \left\{ {\nu_{k} , d_{k} } \right\} \left( {1 \le k \le K} \right) \) are used to generate both the components and total shape spectra in the non-derivative parametric FPT. Finally, the analytical expressions from (2.25) for the spectra of components are subjected to the derivative operator \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \) for \( 1 \le m \le 50, \) and their sums from (2.24) are the envelopes from the parametric dFPT. This yields the derivative components and the derivative total shape spectra in the parametric dFPT. It is the plots of the latter two types of lineshapes that are to be compared with their counterparts from the non-parametric dFPT, which explicitly computes only the envelopes (total shape spectra) and not the components, as it never addresses the quantification problem. The first check along these lines is to restrict the comparisons to the envelopes alone. By doing this test, we verified that exactly the same results are obtained in the parametric and non-parametric dFPT for the derivative total shape spectra with \( 1 \le m \le 50. \) This conclusion has been checked to also hold true for \( m > 50. \) Such a finding fully coheres with the well-known fact that the parametric and non-parametric non-derivative (\( m = 0 \)) FPT yield identical envelope spectra [24, 25]. The second check is to extend the comparison to the critical case which juxtaposes the component spectra from the parametric dFPT alongside the envelopes from the non-parametric dFPT.

This is carried out in Figs. 6, 7, 8 and 9 that compare the lineshapes for the components (parametric dFPT) with the envelopes (non-parametric dFPT) for two increments or rates \( \Delta m \) of \( m \), that are set to \( \Delta m = 1 \) and \( \Delta m = 4. \) The smallest grid \( \Delta m = 1 \) is used to more closely monitor the changes of the derivative lineshapes. Some features of the profile alterations might be overlooked for larger values of the increment \( \Delta m \) of the derivative order. Thus, Figs. 6, 7 and 8 are for the rate \( \Delta m = 1 \) of the derivative order: \( m = 1\left( 1 \right)21 \) or \( m = 1,2,3, \ldots ,21. \) Figure 9 is for \( \Delta m = 4 \) between \( m = m_{ \hbox{min} } \left( {\Delta m} \right)m_{ \hbox{max} } = 24\left( 4 \right)48 \) or \( m = 24,28, \ldots ,48 \). On each panel (a) of Figs. 6, 7, 8 and 9, we plot the pure absorption un-normalized lineshapes for the non-derivative \( \left( {m = 0} \right) \) envelope (coincident for the non-parametric and parametric FPT) and its two components, respectively. This serves as a reminder of the starting lineshapes, so as to comparatively follow the changes of the spectral profiles with increased derivative order \( m \) in the dFPT.

Comparison of the parametric and non-parametric un-normalized spectra for components and envelopes, respectively, in the dFPT for lower derivative orders \( m = 1\left( 1 \right)7 \) using for a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \). The reference non-derivative (\( m = 0 \)) absorption in the FPT is on (a), while (b–h) are for magnitudes. Envelopes for m = 0(1)6 on (a–g) do not separate phosphocholine, PC, from phosphoethanolamine, PE, that are, however, clearly visualized for m = 7 on (h) with an approximate match of their lineshapes (Color figure online)

Comparison of the parametric and non-parametric un-normalized spectra for components and envelopes, respectively, in the dFPT for higher derivative orders \( m = 8\left( 1 \right)14 \) using a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \). The reference non-derivative (\( m = 0 \)) absorption in the FPT is on (a), while (b–h) are for magnitudes. For systematically increased \( m \), the envelopes from the non-parametric dFPT are seen to gradually tend to the components from the parametric dFPT, with the best outcome for \( m = 14 \) on (h) (Color figure online)

Comparison of the parametric and non-parametric un-normalized spectra for components and envelopes, respectively, in the dFPT for higher derivative orders \( m = 15\left( 1 \right)21 \) using a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \). The reference non-derivative \( \left( {m = 0} \right) \) absorption in the FPT is on (a), while (b–h) are for magnitudes. For systematically increased \( m \), the envelopes from the non-parametric dFPT are seen to gradually tend to the components from the parametric dFPT, with the best agreement (in fact, coincidence) on (g) and (h) for \( m = 20 \) and \( m = 21, \) respectively (Color figure online)

Comparison of the parametric and non-parametric un-normalized spectra for components and envelopes, respectively, in the dFPT for higher derivative orders \( m = 24\left( 4 \right)48 \) using a noisy FID \( \left( {\sigma = 0.0289 {\text{RMS}}} \right) \). The reference non-derivative (\( m = 0 \)) absorption in the FPT is on (a), while (b–h) are for magnitudes. For systematically increased \( m \), the envelopes from the non-parametric dFPT are seen to be in perfect agreement with the components from the parametric dFPT. The components are fully confluent with the envelopes and all that is seen from the former are their zero-valued tails. This implies that the non-parametric dFPT with its lineshape estimation alone is capable of exact quantification, on a component-by-component basis (Color figure online)

It is immediately clear from Fig. 6 that for increasing \( m \) the component peak heights are systematically enhanced with a concomitant narrowing of the peak widths. To emphasize the increase of the peak heights with augmentation of the derivative order \( m, \) all the component and envelope lineshapes are plotted in their un-normalized form on Figs. 6, 7, 8 and 9.

It is seen on panels (a) and (b) in Fig. 6 for \( m = 0 \) and \( m = 1, \) respectively, that the PC peak is completely masked in the envelopes by the long tail of the dominant PE resonance. This explains the fact that the sum of the PE and PC resonances appears as a single, unresolved peak. However, with the gradual increase of the derivative order \( m \), the widths of the individual resonances in the envelopes are systematically narrowed, the peak heights augmented and the wings of the tail lowered towards the baseline level. This three-fold trend evidently gears toward the peak separation with a significant diminution in the waveform overlaps between the PC and PE lineshapes. In Fig. 6, such a trend of the envelopes becomes most pronounced on panel (h) for the derivative order \( m = 7. \)

The lineshapes for the three groups with \( \Delta m = 1 \), each referring to seven derivatives, are shown in Fig. 6 (\( m = 1,2, \ldots ,7 \)), Fig. 7 (\( m = 8, 9, \ldots ,14 \)) and Fig. 8 (\( m = 15, 16, \ldots ,21 \)). In Figs. 6, 7, 8 and 9, the component spectra from the parametric dFPT are drawn as red curves, whereas the blue curves are reserved for the envelopes from the non-parametric dFPT. Returning to panel (b) on Fig. 6, we see that the 1st derivative envelope in the non-parametric dFPT is overall quite similar to the PE component peak from the parametric dFPT. In other words, this former processor in conjunction with the 1st derivative alone gives no hint of the presence of a smaller peak in Fig. 6b, other than a slight change in the curvature of the envelope at the chemical shift 3.220 ppm of PC. Such an observation is modified already with the 2nd derivative lineshapes shown in panel (c) of Fig. 6. Therein, the non-parametric envelope exhibits a protrusion which is pushed a bit away from the location 3.220 ppm of PC.

With increase of the derivative order from \( m = 3 \) to \( m = 5 \) on panels (d) to (f), respectively, an oscillatory pattern is observed in Fig. 6. Namely, for these low-order derivative envelopes in the non-parametric dFPT, the structure near the position of PC is seen as fluctuating, and it even completely disappears for \( m = 4 \) on panel (e). Had we not taken the smallest step \( \Delta m = 1, \) this fluctuation would not have been noticed. When the PC structure begins again to emerge for \( m = 5 \) on panel (f) of Fig. 6, still the ensuing envelope has merely a slight shoulder (near the PC location), which is much less delineated than the bump for \( m = 2 \) on panel (c) in the same figure. Such an instability of the secondary structure (i.e. that related to PC) in this narrow frequency range [3.219, 3.222] ppm for \( m \le 5 \) indicates that low-order derivative envelopes are unreliable if the motivation for their use is to separate the overlapping resonances. Eventually, the secondary structure near PC enters a more stable regimen for \( m = 6 \) and \( m = 7 \) on panels (g) and (h), respectively, in Fig. 6. This emergence of stabilization of the PC peak is steadily continued throughout panels (b) – (h) in Figs. 7 and 8. Taken together, Figs. 6, 7 and 8 with the minimal increment \( \Delta m = 1 \) of \( m \) indicate the need for higher-orders of derivatives of envelopes in the non-parametric dFPT. Further, the peak stabilization of the PC resonance is maintained for a larger increment \( \Delta m = 4 \), as illustrated for \( m = 24\left( 4 \right)48 \) or \( m = 24,28, \ldots ,48 \) in Fig. 9.

Regarding Figs. 6, 7, 8 and 9, the word “stabilization” refers to a stable maintenance of the well-delineated, isolated peaks. Since these lineshapes are not normalized, they are seen in Figs. 6, 7, 8 and 9 to undergo a huge increase in the peak heights, as per theory [2], when the derivative order \( m \) is augmented from \( m = 1 \) to \( m = 48 \). The outlined remarks on the stability of the PC peak with respect to the value of the derivative order \( m \) for the envelope in the non-parametric dFPT need to be quantitatively confirmed by the spectra that serve as the gold standard. This is critically important for Figs. 6, 7, 8 and 9, where the input peak heights \( h_{k} \) from (3.1) are not plotted alongside the un-normalized derivative lineshapes with gigantic maximae. The needed gold standard is the collection of the un-normalized derivative component spectra reconstructed by the parametric dFPT. Comparisons of these latter components with the envelopes provide the quantitative value of the non-parametric dFPT. To this end, it suffices to see how closely the blue (non-parametric dFPT) and red (parametric dFPT) curves for the envelopes and components, respectively, match each other for varying differentiation order \( m \). It was seen in Figs. 6 and 7 that the discrepancy between the spectra from the non-parametric dFPT (envelopes) and parametric dFPT (components) is most noticeable at lower derivative orders \( 1 \le m \le 8 \). However, it is evident that the agreement between the envelopes and components is systematically improved from \( m = 9 \) to \( m = 14 \) on panels (c) to (h) in Fig. 7, as well as on panels (b)–(h) for m = 15–21 in Fig. 8, to become perfect on panels (b)–(h) in Fig. 9 (\( m = 24, 28, \ldots , 48 \)). This conclusively indicates that the non-parametric dFPT for estimation of the envelopes alone is fully capable of exactly reconstructing all the constituent components with their entire lineshapes. This is what we have set to prove in the present study. Due to the uniqueness of the converged Padé-based reconstructions, coincidence of the envelopes and components from the non-parametric and parametric dFPT, respectively, must necessarily also yield the same quantitative data (complex fundamental frequencies \( \left\{ {\nu_{k} } \right\} \) and amplitudes \( \left\{ {d_{k} } \right\} \) retrieved by these two types of estimations. We have explicitly verified in the present study that indeed the pertinent numerically determined equalities between the fundamental pairs, \( \left\{ {\nu_{k} , d_{k} } \right\}_{{{\text{non}} - {\text{parametric dFPT}}}} = \left\{ {\nu_{k} , d_{k} } \right\}_{\text{parametric dFPT}} \) do hold true. Here, the set \( \left\{ {\nu_{k} , d_{k} } \right\}_{{{\text{non}} - {\text{parametric dFPT}}}} \) has been presently determined by the procedure from Ref. [2]. On the other hand, the gold standard set \( \left\{ {\nu_{k} , d_{k} } \right\}_{\text{parametric dFPT}} \) has been reconstructed by explicitly solving the quantification problem. This verification has been made in the \( {\text{dFPT}}^{\left( - \right)} \) and \( {\text{dFPT}}^{\left( + \right)} \), with the findings \( \left\{ {\nu_{k}^{ - } ,d_{k}^{ - } } \right\} = \left\{ {\nu_{k}^{ + } ,d_{k}^{ + } } \right\} \). These latter two fundamental pairs reconstruct the corresponding input data \( \left\{ {\nu_{k} , d_{k} } \right\} \) from (3.1), so that \( \left\{ {\nu_{k}^{ - } ,d_{k}^{ - } } \right\} = \left\{ {\nu_{k}^{ + } ,d_{k}^{ + } } \right\} = \left\{ {\nu_{k} , d_{k} } \right\} \) as in (2.18a).

4 Clinical relevant of derivative magnetic resonance spectroscopy

In this work, all the expounded features of derivative signal processing are exemplified in the case of magnetic resonance spectroscopy, MRS, for the problem of breast cancer. In this particular problem of major public health concern, clinically reliable quantitative identification of cancer biomarkers is the prerequisite for diagnostics and treatment of patients. Among several diagnostically-informative metabolites, phosphocholine, PC, stands out as a recognized breast cancer biomarker. It is, therefore, of utmost importance to properly visualize and quantify this particular metabolite. Recall that PC is also a recognized biomarker of other malignancies, including ovarian cancer and brain tumors [4, 5]. Our systematic focus on various cancer biomarkers, including PC, and the development of data analytical methods aims to help physicians incorporate molecular imaging by way of MRS into their routine armamentarium. In particular, radiologists need firmly established procedures with clinically reliable information extracted from the tissue scanned by MRS. These conditions translate into the demand for accuracy, stability, robustness and trustworthy cross-validation of all the reconstructions during data analysis by way of signal processing. We have consistently shown that the customary (non-derivative) fast Padé transform, FPT, answers most favorably to this multi-faceted requirement [17, 21, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. Moreover, the FPT has been demonstrated to be of practical use, since its high-resolution capabilities can considerably reduce the scanning time, and this would make MRS more cost-effective, as well as more comfortable for patients [36]. Nevertheless, despite these achievements, there is always room for further strengthening of signal processing.

Our most recent advances [1, 2] have initiated yet another upgrade to break the ground in combating anew the two major stumbling blocks that hamper a deservedly faster entry of MRS into the standard diagnostic protocol for everyday clinical practice. As mentioned, these stumbling blocks are resolution improvement (with a particular emphasis on separating overlapping peaks) and noise suppression. For example, many cancer biomarkers do not show up in spectral envelopes as isolated resonances that are easy to spot and then parametrize. Quite the contrary, many cancer biomarkers, including PC, are hidden in the crowded parts of total shape spectra. The art is to single them out by confidently peering into the substructure of spectral complexes with a number of tightly packed and indiscernible resonances. Moreover, the unavoidable presence of noise additionally masks weak resonances near the background baselines where cancer biomarkers often reside. Faced with these obstacles, physicians are, of course, well aware that mathematical methods are needed to solve such problems.

Yet, physicians are in need of a clear-cut “bottom line” of the mathematical descriptions for any level of the involved sophistication. In other words, what would clinically be most practical is to glean the gestalt of the output of signal processing in a straight-forward, visual, yet quantitative representation with all the relevant information available on a single screen display. Alongside the initial, crowded and noise-corrupted envelope spectrum, the radiologist would optimally like to see the processed data with clearly delineated lineshapes of components for metabolites, particularly those assigned to recognized cancer biomarkers such as PC, superimposed on a maximally flattened background baseline. Additionally, this visualization should be completed with the numbers (chemical shift, \( T_{2k}^{*} \), and concentration for the kth metabolite) above the peaks, especially for the most diagnostically informative molecules. This maximal level of optimization is precisely what the newly-developed non-parametric high-order derivative fast Padé transform, dFPT, offers, as per the presently expounded theory and results sections. These findings fully justify proceeding to applications of the dFPT to encoded data from in vivo MRS as envisaged in our work to be reported soon.

5 Discussion and conclusions

Derivative signal processing by the Padé-based non-parametric, or equivalently, shape estimation alone is presently established as a special and far-reaching quantitative analysis of spectra. The gist of the matter with this unique strategy in the estimation processes is that the quantification problem itself is never set to be solved. This is indicated already in the adjective “non-parametric” for the term estimation, meaning that from the onset, the peak parameters (position, widths, heights) are not envisaged to be found by explicitly solving the spectral analysis problem, which is how the quantification problem is called in the mathematical literature. Yet, at the end of the estimation of the mere envelope lineshapes, the non-parametric derivative fast Padé transform, dFPT, provides the profiles of the underlying components. We show that for high-order \( m \) of the frequency-dependent derivative operator \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \) these components from the non-parametric dFPT coincide exactly with the corresponding components reconstructed by the parametric dFPT (which explicitly solves the quantification problem). With such a finding, we have confirmed that the entire high-order derivative component lineshapes from the envelopes given by the non-parametric dFPT are, in fact, exact.

To make the non-parametric dFPT quantitative, two additional conditions need to be fulfilled, besides predicting the exact component lineshapes. First, the non-parametric dFPT should give the numerical values of the peak signatures (position, width, height) of each component for any \( m > 0 \) (especially for higher derivative orders \( m \) for which the components from exclusively a shape estimation are most accurate). Second, in the most convenient phase-insensitive spectral magnitude mode, the mentioned peak locations, widths and heights for any fixed order \( m > 0 \) in the non-parametric dFPT should uniquely relate to the corresponding peak parameters of the true absorptive profile for the non-derivative \( \left( {m = 0} \right) \) parametric fast Padé transform, FPT. This is verifiable by computing the envelopes in the non-parametric dFPT at the finest grid of the freely-chosen set of sweep frequencies. This permits precise reconstruction of the numerical values of the peak positions, widths and heights of every resonance. Our algorithm simply reads off these parameters from the evaluated derivative envelopes that we now know are, in fact, the components. It is in this straight-forward way that we extract the peak parameters of the derivative lineshapes by the non-parametric dFPT for any order \( m > 1 \). This fulfills the first of the two cited supplementary conditions. To satisfy the remaining, second condition, the parametric dFPT can be used to show that the peak parameters of the components reconstructed by this estimation for any \( m > 1 \) are uniquely related to the pure absorptive components in the non-derivative \( \left( {m = 0} \right) \) parametric FPT. On the other hand, the envelopes of the non-parametric dFPT for high derivative order \( m \) coincide with the components from the parametric dFPT. Such a coincidence implies that the peak parameters of the envelopes with high order \( m \) from the non-parametric dFPT are also explicitly related to the absorptive peak positions, widths and heights of the components retrieved by the non-derivative \( \left( {m = 0} \right) \) parametric FPT. This settles the issue for the second supplementary condition. As such, although being only a shape estimator at the onset of the analysis, the non-parametric dFPT becomes in the end a quantification-equipped estimator of proven validity.

Hence, the non-parametric high-order derivative fast Padé transform is a stand-alone, quantifying processor for robust and accurate estimation of the peak positions, widths and heights of all the physical resonances. It expediently solves the two major problems in signal processing: resolution improvement and noise suppression. It increases the frequency resolution of envelopes to such an unprecedented extent that they fully coincide with their component spectra. At the same time, it is a powerful noise filter with complete removal of all the unphysical information from the processed data. The unique mechanism for solving the two mentioned main problems is in the occurrence that the derivative operator \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \) narrows the peak widths, enhances the peak heights of the physical resonances and diminishes those that are unphysical. This simultaneously improves the resolution and signal-to-noise ratio. This remarkable feature of the derivative transform \( \left( {{\text{d}}/{\text{d}}\nu } \right)^{m} \) runs counter to the integral transform. An integral transform is a smoothing operator: it averages over the fine structural details. A derivative transform does just the opposite: it performs a “differential diagnosis”, so to speak, by disentangling and emphasizing the inner substructures. Nowhere does the occasionally-used alternative name “anti-integrals” for derivatives become more relevant and astoundingly evidenced as in the present application of the dFPT to diagnostics of breast cancer by magnetic resonance spectroscopy, MRS.

Notes

Acknowledgements

This work was supported by King Gustav the 5th Jubilee Fund, The Marsha Rivkin Center for Ovarian Cancer Research and FoUU through Stockholm County Council to which the authors are grateful.

Copyright information

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

This article is published under an open access license.
Please check the 'Copyright Information' section for details of this license and
what re-use is permitted.
If your intended use exceeds what is permitted by the license or if
you are unable to locate the licence and re-use information,
please contact the Rights and Permissions team.