In vivo derivative NMR spectroscopy for simultaneous improvements of resolution and signal-to-noise-ratio: Case study, Glioma

The theme of this study is derivative nuclear magnetic resonance (dNMR) spectroscopy. This versatile methodology of peering into the molecular structure of general matter is common to e.g. analytical chemistry and medical diagnostics. Theoretically, the potential of dNMR is huge and the art is putting it into practice. The implementation of dNMR (be it in vitro or in vivo) is wholly dependent on the manner in which the encoded time signals are analyzed. These acquired data contain the entire information which is, however, opaque in the original time domain. Their frequency-dependent dual representation, a spectrum, can be transparent, provided that the appropriate signal processors are used. In signal processing, there are shape and parameter estimators. The former processors are qualitative as they predict only the forms of the lineshape profiles of spectra. The latter processors are quantitative because they can give the peak parameters (positions, widths, heights, phases). Both estimators can produce total shape spectra or envelopes. Additionally, parameter estimators can yield the component spectra, based on the reconstructed peak quantifiers. In principle, only parameter estimators can solve the quantification problem (harmonic inversion) to determine the structure of the time signal and, hence, the quantitative content of the investigated matter. The derivative fast Fourier transform (dFFT) and the derivative fast Padé transform (dFPT) are the two obvious candidates to employ for dNMR spectroscopy. To make fair comparisons between the dFFT and dFPT, the latter should also be applied as a shape estimator. This is what is done in the present study, using the time signals encoded from a patient with brain tumor (glioma) using a 1.5T clinical scanner. Moreover, within the dFPT itself, the shape estimations are compared to the parameter estimations. The goal of these testings is to see whether, for in vivo dNMR spectroscopy, shape estimations by the dFPT could quantify (without fitting), similarly to parameter estimations. We check this key point in two successive steps. First, we compare the envelopes from the shape and parameter estimations in the dFPT. The second comparison is between the envelopes and components from the shape and parameter estimations, respectively, in the dFPT. This plan for benchmarking shape estimations by the dFPT is challenging both on the level of data acquisition and data analysis. The data acquisition reported here provides encoded time signals of short length, only 512 as compared to 2048, which is customarily employed. Moreover, the encoding echo time was long (272 ms) at which most of resonances assigned to metabolites with shorter spin-spin relaxations are likely to be obliterated from the frequency spectra. Yet, in face of such seemingly insurmountable obstacles, we are looking into the possibility to extract diagnostically relevant information, having particularly in focus the resonances for recognized cancer biomarkers, notably lactate, choline and phosphocholine. Further, we want to see how many of the remaining resonances in the spectra could accurately be identified with clinical reliability as some of them could also be diagnostically relevant. From the mathematical stance, we are here shaking the sharp border between shape and parameter estimators. That border stood around for a long time within nonderivative estimations. However, derivative shape estimations have a chance to tear the border down. Recently, shape estimations by the dFPT have been shown to lead such a trend as this processor could quantify using the time signals encoded from a phantom (a test sample of known content). Further, the present task encounters a number of additional challenges, including a low signal-to-noise ratio (SNR) and, of course, the unknown content of the scanned tissue. Nevertheless, we are determined to find out whether the nonparametric dFPT can deliver the unique quantification-equipped shape estimation and, thus, live up to the expectation of derivative processing: a long-sought simultaneous improvement of resolution and SNR. In every facet of in vivo dNMR, we found that shape estimations by the dFPT has successfully passed the outlined most stringent tests. It begins with transforming itself to a parameter estimator (already with the 3rd and 4th derivatives). It ends with reconstructing some 54 well-isolated resonances. These include the peaks assigned to recognized cancer biomarkers. In particular, a clear separation of choline from phosphocholine is evidenced for the first time by reliance upon the dFPT with its shape estimations alone.

challenging both on the level of data acquisition and data analysis. The data acquisition reported here provides encoded time signals of short length, only 512 as compared to 2048, which is customarily employed. Moreover, the encoding echo time was long (272 ms) at which most of resonances assigned to metabolites with shorter spin-spin relaxations are likely to be obliterated from the frequency spectra. Yet, in face of such seemingly insurmountable obstacles, we are looking into the possibility to extract diagnostically relevant information, having particularly in focus the resonances for recognized cancer biomarkers, notably lactate, choline and phosphocholine. Further, we want to see how many of the remaining resonances in the spectra could accurately be identified with clinical reliability as some of them could also be diagnostically relevant. From the mathematical stance, we are here shaking the sharp border between shape and parameter estimators. That border stood around for a long time within nonderivative estimations. However, derivative shape estimations have a chance to tear the border down. Recently, shape estimations by the dFPT have been shown to lead such a trend as this processor could quantify using the time signals encoded from a phantom (a test sample of known content). Further, the present task encounters a number of additional challenges, including a low signal-to-noise ratio (SNR) and, of course, the unknown content of the scanned tissue. Nevertheless, we are determined to find out whether the nonparametric dFPT can deliver the unique quantification-equipped shape estimation and, thus, live up to the expectation of derivative processing: a long-sought simultaneous improvement of resolution and SNR. In every facet of in vivo dNMR, we found that shape estimations by the dFPT has successfully passed the outlined most stringent tests. It begins with transforming itself to a parameter estimator (already with the 3rd and 4th derivatives). It ends with reconstructing some 54 well-isolated resonances. These include the peaks assigned to recognized cancer biomarkers. In particular, a clear separation of choline from phosphocholine is evidenced for the first time by reliance upon the dFPT with its shape estimations alone.

Mathematical modeling of complicated systems and the notion of critical parameters
Studying large and complicated systems in chemistry, physics, biology and beyond need not necessarily require overly complicated mathematics. What is essential, however, is to reconstruct the so-named "critical parameters" capable of capturing the main features of the dynamics of the examined system. It is mathematical modeling that can extract these critical parameters from generic systems. Such parameters are called "critical" (or "fundamental" in signal processing or "characteristic" or "eigen-parameters" in quantum physics) because they can reveal or unfold the unknown inner fabric of the examined system. As such they are of critical importance in the adequate and reliable descriptions, characterizations as well as interpretations of measured data, including time signals in nuclear magnetic resonance (NMR) spectroscopy. Importantly, critical parameters quantify the experimentally recorded information, thus making it fully usable in practice for decision making (like medical diagnoses using in vitro and in vivo NMR spectroscopy).
To provide the set of critical parameters, mathematical models rely on certain simplifying assumptions. This is actually the first meeting point of theory and measurements. The reason is that, in fact, every measurement also resorts to some simplifying assumptions. Simplifying assumptions set up a working hypothesis for one or more underlying mechanisms that are assumed to govern the dynamics of the system. Remarkably, a dominant feature of modern applied mathematics is that, most frequently, merely a few assumed simple mechanisms may suffice for accurately describing and interpreting highly complex behavioral patterns observed in measurements.
An excellent illustration of this latter feature of mathematics in e.g. physics, chemistry and biomedicine is NMR spectroscopy . In NMR spectroscopy, the dynamical behavior of a complicated system is investigated by measuring time signals, or equivalently, free induction decay (FID) curves. These FIDs are the answer of the given system to an external perturbation caused by three fields from the electromagnetic spectrum (static as well as gradient magnetic fields and radio-frequency pulses). The goal of NMR spectroscopy is to uncover the hidden structure of molecules from the examined sample, excited externally. This powerful and versatile method of elemental analysis of matter on all levels, including human tissue, can identify and quantify the intrinsic constituents or components of molecules and chemical compounds.
It is a non-invasive, i.e. non-destructive procedure after which the sample remains intact. Advantageously, electromagnetic fields in NMR spectroscopy belong to the category of non-ionizing radiations. This is due to the fact that the applied electromagnetic fields cannot disrupt the internal structure of cellular molecules. Specifically, energies imparted to molecules are far below their excitation and ionization thresholds. Despite the word "nuclear" in NMR spectroscopy, this method is not within the realm of nuclear medicine. To avoid such a potential association (and the unjustified patient's fear from the actually nonexistent nuclear radiation in NMR), the alternative name, magnetic resonance spectroscopy (MRS), is used in medical diagnostics. The main molecules of interest to MRS diagnostics are called metabolites. These are molecules that participate to the various pathways of cellular metabolism.
When plotted, each measured FID appears as a tightly packed set of exponentially damped oscillation waveforms. The most intense oscillations are at the beginning of data acquisition. The weakest oscillations are at the FID's tail near the end of the total acquisition time T where noise prevails. Such oscillations are reminiscent of attenuated sinusoidal and cosinusoidal trigonometric functions. This observation leads straight to the main simplifying assumption as a working hypothesis in NMR spectroscopy: the introduction of a mathematical model, which describes the given encoded FID as a sum of complex damped exponential functions multiplied by their complex amplitudes.
The appearance of the sinusoidal and cosinusoidal trigonometric response functions to the applied electromagnetic fields is physically plausible and, in fact, expected. The reason is that the electromagnetic field itself is comprised of two unattenuated trigonometric waveforms (sine, cosine). These are the two component fields (electric, magnetic) of the electromagnetic field. Such two components oscillate in two orthogonal planes in accordance with the ∕2 phase difference between the sine and cosine functions.
Pure sinusoidal sin (2 f k t) and cosinusoidal cos (2 f k t) oscillations with real linear resonance frequencies f k are not attenuated as their maxima and minima keep on periodically appearing and are always bound by the values ±1. Conveniently, the Euler formula combines these latter two trigonometric functions into a complex undamped exponential, e 2 if k t = cos (2 f k t) + i sin (2 f k t). The latter expression to paraphrase Feynman is 'our jewel as the most beautiful mathematical formula'. However, in every realistic application of transient phenomena, these oscillations are attenuated when used to describe a time evolution of natural systems.
Attenuation occurs because, after a certain finite period, an evolution process decays to zero for every naturally occurring non-stationary phenomenon. According to quantum mechanics, the time evolution operator U is an exponential operator, U(t) = e −iΩt , where Ω is the generator of the system's dynamics (a Hamiltonian operator in physics as a sum of the kinetic and interaction energy operators). As a result, the decay probability of the k th state of a system follows the exponential law, e − k t , where k > 0 is the transition rate. The same holds true also for relaxation of the bulk magnetization of a sample in every NMR measurement.
Here "relaxation" means return of the bulk magnetization to the equilibrium after excitation of a sample by radio-frequency pulses. As such, encoded times signals are exponentially attenuated. This is how a complex damped exponential (called complex harmonic) makes its foothold in the NMR phenomenon, e 2 if k t− k t ≡ e i k t . The complex angular k and linear k resonance frequencies are connected by the relation k = 2 k , where the real Re( k ) and imaginary Im( k ) parts of k = Re( k ) + iIm( k ) are Re( k ) = f k and Im( k ) = k ∕ (2 ). In this way, one is already accustomed with the two real-valued critical parameters in NMR spectroscopy, the resonance or fundamental frequency (chemical shift) f k and the relaxation time, 1∕ k . The time constant T ⋆ 2k = 1∕ k is the lifetime of the metastable resonant state of the spin system. The main contributor to T ⋆ 2k is the dipolar interaction of the spin-active nuclei (nuclei with non-zero spins). This interaction (J-coupling of nuclear spins) causes resonance peaks in a frequency spectrum to split into their multiplets, an occurrence which is transparent in the specialized name, NMR J-spectroscopy. Resonant frequencies (frequencies at which the system resonates with an external perturbation) vary from one molecule to another. 1 It is dependent on a local shielding of the static magnetic field by the electronic cloud (in the host atoms/molecules) surrounding the 1 3 Journal of Mathematical Chemistry (2021) 59:2133-2178 spin-active nuclei. Thus, reconstructing k offers a chance to unfold the molecular structure.

Parameter estimators
The mentioned mathematical model for a time signal also follows directly from the quantum-mechanical nature of the NMR effect without the need to visually inspect any FID plot. Namely, an FID given by a linear combination of complex attenuated harmonics, would be the result of an idealized measurement with no corruption (no systematic errors, no noise, etc.). This would be in full harmony with the quantummechanical auto-correlation functions C(t) that are sums of complex damped harmonics. Such auto-correlation functions are the matrix elements (Φ 0 |U|Φ 0 ) of the quantum-mechanical time evolution operator U(t) = e −iΩt , where Φ 0 is the ground state of a dissipative system.
In practice, no encoding could be ideal (error-free, noiseless, … ) since various experimental imperfections always introduce disturbances in measured time signals causing them to deviate from the sums of complex damped harmonics. Nevertheless, the underlying model is physically as well as mathematically plausible and adequate and, therefore, it can be viewed as an appropriate starter. To handle measuring imperfections (noise, etc.), this starter needs to be refined by some robust computational algorithms.
Attractive as it surely is, an NMR spectrometer/scanner, being a measuring instrument, cannot quantify the encoded time signals on its own. Quantification is carried out by mathematical methods. The system (scanned sample) is parametrized through quantification. Quantification of the encoded NMR data consists of retrieving the sought parameters, the resonant fundamental frequencies from the damped harmonics whose multiplying pre-factors are the amplitudes d k = |d k |e i k of each component harmonic of the given FID. Since both the resonant frequency k and the associated amplitude d k are complex numbers, each harmonic is described by four real-valued parameters, two for k = Re( k ) + iIm( k ) ≡ {Re( k ), Im( k )} and two for d k = Re(d k ) + iIm(d k ) ≡ {Re(d k ), Im(d k )}. Thus, the critical NMR parameters can here be specified as: {real resonant frequency Re( k ) ≡ f k , imaginary resonant frequency Im( k ) ≡ k ∕(2 ), magnitude (modulus |d k | ) of amplitude and phase k of amplitude} per harmonic. An extraction of these critical parameters from the encoded FID would uncover the system's dynamics. This is called the quantification problem in the field of signal processing (in mathematics, it is known as the spectral analysis problem).
Real resonant frequencies f k , as the eigen-frequencies (or characteristic frequencies) of the given molecule, inform about the electronic environment of resonating nuclei. They are equal to the discrete energies (in the atomic or natural units) at which the system actually exists. The reciprocals of imaginary frequencies k are the lifetimes of non-stationary states of the system. These decaying or meta-stable states of a dissipative system are driven or governed by a non-Hermitean generator Ω of the dynamics ( Ω † ≠ Ω). The magnitudes of complex amplitudes are proportional to the number or the abundance (concentration) of resonating nuclei. The phase of a complex amplitude represents the phase of the sinusoidal/cosinusoidal component in the FID.
Ideally, all these phases should be equal to zero in the absence of various measuring errors e.g. a time delay (dead time) as the difference between the end of the excitation pulse and the beginning of data acquisition (encoding). Overall, we see that the physical basis of the NMR phenomenon itself dictates the most natural mathematical model for the waveform of an encoded time signal as a linear combination of attenuated complex harmonics. Generally, there could be no trustworthy data analysis/interpretation in NMR spectroscopy without a solid quantum physics platform on which to build the robust algorithms for explicit numerical computations. This is how quantum mechanics sees the field of signal processing.
Moreover, a time signal as a sum of complex damped exponentials is not limited to NMR spectroscopy alone. It is universally applicable to time signals in vastly different research branches. The basis of this versatility is the equivalence of time signals with quantum-mechanical auto-correlation functions. The universality of the discussed mathematical model for general time signals emerging from measurements is rooted in the circumstance that the entire information obtainable from measurements is predictable by the Schrödinger equation of quantum physics. This is the content of the well-known quantum-mechanical "completeness relation" via the exact decomposition of the unity operator (associated with maximal probability for an outcome of a measurement) in terms of the eigen-projection operators built from the complete set of the Schrödinger eigen-wave functions [1,3].
It is then quantum-mechanical signal processing that can achieve the exact parametrization of general systems studied by NMR spectroscopy. This pathway consists of solving the eigenvalue problem with the Schrödinger time evolution operator U. Diagonalization of U in e.g. the basis of the Schrödinger time-dependent wave functions Φ(t) is possible without the explicit knowledge of the generator Ω of the system. All that is needed for diagonalization of any operator are its matrix elements. Paradoxically, the matrix elements of an unknown operator can be computed. Here too, the answer is in the equivalence of the auto-correlation functions {C n } and time signals {c n }. It is this equivalence that makes it possible to explicitly express all the matrix elements of U by means of the given input time signal data points, {c n } (0 ≤ n ≤ N − 1), without ever knowing the operator Ω itself. Here, N is the total signal length, which is connected to the total acquisition time T by T = N where is the sampling or dwell time, which is the reciprocal of the bandwidth (BW) or sweepwidth (SW).
Technically, the discussed eigen-problem is, in fact, a generalized eigen-problem due to non-orthogonality of the discretized Schrödinger basis set functions {Φ n } . The solutions of this quantification problem are the eigen-frequencies and the stationary eigen-state wave functions Ψ k . Here, the obtained eigen-frequencies are the sought resonant complex frequencies { k }. On the other hand, the squared projections of the stationary Schrödinger eigen-state wave functions Ψ k on to the ground state Φ 0 of the system are the complex amplitudes {d k } of the components of the time signal. This is the essence of quantification by quantum-mechanical signal processing, where the critical parameters are the eigen-parameters (eigenfrequencies k and eigen-amplitudes d k ). This particular kind of signal processing involves simultaneously both the stationary (time-independent) and non-stationary (time-dependent) Schrödinger equations. Such an estimation does not necessarily use the spectrum. However, the spectrum can be generated afterward from the obtained critical parameters [1,23].
An FID and a spectrum are the two sides of the same coin. Viewed together, they constitute a dual representation: the auto-correlation function in the time domain is associated with a spectrum in the frequency domain. The two representations are connected by the direct and inverse Fourier integrals. The role played by the evolution operator U(t) = e −iΩt in the time domain is taken up by the resolvent or Green operator G( ) = 1∕( − Ω) in the frequency domain. The operators U(t) and G( ) are also inter-connected by the direct and inverse Fourier integrals in terms of the time t and frequency variables, respectively [1,23].
In general, while the structure of the system may be totally opaque (invisible) in an FID, it could become at least partially transparent (visible) in the frequency domain by some adequate, mechanistically grounded mathematical models, the prime example of which is the universal framework of quantum-mechanical signal processing [1]. This forms the basis of the so-called parameter estimators in the field of signal processing. In practical applications, for super-resolution and the highest accuracy in robust and self-correcting estimations of the critical parameters, the fast Padé transform (FPT) holds promise of qualifying as the method of choice not only for NMR spectroscopy, but also in the whole area of signal processing. This processor gives the envelope (total shape spectrum) as the unique quotient of two polynomials extracted from the given time signal. Such a rational response function (polynomial ratio) can be computed at any sweep frequency.
The mentioned spectral transparency is manifested by the emergence of a number of peaks or resonances in the given graphed frequency response function. Such a frequency representation contains the same parametrization of the system as encountered in the associated FID. Here, however, we are talking about the peak parameters. To understand their meaning, it is useful to relate them to the eigencharacteristics of the components of the already analyzed FID. Thus, the real and imaginary parts of the complex eigen-frequency k in the given harmonic from the studied FID are directly connected to the peak position and the peak width, respectively. In a spectrum, the full width at half maximum (FWHM) is usually used while referring to the peak width. The magnitude |d k | of the complex amplitude d k of a harmonic in the FID is the doubled peak area of a purely absorptive Lorentzian spectrum.
A Lorentzian is a perfectly symmetric spectral bell-shaped line profile (the Cauchy distribution in mathematics) when the phase of the given component of the FID is zero. For nonzero phases, as is always the case with encoded FIDs, the spectral lineshapes are asymmetrical and deformed. In such a case, both the real and imaginary parts of a complex envelope contain mixtures of absorptive and dispersive lineshapes. This causes e.g. the real part of a complex envelope to take on both positive and negative values. To somewhat mitigate this lack of positive-definiteness, usually the encoded FID is multiplied by a factor of the form e i 0 , where 0 is a zero-order phase correction. Of course, a single exponential e i 0 , cannot make the entire envelope positive-definite. Although the ensuing phased lineshape profile might appear as being more symmetric and less deformed, it is most frequently unamenable to quantification.
Overall, the stated dual representation translates the inquiry of the critical parameter structure of the time signal to the quest for the true components of the corresponding frequency envelope spectrum. This is because the genuine components of the given envelope are generated exclusively from the same critical parameters that are the building blocks of the corresponding time signal. Consequently, it is immaterial in which analysis domain (time, frequency) the quantification problem is solved since the results of a proper implementation ought to be the same. In fact, as explained, to solve the quantification problem exactly, quantum-mechanical signal processing uses simultaneously both the time and frequency domain observables by way of intertwining the non-stationary and stationary Schrödinger equations [1,23]. This is to be contrasted to all the available fitting techniques that artificially separate the quantification problem in the time and frequency domains [24,25]. Invariably, however, the latter obtained estimations in the two equivalent domains are unequal as they stem from independent and unrelated fittings of the given time signal and the associated frequency envelope spectrum.

Shape estimators
In the preceding sub-section, we discussed parameter estimation by the parametric version of the FPT. The nonparametric variant of the FPT can perform shape estimation as soon as the frequency-dependent rational polynomial P K ∕Q K becomes reconstructed from the given time signal. This FPT, as a shape estimator, cannot quantify on its own. Shape estimators deal with reconstructions of spectral lineshapes alone with no autonomous capability for solving the quantification problem (determining the peak position, width, height and phase of each physical resonance). However, a purely shape estimation by the derivative fast Padé transform (dFPT) can autonomously perform quantification. The dFPT consists of applying the differential operator of the m th order, D m = (d∕d ) m , to the complex spectrum P K ∕Q K . The principal advantages of the dFPT (derivative) over FPT (standard, nonderivative) is simultaneous improvements of resolution and signal-to-noise ratio (SNR).
The mechanism for this twofold achievement is in the linewidth narrowing (coupled with augmentation of the associated peak heights) and flattening of the background baseline. Increased peak heights themselves indicate improved SNR. These theoretically grounded concepts of derivative signal processing have been implemented in NMR spectroscopy for using different types of time signals: synthesized (noise-free, noise-corrupted) [16][17][18] and encoded both from phantoms [21][22][23] and patients [20]. The obtained results were astounding as they put into practice all the mentioned facets of the dFPT by demonstrating simultaneous improvements of resolution and SNR.
The present study applies the FPT and dFPT to time signals encoded from the brain of a patient with glioma. The goal is to determine how much of all the mentioned capabilities of the dFPT can be confirmed in the case of a very demanding short time signal (only 0.5 KB) encoded at a clinical scanner of a low magnetic field strength (1.5 T) with a partial suppression of water in the course of measurement. The clinical goal is to find out whether cancer biomarkers can be reliably identified and quantified (e.g. separation of phosphocholine from choline), a task which has never been accomplished before by any shape estimator for encoded FIDs.
Note that, throughout, we employ the so-called 'natural units' (also called 'atomic units') with ℏ ≡ h∕(2 ) = 1, c = 1, e = 1 and m e = 1 , where h is Planck's constant, c is the speed of light in a vacuum, e is the electron charge and m e is the electron mass [26]. In these units, energy (E) and frequency ( ) are synonyms, because for ℏ = 1, the defining relation E ≡ ℏ becomes E = . Moreover, in the natural units, every physical quantity is either dimensionless (like e and m e ) or its dimension is expressed in some powers (integers or rational numbers) of the unit of length (cm), e.g.
[E] = cm −1 . In spectroscopies (be they in physics, chemistry, etc.), energy is often given in units cm −1 [27]. To return to the actual dimension of a given quantity, its dimension in the natural units should be multiplied by the requisite powers of ℏ and c.

Quantum-mechanical solution of the quantification problem
In the Introduction, we emphasized that the time evolution of a general system is described quantum-mechanically by the linear evolution operator in the form of a complex exponential U(t) = e −iΩt . Moreover, as stated therein, this automatically implies that the time-dependent auto-correlation function C(t) ≡ (Φ 0 |U|Φ 0 ) of that system is given by a linear combination of the complex harmonics: where C(n ) ≡ C n and 0 ≤ n ≤ N − 1. Moreover, we have Im( k ) < 0 to secure an exponential decay of C(t) to zero at t → ∞. Here, is the sampling time for the discretized version t n of the continuous time variable t, with t n = n . This parametrization of the state of the system is given in terms of the scalar quantities { k , d k }, the fundamental or characteristic resonance angular frequencies and the corresponding amplitudes (both complex), respectively.
Quantum mechanics finds the stationary Ψ k and nonstationary Φ(t) states of any system (including a human tissue) by solving the respective time-independent and time-dependent Schrödinger equations: Recall that we are using the mentioned natural units and, therefore, the reduced Planck's constant ℏ is set to unity in (2.1), (2.2) and elsewhere. Throughout, we assume that all the eigen-values { k } of Ω are discrete. The boundary condition at the onset of time evolution (t = 0) to the second of these two equations is Φ(0) = Φ 0 , where Φ 0 is the known initial state of the system. For a dissipative system, the dynamic operator Ω is time-independent. This circumstance solves the time-dependent Schrödinger equation as: The main quantum-mechanical working postulate is that the entire information about the system is contained in the state wave functions. This is expressed as the completeness relation: where 1 and k are the unity and the projection operator, respectively. This equation is the spectral decomposition of the unity operator, 1 . Rewriting U as the identity U(t) = U(t)1 and using therein the stationary Schrödinger equation from (2.2) and the completeness relation (2.4) would yield the spectral decomposition of the evolution operator: Here, it is seen that the projection operator k has a clear physical meaning as it represents the operator-valued amplitude, the strength of the harmonic e −i k . In a similar manner, it follows that the exact wave functions Ψ k and Φ(t) also have their spectral decompositions, both in terms of the fundamental parameters { k , d k } ∶ In deriving Eq. (2.8), we used the definition (2.4) for the projection operator k . The complex amplitude d k from Eq. (2.8) is seen as the matrix element of the projection operator k taken over Φ 0 . In Eq. (2.8), this matrix element, (Φ 0 | k |Φ 0 ), is further reduced to the squared inner (or scalar) product of Ψ k and Φ 0 via (Ψ k |Φ 0 ) 2 . Such a scalar product is symmetric in the sense that it has no usual complex conjugation, i.e. (a|b) = (b|a), to reflect the lack of hermiticity of the system's 'Hamiltonian' operator, Ω † ≠ Ω.
The quantum-mechanical solutions (2.6) and (2.7) are formal because the dynamic operator Ω of the system is unknown in any inverse problem, including the harmonic inversion, i.e. the quantification problem in NMR spectroscopy. However, (2.6) this would be of no concern if the matrix elements of Ω or U were given. They are. In particular, all the needed matrix elements of U( ) in the Schrödinger basis set function {Φ n } (0 ≤ n ≤ M − 1) (also called the Krylov basis in the field of signal processing) are explicitly expressed in terms of the input data, the known auto-correlation functions: where, as before, N is the total length of the sequence of the auto-correlation functions {C n } (0 ≤ n ≤ N − 1). The only function f (x) having the property f (x)f (y) = f (x + y) is the exponential function (scalar, operator). This feature implies that even the s th power of U is also reduced to an auto-correlation function: The special cases with s = 0 and s = 1 give respectively the overlap matrix S m,n = (Φ m |Φ n ) and the matrix elements U n,m = (Φ m |U( )|Φ n ) of the time evolution matrix = {U n,m }. The basis set {Φ n } is non-orthogonal as the overlap matrix = {S n,m } is not zero. This means that in quantum mechanics for harmonic inversion with the Schrödinger expansion functions, the fundamental parameters are not obtained by solving the ordinary eigen-value problem: Instead, the corresponding generalized eigen-value problem with the presence of the overlap matrix needs to be solved: Here, the elements {A n,k } of the column matrix k = {A n,k } are the expansion coefficients from the development of Ψ k in the Schrödinger expansion basis set, Once the set of the solutions {z −1 k , A n,k } becomes available by solving (2.12), the fundamental frequencies and amplitude follow as:

Quantum-mechanical frequency-domain response function and the fast Padé transform
The quantum-mechanical dual representation (time, frequency) of the system's state is provided by the Fourier integrals. The Fourier integrals map e.g. the timedependent evolution operator U(t) to the frequency-dependent Green operator, G( ), and vice verse (direct and inverse Fourier integrals). The Green operator G( ) = ( − Ω) −1 , is also known as the resolvent operator which is, by definition, The role of the auto-correlation function C(t) in the time domain is taken up by the Green function G( ) in the frequency domain. The Green function is a frequency-dependent spectrum, a scalar quantity given by the matrix element of operator G( ) from Eq. (2.16) over the two copies of initial state Φ 0 of the system: Here, the term d k ∕( − k ), as the k th partial fraction, represents the k th component spectrum, which is a scalar version of Padé approximant (PA) in its paradiagonal form: is called the Heaviside partial fraction decomposition of the spectrum G( ). Thus, the sum of all the K components g k ( ) represents the envelope spectrum, G( ) = ∑ K k=1 g k ( ). The sum of the PAs for the components g k ( ) in ∑ K k=1 g k ( ) is also a PA, a rational polynomial, a ratio of two polynomials. In this latter rational polynomial, the degrees of the numerator and denominator polynomials are K − 1 and K, respectively, P K−1 ( )∕Q K ( ). Hence, if the summation over the K components g k in the Heaviside partial fraction representation of the envelope G( ) is explicitly carried out, the following paradiagonal PA would be obtained: Journal of Mathematical Chemistry (2021) 59:  This means that all the processors, in fact, compute the Padé approximant, if they use the rhs of Eq. (2.18) to express the envelope spectrum by the Heaviside sum of partial fractions. Moreover, this Heaviside sum must be obtained if the auto-correlation function C(t), or equivalently, the time signal c(t), is represented by the sum (2.1) of complex damped exponentials, as is the ubiquitous practice throughout NMR spectroscopy and beyond. In signal processing, the alternative name for the Padé approximant is the fast Padé transform, FPT [1].

Probabilistic interpretation of the auto-correlation function and ergodic hypothesis
The name auto-correlation function C(t) comes from connecting the same wave function Φ(t � ) of the system at two different times, t � = 0 (initial) and t � = t (a subsequent time instance). In the discretized representation, this is written as If we now insert the system's discrete wave function Φ n from (2.7) into this latter definition of C n , we would have: where d 1∕2 (Φ 0 |Ψ k ) = d k as in (2.8). This is the decomposition of the auto-correlation function in terms of the eigen-frequencies and eigen-amplitudes { k , d k }, in accord with Eq. (2.1). As such, the quantity C n from Eq. (2.1) is not a mathematical model introduced ad hoc for describing the auto-correlation function. Rather, it is derived from quantum physics with no approximation. Alternatively, as detailed in Ref. [23], we can derive the same result (2.1) by expressing C n as a contour integral over the system's frequency response function, G( ). As per Eq. (2.18), this latter function, as the Green's function, or equivalently, the resolvent function, has all its singularities as poles located at the eigen-frequencies { k } of the dynamic operator Ω. Therefore, the mentioned contour integral can be computed exactly in terms of the Cauchy residues at these poles. The result is again C n = ∑ K k=1 d k e −i k n as in (2.1) or (2.21). The physical picture conveyed by the form (2.1) of the auto-correlation function C n is plausible [28]. This function describes the system as distributing itself among the eigen-states Ψ k of the dynamic operator U( ), with the weight factors given by the squared overlap (Φ 0 |Ψ k ) 2 between the initial state Φ 0 and the eigen-state Ψ k . That overlap, as per Eq. (2.8), is the amplitude d k or the strength (intensity) of the transition from Φ 0 to Ψ k .
Since the auto-correlation function C(t) from Eq. (2.1) is a probability amplitude function, in order to be a physical quantity, it must tend to zero as t → ∞. This indeed occurs due to Eq. (2.1), where lim t→∞ e −iRe( k )t+Im( k )t = 0 because Im( k ) < 0. Such a situation amply illustrates the importance of complex valuedness of the fundamental frequencies { k }. In complex characteristic frequencies { k }, the condition Im( k ) < 0 secures the physical role the complex harmonic e −i k t in Eq. (2.1). This harmonic is not an ordinary, i.e. unattenuated oscillating complex exponential, which has no limit at t → ∞. Rather, it is a damped complex exponential, This latter feature gives the physical meaning of a proper probability amplitude function for the auto-correlation function C(t). It guarantees that, for Im( k ) < 0, a true transition out of the initial state Φ 0 can occur in the sense that C(t) → 0 as t → ∞. At a fixed finite time t < ∞, the non-zero probability of such a transition is provided by the squared magnitude |C(t)| 2 of the complex probability amplitude There is yet another notable property of every datum C n in the assembly {C n } (0 ≤ n ≤ N − 1) as dictated by the geometric progression form (2.1) itself. For each number n, this form is a linear combination of the weighted harmon- . We then see that the collection of such K weighted harmonics contains the entire set of the fundamental solutions { k , Ψ k } (1 ≤ k ≤ K) of the quantification problem (2.2), the Schrödinger eigenvalue problem. In other words, every auto-correlation function point C n from Eq. (2.1) is a representative of the entire system in the sense of encapsulating all the critical parameters { k , d k } of the system. 2 However, even though every auto-correlation function point holds the whole fun- , the impact of each datum C n on shaping the system's time-domain response (2.1) to an external perturbation is not the same. The reason is in the nature of the geometric sequence (2.1), where each harmonic e −i k is raised to a power of n to become e −i k n = e −iRe( k )n +Im( k )n . For Im( k ) < 0, this implies more exponential damping with the passage of time n , meaning that the earlier sampled data points C n will play a more significant role than those toward the tail of the auto-correlation function.

From theory of quantum physics to the field of signal processing
To pass from quantum physics to the field of signal processing, it suffices to note that the auto-correlation functions C(t) and C n from (2.1) are equivalent to the time signals c(t) and c n , respectively [1]. This occurs because the time signal c(t) is also given by a linear combination of damped complex exponentials, and so is c n ∶ The only difference is in using two different harmonics e −i k t and e i k t with the opposite restrictions Im( k ) < 0 and Im( k ) > 0 for the auto-correlation function (2.1) and the time signal (2.22), respectively. This is immaterial, as it is only a convention. Due to the equivalence between C n and c n , the quantum-mechanical probabilistic interpretation of the auto-correlation function automatically holds true also for the time signal and, thus, need not be repeated.
After the just expounded realization that the best physics theory, quantum mechanics, predicts the exact spectral envelope in the form of the Padé approximant, P K−1 ( )∕Q K ( ), there is no need to reconstruct the fundamental frequencies and amplitudes, { k , d k }, by solving the generalized eigen-value problem (2.12). The alternative is to use directly one of the algorithms from the FPT. In principle, the FPT can be employed as the FPT (−) and FPT (+) in terms of the independent variable z −1 and z, respectively. Two different analytical properties of the FPT (+) and FPT (−) , inside ( |z| < 1 ) and outside ( |z| > 1 ) the unit circle, respectively, in the complex z−plane influence notably their overall performances.
In the past, we examined their performances by reporting detailed illustrations [4][5][6][7][8][9][10][11][12][13]. In our earlier investigation from Ref. [20], the FPT (+) has been applied to a time signal encoded from a pediatric patient with cerebral asphyxia. In the current work, the FPT (−) will be applied and denoted simply as the FPT. Likewise, for brevity, no minus superscript will be used in any of the related quantities (e.g. the numerator/denominator polynomials, their expansion coefficients, reconstructed frequencies/amplitudes, Padé spectra, etc.). Moreover, from now on, only the diagonal FPT will be utilized, Typically, the FPT begins by reconstructing the envelope spectrum, e.g. in its diagonal form P K (z −1 )∕Q K (z −1 ) using only the intact input time signal {c n } (0 ≤ n ≤ N − 1) without finding first the fundamental resonance parameters { k , d k }. In other words, initially, the FPT does not compute the envelope P K ∕Q K by the Heaviside partial fraction summation of the type (2.20), as the critical parameters { k , d k } are still unknown at this stage. Rather, the FPT offers two pathways, one nonparametric, and the other parametric.
The nonparametric FPT completes its task as soon as the expansion coefficients {p r , q s } of the numerator and denominator polynomials {P K , Q K }, respectively, have been extracted. This is sufficient to compute the sought envelope P K (z −1 )∕Q K (z −1 ) at any desired mesh of the sweep frequencies f = Re( ), which is entrenched in the harmonic z −1 = e −2 i . Therefore, this version of the FPT represents a shape estimator as it predicts only the forms of the lineshape profiles in the chosen mode (magnitude, power, real or imaginary part, Argand plot, … ). It does not autonomously generate the resonance peak parameters { k , d k }. Algorithmically, the expansion coefficients {p r , q s } are found from the definition of the nonparametric FPT (i.e. from the difference, Input-Model = Error): where the lhs is the exact spectrum in the form of the MacLaurin polynomial, the finite z-transform [1]. The MacLaurin polynomial ∑ N−1 n=0 c n z −n itself is the truncated MacLaurin series expansion for the exact Green function. The error O(z −2K−1 ) in the matching condition (2.23) is a series in powers of z −1 beginning with z −2K−1 . To proceed, the expansion representations of {P K , Q K } are used: We omit O(z −2K−1 ) and multiply Eq. (2.23) by Q K (z −1 ) to write the relation . In other words, from the mathematical standpoint, the nonparametric FPT is a completely algebraic processor. The parametric FPT needs only one numerical step beyond the nonparametric FPT to solve the quantification problem, i.e. to determine the fundamental eigen-set { k , d k } and, thus, to reconstruct the resonance parameters (peak position, widths, heights, phases). The concluding step in the parametric FPT is rooting the denominator polynomial, Q K (z −1 ) = 0. This is the so-called characteristic or secular equation whose K roots give the fundamental frequencies k = (i∕ )ln(z −1 k ), as in Eq. (2.13). A polynomial rooting is a nonlinear numerical operation. However, even this last nonlinear computational step in the FPT is avoided altogether by diagonalizing the equivalent Hessenberg matrix (an extremely sparse matrix with the expansion coefficients {q s } on the first row and zeros elsewhere) [1]. Finally, the amplitudes {d k } are obtained from the analytical expression for the Cauchy residue at z −1 = z −1 k of the envelope spectrum P K (z −1 )∕Q K (z −1 ), which is available from the nonparametric FPT. The result is given by: With the availability of the fundamental parameters { k , d k }, we can construct the envelope also in the parametric FPT. The result is the Heaviside partial fraction representation for the envelope, which for the parametric diagonal FPT, takes the form: On the rhs of this equation, the k th term d k z −1 ∕(z −1 − z −1 k ), lying under the summation sign, represents the k th component spectrum, the complex-valued lineshape profile of the k th resonance.
As to the spectral modes, only the magnitude |P K ∕Q K | and the real part Re(P K ∕Q K ) of the complex envelope P K ∕Q K will be presented. For brevity, these will written as: When plotting both an envelope (nonparametric FPT) and the corresponding components (parametric FPT) on the same graph, for their clearer distinction, we shall use the subscripts 'Tot' and 'Comp', respectively: We see how the nonparametric and parametric FPT work in concert. First, we reconstruct the envelope P K (z −1 )∕Q K (z −1 ) by using the nonparametric FPT. Then from this point and on, we pass to the parametric FPT. The above outlines show that the parametric FPT is ideally suitable for solving the quantification problem since only a minimal computational effort of entire linear algebra is necessary to exactly determine all the physical resonance parameters. These are (i) the complex fundamental frequencies k that give the peak positions Re( k ) = f k as well as the peak widths Im( k ) = k and (ii) the corresponding complex fundamental amplitudes yielding the magnitudes |d k | (proportional to the peak areas) as well as phases k = arg(d k ).
The derivative fast Padé transform, dFPT, is of primary interest when the nonparametric version is used for shape estimation. When the derivative operator of order m, defined as D m = (d∕d ) m , is applied to the nonparametric envelope P K (z −1 )∕Q K (z −1 ), the result is either a general analytical formula given by the Bell polynomial [29,30] or a stable recursion [22]. In practice, earlier and presently, we employ the recursive algorithm. In the case of the parametric dFPT, also the analytical expressions are available for both components and envelopes [22]. For instance, the envelope in the parametric dFPT is reduced to the Eulerian polynomial [22,31]. The convention (2.31) for the nonderivative FPT, will similarly be used also for the derivative dFPT (to be presented in the magnitude mode alone) as follows: In the Sect. 3, the obtained envelopes from the nonparametric dFPT are benchmarked in comparisons with the parametric dFPT. Such comparisons are made on two different levels. First, we compare the envelopes from the nonparametric and parametric dFPT. Subsequently, we compare the envelope from the nonparametric dFPT with the components from the parametric dFPT. This is ultimately the most important test to check whether the nonparametric dFPT is capable of doing quantification (i.e. finding the peak parameters despite being a shape estimator at the onset of the application). Our previous studies [16][17][18] (synthesized noiseless as well as noisy FIDs reminiscent of time signals encoded by in vitro proton MRS from malignant breast tissue) and [21][22][23] (FIDs encoded from the Philips polyethylene phantom [32,33]) have successfully implemented this most stringent benchmarking of the nonparametric dFPT. We will presently use the same verification procedure and the findings will be given in the Sect. 3.

Encoding time signals or FIDs
The standard encoding protocol for in vivo MRS was used to encode time signals or FIDs from a brainstem glioma in a four year-old patient. Encodings have been made at the Astrid Lindgren Children's Hospital in Stockholm, employing a 1.5T General Electric (GE) clinical scanner. Proton MRS was applied with single-voxel point-resolved spectroscopy sequence (PRESS). During encoding, water was partially suppressed by means of the conventional inversion recovery procedure [32,33]. The parameters of this acquisition are: the Larmor frequency Ethics committee approval has been obtained through a formal application (Dnr # 2007/708-31/1) to the Regional Ethics Committee at the Karolinska Institute. The Committee's evaluation and statement were that there are no ethical issues whatsoever that would preclude the implementation of this research.

Challenges in reconstruction of spectra
The acquisition parameters identical to the just outlined list have also been used in our most recent study [21][22][23] on 1 H MRS applied to the Philips polyethylene phantom 3 [32,33]. Therein, the overall success of the dFPT was shown to be striking. Therefore, it would be very important to find out how much of the unprecedented performance of the dFPT [21][22][23] could presently be replicated using the encoded in vivo MRS time signals.
The achievements from Refs. [21][22][23] for the polyethylene phantom [32,33] were in super-high resolution of the J-coupled resonances in two different chemical groups of ethanol with the proper relative intensity ratios of the lineshapes. This has been accomplished despite utilizing the short FIDs of only 512 data points and applying a weak static magnetic field, 1.5T. Moreover, the reconstructions by the dFPT from Refs. [21][22][23] were of the same spectral quality for two diametrically opposite FIDs, one encoded with and the other without water suppression.
The present work is on the FIDs encoded with water suppression in the measuring process as was indeed the case in Ref. [21]. Therefore, the findings in the work [21] are the most appropriate reference data for judging the validity of the present reconstructions by the dFPT.
Despite the identical acquisition parameters for the polyethylene phantom [21] and the present FIDs encoded by in vivo MRS from a patient, a fundamental difference nevertheless exists. It is the type of the encoded time signals that are substantially different. For the polyethylene phantom, SNR is high. Most importantly, the chemical content as well the actual concentrations of the constituent molecules from the polyethylene sphere are known [32,33]. By contrast, the present in vivo MRS time signals are of relatively low SNR. Moreover, prior to the analysis, no quantitative information is available on the metabolites of the scanned tissue (their number, kind, abundance, etc.).
In other words, for the polyethylene phantom data, the answer from an inquiry by MRS is, in principle, known and the task of the applied signal processor is to determine its ability to match the quantitative input data (the given volumes or concentrations of the chemicals in the sample). This task is still very demanding for any signal processor, as detailed comparisons between the dFPT and the derivative fast Fourier transform (dFFT) have clearly shown [21][22][23]. The net outcome of this competition promoted the dFPT as a clear winner, whereas the dFFT failed flagrantly.
On the other hand, for the in vivo MRS problem under study here, there are multiple challenges beyond those in the polyethylene phantom. The sharpest challenge is that the answer is completely unknown in every facet of the inquiry by in vivo MRS. Faced with this adversity, we resort to the intrinsic cross-validation of the reconstructions by the dFPT. Note that the dFPT mentioned throughout this sub-section is the nonparametric dFPT which provides only the qualitative estimates of the forms of the envelope lineshapes, and not the resonance signatures (peak positions, widths, heights, phases). Similarly to Refs. [21][22][23], the intrinsic benchmarking of the nonparametric dFPT will be carried out by comparisons of its envelopes with the envelopes and components reconstructed by the parametric dFPT.
In particular, if the envelopes and the components from the shape and parameter estimations by the dFPT, respectively, coincide with each other, this would constitute the benchmarking validation of the nonparametric dFPT. Such a stringent test has most successfully been passed by the dFPT for the polyethylene phantom [21][22][23]. The task now is to benchmark the performance of the nonparametric dFPT in order to see how much of the quality of the reconstructions for the polyethylene phantom can eventually be achieved for the presently examined in vivo MRS time signals.
Moreover, there is also the clinical realm of our multi-disciplinary approach to NMR spectroscopy: from physics via analytical chemistry to medicine by way of mathematics. Namely, our ultimate goal actually surpasses the just-stated inquiry as we are determined to find out the clinical significance of the output spectral data from the dFPT for the encoded input FIDs using in vivo MRS in cancer medicine. In particular, it is high time to ask whether would it be possible, using exclusively the shape estimation by the dFPT, to unequivocally identify and quantify the potential as well as the recognized cancer biomarkers (phosphocholine, choline, lactate, … ).
The clinical importance of such findings would be in empowering the physician in helping make differential diagnosis: malignant versus benign tissues of the patient examined by MRS. Such decision making would facilitate the entry of MRS into the everyday diagnostic armamentarium in hospitals, a long overdue goal. In vivo MRS made its debut several decades ago and yet, due to the lack of the clinically reliable data analyses, this is still not a standard, routine diagnostic modality. This is very unfortunate because in vivo MRS definitely has its tremendous potential particularly for early tumor diagnostics on a molecular level (no exposure to ionizing radiation, no invasive procedures), when chances for control and ultimately helping cure the disease are the most notable.
Our previous detailed study [14] on the same in vivo MRS time signals as those used presently has been focused on the nonderivative FPT (both nonparametric and parametric). Therein, both interesting and important findings have been reported on the component spectra, signal-noise separation (SNS) by pole-zero cancellation (Froissart doublets) and the like matters [12,13]. Nevertheless, the principal reason for revisiting this problem here is that we shall address it now from a very different perspective, the derivative signal processing. The motivation for such an undertaking is rooted in the astoundingly high performance of the dFPT for the two substantially different types of FIDs, one fully controlled (synthesized: noiseless, noisy) [15][16][17] and the other semi-controlled (encoded from a phantom) [21][22][23].
In theoretically-generated FIDs (simulated, synthesized), all the input peak parameters are known in advance and, thus, such time signals are evidently fully controlled. Time signals encoded from phantoms are semi-controlled in the sense that the chemical content of the examined sample is given, i.e. fixed. Here, control is not full as it depends on both data acquisition and data analysis. On the other hand, FIDs encoded by in vivo MRS are uncontrollable in the sense that no prior knowledge is available on the actual chemical content and metabolite abundance in the scanned tissue. The controlled and semi-controlled FIDs are crucial for MRS as they can both furnish robust benchmarkings of signal processing methods and establish the norms. Guided by such rules, one is in a much better position to analyze the FIDs encoded by in vivo MRS and that is what we are set to do here with the sole focus of derivative signal processing using the dFPT as a shape estimator. The dFPT as a parameter estimator will be used to intrinsically benchmark the nonparametric dFPT, as stated.

Logistics of using encoded FIDs for reconstructions
This sub-section reports on all the MRS data (input, output) envisaged for the present study at TE = 272 ms. The input data are the encoded time signals or FIDs. The output data are the reconstruction findings from computations (the latter data employ the averaged time signal, as mentioned). The fast Fourier and fast Padé transforms are used in two different forms, one is nonderivative (FFT, FPT), whereas the other is derivative (dFFT, dFPT). Fourier processing is qualitative as it provides only lineshape estimations of envelopes (no components). However, Padé processing, can be both qualitative (envelopes) and quantitative (components). To exhibit these two features, we apply the Padé nonparametric and parametric estimations in their nonderivative and derivative versions. The nonparametric and parametric Padé variants cross-validate each other as an internal benchmarking.
Fourier processing is customarily carried out with a number (usually 1-4) of zero-fillings of the input FIDs. We checked here that no substantial difference occurs in the Fourier envelopes by using 1 or 2-4 zero-fillings of the encoded FID. Therefore, we opt to zero-fill the encoded FID only once. For consistency in comparison between the Fourier and Padé processing results, to this same once-zero-padded time signal, we apply the FFT, dFFT as well as the FPT and dFPT. No apodization is employed with the encoded FID, i.e. the time signal is not multiplied by an exponential e − t ( > 0) to force a faster decay to zero. Regarding the chemical shift axis, water resonance is taken to be located at 4.61 ppm (parts per million) to match the same convention from our past study on the same problem [14].
For the real parts of complex spectral envelopes, the zero-order phase correction is utilized, i.e. the time signal is multiplied by an exponential e i 0 . Of course, this phase factor is immaterial for any magnitude mode spectrum. The main emphasis in derivative signal processing is in the phase-insensitive magnitude mode of spectra. The more important reason for preferring the magnitude mode of derivative envelopes is in effectively suppressing the side lobes. These side lobes accompany the given main peak and complicate the interpretation and attempts at quantification when the real part of a complex spectrum is used. For all the present illustrations, the Padé spectra are computed with a fixed model order, K = 258, which exhausts about 1/2 of the encoded FID whose full length is N = 512. With these remarks, we can now pass to the presentation of the illustrations by way of the most salient figures (5 in total). Each figure has several panels ranging from 5 to 10. Figure 1 is the most all-inclusive as it combines the input and output data. It plots the time signals and spectral envelopes whose real parts and magnitudes are displayed as the nonderivative, derivative, phase-uncorrected and phase-corrected lineshape profiles. Panels (a-d) are on the time signals (for simplicity no zero-filled parts are graphed here), whereas panels (e-h) are on spectral envelopes from the FPT (nonderivative) and the first derivative within the dFPT. The FIDs without and with phase correction are on the top panels: (a, b) for the former, and (c, d) for the latter time signal. The real parts of the FIDs are on panels (a, c), whereas the corresponding imaginary parts are on panels (b, d).

The most salient illustrative figures
The zero-order phase 0 is usually chosen manually. Instead, we obtain It is observed in Fig. 1e that the lineshape of the envelope Re(FPT), reconstructed with the unphased FID ( 0 = 0), is highly asymmetrical, deformed and far from positive-definiteness. As such, this spectral profile is not amenable to any attempt aiming at quantification. On the other hand, the lineshape of the envelope Re(FPT) from Fig. 1f, retrieved with the phased FID ( 0 = −1.4657 rad) appears to be more symmetrical with respect to the chemical shift axis, less deformed and closer to positivedefiniteness, at least in some parts of the spectrum.
Nevertheless, despite such a qualitative 'improvement', the chances are still very slim (if any) for a meaningful quantification e.g. by integration or fitting through Lorentzians, Gaussians or their combinations (Voigtians). Integration would not be practical for the great majority of the peaks because neither the upper nor the lower integration limits (bounds) can confidently be determined. Fitting by Lorentzians, Gaussians or Voigtians would not be a better alternative either, since it would invariably yield some false metabolites (over-modeling) and/or fail to recover some true metabolites (under-modeling), both amounting to the ubiquitous nonuniqueness.
Precisely the same type of observations would also be made if instead of the FPT, we used the FFT or any other processors to reconstruct the envelope for the unphased (panels a,b) and phased (panels c,d) time signals in Fig. 1. This is the 1 3 common feature shared by all the nonderivative processors for computations of various modes of the given total shape spectrum. It is for these reasons (particularly for the predominantly fitted FFT total shape spectra) that in clinical applications of in vivo MRS to the human brain, one invariably finds the statements of the type: merely a few (3)(4)(5)  The reason for focusing almost exclusively on the real parts of complex envelopes in clinical MRS is in the hope that by phasing the encoded FIDs an approximately absorptive mode of the spectrum could be produced. An absorptive Lorentzian lineshape can most readily be interpreted in terms of the resonance location, FWHM and the peak area, which is proportional to the number of the resonating protons (and, hence to the concentration of the associated metabolite -the most relevant information for diagnostic purposes). However, as encountered throughout the MRS literature as well as in our Fig. 1f, these hopes appear to be on shaky ground, as therein the real part mode is anything but absorptive for the great majority of the peaks. This occurs because a single phase 0 is unable to mitigate the unavoidable admixture of the dispersive lineshapes in the real parts of the given complex spectrum. Not much improvement would follow when employing e.g. an additional first-order phase 1 for yet another correction of the encoded FIDs.
Given this unenviable prospect of MRS in medical diagnostics with clinical scanners (1.5 and/or 3T), one could provisionally hope that at least one trouble (the lack of positive-definiteness) could be avoided altogether by resorting to the magnitude mode envelope which is given in Fig. 1g. Alas, this attempt has not met with the success either because the magnitude envelope |FPT| (g) increases the linewidths by √ 3 relative to Re(FPT). This is evident from panels (g) and (f). Any linewidth broadening in |FPT| leads to more overlaps of adjacent resonances and this amounts to resolution deterioration and lowering of SNR.
However, the goal of MRS is precisely the opposite: increased resolution and SNR. From all the evidence available thus far in the MRS literature such a goal is definitely unattainable by nonderivative conventional shape estimations. Therefore, some other alternatives should be tried by preserving the positive-definiteness of the magnitude mode envelope |FPT| (g), but hopefully without linewidth broadening. In other words, at least in an assumed first-order improvement of the envelope, the aim is to merge together the best parts of the two 'worlds', |Re(FPT)| (f) and |FPT| (g). If that would be feasible, we would then have a solid platform onto which to build some systematically improved higher-order spectra associated with the encoded FIDs.
To follow this initiative towards the simultaneous improvement of resolution and SNR, good hints can be found in perusing the ideal Lorentzian complex spectrum with a single resonance, L( ) ≡ d r ∕( − r + i r ) corresponding to a synthesized FID whose amplitude d r ≡ |d r |e i r is a constant real number ( r = 0 ). With no phase in d r , the real part Re(L) of L( ) is a perfect bell-shaped absorption with no admixture of the dispersive part Im(L). The key question is then: which additional Journal of Mathematical Chemistry (2021) 59:2133-2178 transform could be applied to L( ) so that the magnitude of the resulting spectrum, denoted by e.g. |L � ( )|, would have the same width as the pure absorption Re(L)?
The answer is that L( ) should be subjected to a special supplementary transform given by the first derivative, d∕d , in which case the unknown spectrum L � ( ) becomes L � ( ) = (d∕d )L( ).
As such, the sought best parts of Re(L) (the most readily interpretable peak width via the FWHM, the peak area, etc.) and positive-definiteness (phase-insensitiveness) of |L( )| are joined together in |D 1 L( )| ≡ |(d∕d )L( )|. To return to the corresponding FPT, it suffices to note that L( ) itself is, in fact the simplest paradiagonal FPT of the type P K−1 ( )∕Q K ( ) for K = 1 as L( ) = P 0 ( )∕Q 1 ( ), where P 0 = d r and Q 1 ( ) = − r + i r .
Therefore, this prelude is a motivation to try the chance with the magnitude mode envelope |D 1 FPT| of the complex first derivative total shape spectrum D 1 FPT ≡ (d∕d )FPT. Recall that, according to the nomenclature in Eqs. (2.30)-(2.32), the acronym FPT stands for the diagonal fast Padé transform, FPT ≡ P K (z −1 )∕Q K (z −1 ) where, as before, z −1 is the harmonic variable containing the sweep frequency f = Re( ) via z −1 = e −2 i . The first derivative magnitude spectrum |D 1 FPT| is displayed in Fig. 1h. Therein, it is seen that |D 1 FPT| (h) is superior to all the three formerly considered envelopes, Re(FPT) (e, f) and |FPT| (g). In particular, |D 1 FPT| (h) simultaneously improves resolution and SNR relative to Re(FPT) (e, f) and |FPT| (g).
As desired, the linewidths of resonances from Re(FPT) (f) and |D 1 FPT| (h) are the same. Yet, as opposed to |Re(FPT)| (f), it is evident that |D 1 FPT| (h) can achieve a notable splitting of some of the diagnostically most important metabolites, e.g. choline, Cho, and phosphocholine, PC, around 3.2 and 3.22 ppm, respectively. At first glance, it may seem contradictory that between the two envelopes with theoretically the same linewidths of all the resonances, one spectrum |D 1 FPT| (h) is better resolved than the other Re(FPT) (f). This puzzle is solved on account of having a better SNR in |D 1 FPT| (h) than in Re(FPT) (f). The higher level of noise in Re(FPT) (f) increases the overlap of adjacent peaks and this explains the appearance of e.g. only a 'total choline' (tCho) peak as a compound resonance around 3.2 ppm in Re(FPT) (f).
Splitting of the 'total choline' peak (Cho+PC) in |D 1 FPT| (h) is seen not to descend down to the chemical shift axis. Therefore, such a resolving of the formerly overlapped peaks should be viewed only as an indication of a potential trend toward a deeper separation of the Cho and PC resonances. This is precisely a sort of hope which was fostered prior to arriving at Fig. 1h aiming at the proper directionality: try to systematically improve |D 1 FPT| (h) by some higher-order spectral envelopes. Here, the clause a 'higher-order' envelope now can be identified as a 'higher-order derivative' envelope. Nice as it may seem though, |D 1 FPT| is not yet quantifiable because of the still present background baseline which prevents determination of the integration limits for unambiguous finding of the peak areas by a numerical quadrature (e.g. Gauss-Legendre rule or the likes).
Overall, Fig. 1 is a 'soft' introduction to derivative signal processing illustrated in the setting of the Padé versatile methodologies. It is a pedagogically-minded step-wise introduction which singled out only the first derivative Padé spectrum, |D 1 FPT| (Fig. 1h). Pedagogical, because it exposes the drawbacks of the customary, nonderivative estimations and provides the pathway for potential surmounting these obstacles by using the first derivative envelopes.
Step-wise, because it highlights the advantages of |D 1 FPT| (Fig. 1h) relative to Re(FPT) (f) and hints on possible improvements by passing from the first-order to some high-order derivatives. It is deemed that this 'preface' to the full-blown display of all the derivatives of interest |D m FPT| (m = 1, 2, 3, …) is instructive because it brings in a direct contact the first derivative envelopes with the more familiar nonderivative envelopes.
In such a way, the reader acquires the important initial glimpse of the multiple advantages already on the level of the first derivative envelopes. Thus, without being overwhelmed by all the derivative envelopes until the converged spectra have eventually been attained, the reader will hopefully be motivated to proceed further to see what else is lying ahead in the store of derivative estimations. This brings us to Fig. 2. Figure 2 also has an important story to tell in a fashion 'all that glitters is not gold' when it comes to an anticipation that any derivative estimator would outperform its own nonderivative counterpart. That this expectation is false is evidenced in Fig. 2 (ten panels, a-j), which compares two derivative estimators, one based on Padé and the other on Fourier signal processings. The left (a-e) and right (f-j) columns of Fig. 2 are for Fourier and Padé estimations, respectively. The top two panels (a, f) of this figure are for nonderivative processing. Panels (b-e) and (g-j) are for derivative estimations using the derivative operator D m = (d∕d ) m (m = 1-4).
The applied one zero-filling of the FID artificially extends the original full signal length N = 512 to 1024. The result is a trigonometric interpolation in the corresponding Fourier spectral envelope [35][36][37][38]. Its manifestation in |FFT| (Fig. 2a) is in the appearance of many tightly packed sharp-edged wiggles that often shoulder the nearby stronger peaks. Such an effect of zero-filling of the FID is largely smoothed out in the Fourier derivative envelopes, as clear already in |D 1 FFT| (Fig. 2b). This is rather interesting because it is just the opposite to what one would expect. Namely, the derivative operator is theoretically anticipated to enhance more the thin sharpedged peaks than the broader resonances. The derivative operator partially or completely suppresses these wiggles in the Fourier spectra because it treats them as noise-like corruption, i.e. it recognizes their faulty nature (not being a part of the physical, genuine content of the originally encoded FIDs).
In the literature of NMR spectroscopy, it is too often stated that zero filling of FIDs improves resolution in the corresponding Fourier spectra. If this were indeed true, it would mean that zero filling could be credited to the emergence of some additional information. However, this is impossible since the entire information is already contained in the encoded FID of the original total signal length N [38]. The alleged 'improved' spectral resolution with a zero-filled FID is manifested in the wiggles that, however, disappear with the application of the derivative operator.
The overall impression from looking at panels (a-e) in Fig. 2 is that Fourier derivative signal processing has a grand fall. It is at odds with all the facets of the hopedfor expectation from the well-conceived derivative estimation. This occurs because in the dFFT, the linewidths are broadened and SNR deteriorated already with the  1-4). The situation gets much worse for m ≥ 5 (not shown here to avoid clutter, however, in this context see Fig. 8 in Ref. [23]). Even the first derivative |D 1 FFT| (Fig. 2b) broadens some of the resonances that happen to be clinically important, as clear from e.g. the frequency band 1.3-1.4 ppm, which is the interval of the diagnostically very relevant metabolite: lactate, Lac. The other two metabolites of critical diagnostic significance (Cho, PC), lying within a narrow frequency range (3.20-3.22 ppm), are observed to have their distinct 'ups and downs'. Their individual resonances are invisibly glued together in a compound and wide 'total choline' peak (tCho = Cho + PC) in the nonderivative (m = 0) envelope |FFT| (Fig. 2a).
However, they become partially visible (with a clear albeit narrow splitting) in the first derivative (m = 1) envelope |D 1 FFT| (Fig. 2b). Of importance to note (for the subsequent discussion) is that upfield, in the immediate neighborhood of the PC resonance, there is a quite pronounced peak of taurine (Tau). Other weak and strong, narrow and wide resonances superimposed on a lifted background baseline follow in the upfield direction of increased sweep frequencies.
Even a better separation of Cho and PC is seen in the second derivative (m = 2) envelope |D 2 FFT| (Fig. 2c). However this good progress is upset by the adjacent wide and taller taurine peak which begins to overwhelm the PC and Cho pair. This obstacle in Fig. 2c is further exacerbated by the systematic linewidth broadenings especially of the intense resonances in the upfield direction. Moreover, linewidth broadenings also occur downfield in |D 2 FFT| (Fig. 2c). For example, the lactate doublet becomes wider in |D 2 FFT| (Fig. 2c) than in |D 1 FFT| (Fig. 2b). Such trends throughout the envelope are accompanied with the lifted background baseline indicating that |D 2 FFT| begins to decrease SNR. This is not a transient weakness of the dFFT because, eventually, the full breakdown of the Fourier-based derivative estimations becomes evident from the third derivative (m = 3) envelope |D 3 FFT| (Fig. 2d). Herein, the mentioned increase of the taurine peak has almost swamped its neighbors, the Cho and PC resonances. Likewise, in the downfield region (1.0-1.5 ppm), the elevated background baseline increased the overlap of the lactate doublet by diminishing the dip in between these two peaks. All the other resonances are also widened in |D 3 FFT| (Fig. 2d) resulting in a further worsening of SNR relative to |D 1,2 FFT| (b,c).
This dire situation becomes a flagrant failure in the fourth derivative (m = 4) envelope |D 4 FFT| (Fig. 2e). Therein, the wide and taller taurine resonance practically annihilates the Cho and PC peaks. The background baseline is supposed to be flattened with the increased derivative order m. However, just the opposite happens when passing from e.g. m = 2 in |D 2 FFT| (Fig. 2c) to m = 4 in |D 4 FFT| (Fig. 2e). By going from |D 3 FFT| (Fig. 2d) to |D 4 FFT| (Fig. 2e) basically all the physical information is lost, including that of diagnostic relevance (Cho, PC, Lac). In fact, the three derivative envelopes |D 2,3,4 FFT| (c, d, e) are worse than their nonderivative seed spectrum, |FFT| (a).
The reason for the constantly increasing background baseline with the augmented derivative order m in the dFFT is the higher and higher tail of the water resonance located at 4.61 ppm (outside the frequency interval in Fig. 2). The inability of the dFFT to narrow the bottom of the water resonance (i.e. to localize it tightly around the resonance frequency 4.61 ppm) results in the water tail extension downfield. In fact, the right wing of the water peak is partially seen around 4.5 ppm on all the panels on the left column of Fig. 2. This right wing of the water peak is the tallest spectral structure in |FFT| (a). Temporarily, it ceases to dominate in |D 1,2,3 FFT| (b, c, d) only to rise again in |D 4 FFT| (e). In fact, much of the disagreeable features of the dFFT is due to the residual water peak whose widening with the increased derivative order m is detrimental to the quality of the remainder of the spectrum, as is obvious in |D 4 FFT| (Fig. 2e).
Overall, the failure of the dFFT for derivative estimation can be traced back to the way in which the Fourier derivative envelops are reconstructed. They stem from the application of the derivative operator D m = (d∕d ) m (m = 1, 2, …) to the standard nonderivative FFT complex envelope. This application produces the power function, t m = (n ) m , which multiplies the time signal c n . The multiplication term is increased at larger signal numbers n at which noise prevails in all encoded FIDs. Therefore, instead of being suppressed by the derivative operator D m , noise is generally enhanced in spectral envelopes from the dFFT. This leads to more overlaps due to more line broadenings as noise begins to fill in the dips in adjacent peaks. Ultimately, all the physical information is irretrievably lost in the dFFT, as shown most clearly in panels (c-e) from Fig. 2.
The right column (f-j) of Fig. 2 for Padé processing paints a totally different picture of derivative estimations. Spectra |FPT| (Fig. 2f) and |D 1 FPT| (Fig. 2g) have already been fully analyzed with Fig. 1. They are reproduced here for completeness as well as for comparison with |FFT| (Fig. 2a) and |D 1 FFT| (Fig. 2b). The Padé-based nonderivative envelope |FPT| (Fig. 2f) is quite of a similar spectral structure as its Fourier neighbor |FFT| (Fig. 2a). The main difference is the lack of any visible evidence that zero-filling of the FID produces the spectral wiggles, as opposed to |FFT| (Fig. 2a). Should such wiggles be present also in FPT = P K ∕Q K they would be suppressed. The mechanism for such a suppression is in recognizing that these wiggles are a kind of error (in the sense of being some unphysical, spurious part of the content).
It is the rational form P K ∕Q K of the FPT envelope that is a self-correcting response function capable of canceling out any error occurring in P K and Q K . When a zero-filled FID is subjected to the FPT, both sets of the coefficients {p r } and {q s } of P K and Q K , respectively, become corrupted with noise-like wiggles that, in turn, are largely canceled in the quotient P K ∕Q K . Such an occurrence is reminiscent of the ubiquitous finding in many experiments that errors incurred in measuring two observables A and B are significantly reduced in their ratio A/B. This is an empirical evidence, of course.
However, a similar cancellation in the spectrum P K ∕Q K from the FPT is not empirical. Rather, it is based on the theoretically established mechanism called pole-zero cancellation. The spectral wiggles might appear transiently in P K ∕Q K but, due to their Padé-identified spuriousness, they are treated as a noise-like function and, as such, categorized as Froissart doublets. The wiggle-type spuriousness becomes equi-distributed between the numerator P K and denominator Q K polynomials as zeros and poles respectively of the FPT envelope P K ∕Q K . Then pole-zero cancellation takes place because the spectral wiggles (due to the zero-filled FID), being associated with the common roots of P K and Q K , are canceled out in P K ∕Q K . This is the denoising Froissart filter at work [12,13].
To pass to the Padé-based derivative estimations by way of the dFPT, we can now inspect the panels (g-j) of Fig. 2. The Padé first derivative |D 1 FPT| (g) is seen to outperform the corresponding Fourier counterpart |D 1 FFT| (b) regarding both resolution and SNR. As opposed to line broadening in |D 1 FFT| (b), line narrowing is seen in |D 1 FPT| (g) as best evidenced in the lactate region. Importantly, the water peak narrowing is noted to result in its much smaller remnant (more than 3 times smaller) near 4.5 ppm in |D 1 FPT| (g) than in |D 1 FFT| (b). The underlying shrinking of the water peak tail results in flattening of the background baseline in |D 1 FPT| (g).
These favorable features and the overall trend in |D 1 FPT| (g) are further hugely progressed in the second derivative |D 2 FPT| (h). Herein, the great majority of resonances are sharply resolved very closely to the chemical shift axis. This is possible because of the continued flattening of the background baseline since in |D 2 FPT| (h), the right wing of the water peak residual near 4.5 ppm is about 15 times smaller than that in |D 2 FFT| (c). In other words, SNR is hugely improved in |D 2 FPT| (h) with respect to |D 2 FFT| (c).
Several tightly overlapped resonances can still be seen in |D 2 FPT| (h). More precisely, these are five pairs of the overlapped resonances out of 54 resonances in |D 2 FPT| (h). Among these resonances, the most pronounced are the two pairs of the overlapped peaks in the region 3.95-4.10 ppm. The five pairs of overlapped peaks are highlighted by the red star symbols on the envelope (h) for the derivative order m = 2. The purpose is to facilitate visualization of the splittings in each of these resonance pairs when the derivative order is augmented to m = 3 (i) and m = 4 (j).
Note that the separation of PC from Cho near 3.20-3.22 ppm is strongly improved in |D 2 FPT| (h) relative to |D 1 FPT| (g). Moreover, in |D 2 FPT| (h), the taurine metabolite near 3.3 ppm is assigned to a narrow, sharp and tall resonance. This is opposed to |D 2 FFT| (c), where a broad taurine peak begins to mask the PC resonance, as discussed in the analysis of the left column of Fig. 2.
By going to the third derivative |D 3 FPT| (i), three out of five overlapped pairs of peaks become fully resolved down to the chemical shift axis. The implication is that here some 50 resonances reached their converged status with their bottoms (bases) lying on the chemical shift axis and being completely isolated from their immediate neighbors. This means that the integration limits are unambiguously definable to facilitate determination of the peak areas by any of the usual quadrature rules (e.g. Gauss-Legendre, etc.). Notably, however, the 4th and 5th overlapped pairs of resonances around 4 ppm still persist to evade full splitting in |D 3 FPT| (i). As noted earlier, the PC and Cho peaks are diminished in |D 3 FFT| (d) by the dominating broad taurine peak around 3.3 ppm. By contrast, the PC and Cho peaks retain their full separation and visibility in |D 3 FPT| (i), while simultaneously the taurine resonance, close to 3.3 ppm, keeps on getting thinner and sharper.
Since the operator |D 3 | was found to be insufficiently powerful to resolve the two remaining pairs of overlapped peaks near 4 ppm in Fig. 2i, we resorted to the fourth derivative whose action is seen in |D 4 FPT| (j). Both these pairs of the overlapped resonances are now fully resolved all the way down to the chemical shift axis. Moreover, all the other 50 peaks are further refined in |D 4 FPT| (j).
To better exhibit the appearance of all the five overlapped pairs of peaks and especially the 4th and 5th pairs (close to 4 ppm), the dynamic range of the spectral intensities on the ordinate is reduced. This cuts the tops of the tallest seven resonances in the envelope spectrum. However, there is no issue with such a convenient shrinking of the ordinate, since the truncated seven resonances are fully visible through their converged lineshapes in |D 3 FPT| (i). Crucially, the recognized cancer biomarkers (PC, Cho, Lac, … ) are in full display in |D 4 FPT| (j), continuing the pace which is steadily progressing from panels (g) to (j) on the right column of Fig. 2.
Overall, the dFPT from panels (g-j) met with success in all the facets of the expected improvements by derivative signal processing, particularly regarding enhanced resolution and SNR. This is the result of a systematic synergism: linewidth narrowing and the background baseline flattenings. Such an achievement is remarkable on its own right, but also by reference to our recent studies on the polyethylene phantom [21][22][23]. While in Refs. [21][22][23] there were altogether 9 peaks to resolve (seven J-coupled peaks of ethanol, a triplet and a quartet, accompanied by two singlets, methanol and acetate), in the present study some 54 resonance are completely split apart. In other words, the challenge in the current work is far more demanding and yet the performance of the dFPT is seen to be in a fully successful display.
These are the mathematical aspects of the unprecedented power of derivative signal processing within the nonparametric dFPT, especially relative to both dFFT and FPT (nonderivative, nonparametric, Fig. 2f). However, the validation of the dFPT is not exhausted by the sole mathematical aspects. Quite the contrary, it encompasses the clinical aspects of the problem under study. Notwithstanding the remarkable performance of the dFPT for the polyethylene phantom [21][22][23], as a specially important test for an innate sample, the present study is on in vivo MRS time signals encoded from a patient with a brain tumor (glioma).
Of crucial diagnostic significance is to achieve a very clear delineation of the recognized cancer biomarkers (e.g. Lac near 1.3-1.4 ppm as well as Cho and PC around 3.20-3.22 ppm), according to |D 3 FPT| (i) and |D 4 FPT| (j). Never before have Cho and PC been split apart by shape estimations for time signals encoded at clinical scanners (1.5, 3T). This is all the more remarkable given that fact that we used a very short FID (only 512 data points) encoded at a clinical scanner of low magnetic field (1.5T).
The just effected comparison in Fig. 2 between the Padé and Fourier estimations using their nonderivative and derivative versions is important and interesting. Of course, Fourier processing is limited to estimations of envelope lineshapes alone. Therefore, to make a fair comparison in Fig. 2 between the respective pairs {FFT, dFFT} and {FPT, dFPT} , we also used shape estimations by means of the nonparametric versions of the nonderivative and derivative fast Padé transform.
Instructive as it is, comparisons between the two estimator tandems {FFT, dFFT} and {FPT, dFPT} are not the end of the story, as far as thorough testings and benchmarkings of the dFPT are concerned. A deeper level of the validity checking of the nonparametric dFPT is further required. It can be provided by comparisons between the nonparametric and parametric versions of the dFPT. The parametric dFPT gives both the components and envelopes. Therefore, the first comparison, while staying in the home of the Padé methodology, is to juxtapose the two sets of envelopes, one from the parametric and the other from the nonparametric estimations (both nonderivative and derivative), as done in Fig. 3. Then in the second go comes the most stringent of all testings, the detailed comparisons between the components (parametric dFPT) and envelopes (nonparametric dFPT), as given in Fig. 4. By definition, all the component spectra build the envelope spectrum while working with the parametric dFPT. However, the key concern is quite a fundamentally different issue, which is a potential relationship between the envelopes from the nonparametric dFPT and the parametrically reconstructed dFPT components. Now one may wonder why are these two additional benchmarkings (Figs. 3, 4) needed given that all the resonances in Fig. 2j for the nonparametrically retrieved spectrum |D 4 FPT| are completely isolated and fully resolved down to the chemical shift axis? The need for these supplementary checkings is justified to independently confirm the veracity of the total number of resonances as well as their peak areas found by the dFPT as a shape estimator.
Even a slight perusing of Fig. 3 immediately reveals full concordance between the envelopes on the left (a-e) and right (f-j) columns due to the parametric and nonparametric Padé estimations, respectively: nonderivative (a, f) and derivative (b-e, g-j). A more detailed inspection of the envelopes |D 0−4 FPT| (f-j) for the nonparametric Padé processing shows that they totally coincide with their parametric Padé counterparts |D 0−4 FPT| (a-e), respectively, where |D 0 FPT| = |FPT| (D 0 = 1). This holds true for every single spectral resonance, irrespective of its intensity (small, middle, large). The same conclusion also applies to the highlighted five pairs of the overlapped resonances, marked by the blue (c-e) and red (h-j) star symbols on the left and right columns of Fig. 3, respectively. We can then safely conclude that, within both the nonderivative and derivative envelopes, the highest accuracy of the nonparametric pair of processors {FPT,dFPT} is amply confirmed by the corresponding parametric estimators {FPT,dFPT}, respectively.
The maximum possible extent of agreement between the nonparametric and parametric Padé estimations of envelopes, including their respective nonderivative and derivative branches, is gratifying. Nevertheless, it is still a necessary, but not sufficient condition for the full, all-encompassing cross-validation of the nonparametric dFPT. The reason is that agreement seen in Fig. 3 does not guarantee that some of the weak and possibly diagnostically important resonances might be absent from the nonparametric or parametric or both envelopes. Such weak resonances may not be able to perceivably change the lineshape of the envelopes from the left and right column of Fig. 3. This latter observation is precisely the kind of rationale which is a motivation to extend the comparisons by encompassing juxtapositions of the nonparametric envelopes and the components (parametric FPT, dFPT). With such a guidance we are now directed toward Fig. 4. Therein, the parametric and nonparametric Padé reconstructions are on the left (a-e) and right (f-j) columns, respectively. The nonderivative spectra are on the first row (a, f). On the other hand, the derivative spectra are configured on the second-to-the-fifth rows: |D 1−4 FPT| Comp (b-e) versus |D 1−4 FPT| Tot (g-j), respectively. Here, following the convention (2.31), the subscripts 'Comp' and Of course, the nonderivative spectra |FPT| Comp (a) and |FPT| Tot (f) are not directly comparable as there are many component curves on the left and only one curve for the envelope on right hand side. Yet, it is interesting to juxtapose |FPT| Comp (a) and |FPT| Tot (f) to make in evidence why and how the nonderivative envelope hides its constituents. This hiding in |FPT| Tot (f) is a consequence of either destructive and constructive (or both) interference effects of complexvalued quasi-absorptive and quasi-dispersive component lineshape profiles. Naturally, the water resonance peak also takes part in such interactions.
Note that the intense right wing of the water resonance near 4.5 ppm, as the strongest spectral structure in |FPT| Tot (f), cannot possibly result from the components in |FPT| Comp (a) in the same region (i.e close to 4.5 ppm). This is explained by the mentioned occurrence that, to avoid clutter, the lineshapes of the resonances centered at frequencies outside the window 0.35-4.535 ppm, are not graphed in |FPT| Comp (a). The water resonance is one of such remote peaks. Had we plotted e.g. the right wing of the water peak in |FPT| Comp (a), it would be almost as high as the tallest peak at 3.95 ppm.
We now proceed to analyze the derivative spectra by comparing the components |D 1−4 FPT| Comp (b-e) with the envelopes |D 1−4 FPT| Tot (g-j). The most dramatic change occurs in the component spectra when passing from the nonderivative |FPT| Comp (a) to the first derivative |D 1 FPT| Comp (b). In |D 1 FPT| Comp (b), the first gigantic suppression occurs with the three broadest resonances centered near 1.6, 3.8 and 4.25 ppm. In particular, there is no trace left whatsoever from the formerly widest resonance centered near 1.6 ppm.
The other two wide resonances centered close to 3.8 and 4.25 ppm are practically annihilated and only some minuscule, low-lying remainders can barely be noticed. The other formerly less wide peaks centered close to 2.6, 2.8 and 3.05 ppm still have their visible remainders. However, the relatively strong and wide peak near 4.45 ppm from |FPT| Comp (a) has its clear residual in |D 1 FPT| Comp (b).
All told, already the first derivative operator D 1 is capable of hugely diminishing the widest resonances. They might be pushed into or very near the chemical shift axis, but they are not gone for good. They can reappear, by reducing the dynamic range on the ordinate. This would result in cutting the top parts of the tallest peaks and in the emergence of the low-lying resonances. In a similar context, this device has earlier been encountered in Figs. 2j and 3j.
On the right column, the nonparametric envelope |D 1 FPT| Tot (g) begins, at least at some frequencies, to approximately mimic its component counterpart spectra |D 1 FPT| Comp (b). This represents a good step forward relative to the notable dissimilarity between the corresponding nonderivative spectral pair, |FPT| Comp (a) and |FPT| Tot (f). Nevertheless, such a step is far from being definitive, as many questions remain unanswered by |D 1 FPT| Tot (g).
That is why we need to pass to comparisons between the second derivatives, |D 2 FPT| Comp (c) and |D 2 FPT| Tot (h). It is seen here that the formerly wider peak near 4.45 ppm is now a minuscule resonance. This clears the way for the firm emergence of all the narrow and sharp resonances that survived the application of the derivative D 2 operator. The situation on the right column with the nonparametric total shape spectrum |D 2 FPT| Tot (h) has greatly improved with respect to |D 1 FPT| Tot (g).
The improvement is manifested in the fact that in |D 2 FPT| Tot (h) most (if not all) resonances from this envelope are quite reminiscent of the components from |D 2 FPT| Comp (c). Of course, this is only a qualitative reminiscence because of the presence of a number of overlapped resonances in |D 2 FPT| Tot (h). Following the practice from Figs. 2 and 3, the five pairs of tightly overlapped resonances are emphasized. They are marked by the blue and red star symbols on the components (c) and envelope (h), respectively, for the derivative order m = 2. This makes it easier to follow the separations in each of these resonance pairs when the derivative order is increased to m = 3 (d, i) and m = 4 (e, j).
It is the presence of the overlapped peaks in |D 2 FPT| Tot (h) that demands comparison of the third derivative spectra |D 3 FPT| Comp (d) and |D 3 FPT| Tot (i). Herein, it is pleasing to see a significant refinement of the sought correspondence between |D 3 FPT| Comp (d) and |D 3 FPT| Tot (i). Many resonances in the envelope |D 3 FPT| Tot (i) attained their one-to-one correspondence (meaning coincidences of the peak positions, widths, heights and, hence, peak areas) with the components |D 3 FPT| Comp (d). Yet, the task is not completed with the third derivative due to the still present overlapped resonances (four in total) in |D 3 FPT| Tot (i) within the narrow range 3.8-4.05 ppm, as discussed with Fig. 2.
Thus, we are led to compare the fourth derivatives |D 4 FPT| Comp (e) and |D 4 FPT| Tot (j). As before in Figs. 2 and 3, we need not consider the full dynamic range for the lineshape intensities on the ordinates, because especially the tallest seven resonances have fully been isolated in the third derivative envelope |D 3 FPT| Tot (i). As such, the reduced dynamic range suffices for |D 4 FPT| Comp (e) and |D 4 FPT| Tot (j) to exhibit most transparently the weak resonances. It is seen on the fifth row (e,j) that all the overlapped resonances are resolved in the nonparametric envelope |D 4 FPT| Tot (j). This has already been noted in Figs. 2 and 3. However, now there is a new input into the analysis, the gold standard, |D 4 FPT| Comp (e), to scrutinize |D 4 FPT| Tot (j). And that gold standard tells us that there is a complete coincidence of every single resonance from the nonparametrically reconstructed envelope |D 4 FPT| Tot (j) and the associated components |D 4 FPT| Comp (e) from the parametric dFPT.
We can then state with full certainty that what we have seen earlier in Figs. 2j or 3j for |D 4 FPT| Tot is indeed the true answer to the exact quantification by means of the fourth-order derivative exclusively shape estimations in the dFPT. However, this could not be stated before the testings in Figs. 3

and 4 have been completed.
Of paramount importance were the comparisons between the total shape spectra (nonparametric dFPT) and components (parametric dFPT). Complete coincidence of the nonparametric envelopes |D 4 FPT| Tot (j) and the parametrically recovered components |D 4 FPT| Comp (e) represents itself the exact solution of the quantification problem. In other words, the nonparametric fourth-order derivative dFPT solves the quantification problem (i.e. finds the peak positions, widths, heights) of all the physical, genuine resonances without ever explicitly rooting the characteristic polynomial Q K , as opposed to the parametric FPT and dFPT. Figure 5 gives a summary of the analysis as it encapsulates all the key facets of derivative signal processing under the umbrella of the nonparametric dFPT. Herein, there is no need to replot the components from the parametric dFPT as this has already been analyzed in Fig. 4. In Fig. 5, four derivative orders are used again of the operator D m (m = 1-4). This figure is also headed with the starter (m = 0), which is the seed or reference spectrum, the nonderivative envelope, |FPT| (a). Then, the dFPT begins its run as a purely shape estimator, just like its parent predecessor FPT (a). Eventually, the nonparametric dFPT completes the main task of MRS as a parameter estimator by going step-wise from its first (b) to the fourth (e) derivatives.
Integrations and derivatives are two opposite operations. They may annul each other's effect. An elementary example is the indefinite integral of the time power function t m with the result t m+1 ∕(m + 1) plus a constant. On the other hand, the first derivative d∕dt of the latter integration outcome recovers exactly the starting function, (d∕dt)t m+1 ∕(m + 1) = t m . Generally, integrations may smooth out (in an averaging sense) the delicate differences among the adjacent parts of the given functional curve. By contrast, derivatives may sense the slightest curvature changes and, thus, could accentuate the differences among the neighboring segments of the examined curve. Hence, there is a compensatory effect if derivatives are applied after integrations.
In signal processing, what happens to be hidden as a result of smoothing by e.g. the Fourier integration over a time interval 0 ≤ t ≤ T of an FID weighted with the exponential e −2 i t , has a chance to be recovered by a derivative operator D m = (d∕d ) m for some values of the order m. However, this does not happen in the dFFT. The reason is that the dFFT processes the altered time function (−2 it) m c(t), where the D m operator produces the power-type apodization (−2 it) m . Such a particular apodization may annul the expectations from derivative estimations because it weighs heavily the later encoded FID data points at which noise prevails. No similar failure is present in the dFPT, which processes the original time signal c(t) with no apodization whatsoever. This fact enables the D m operator to bring its two main advantages within the dFPT, improved resolution of resonances and reduced SNR.
The dFPT begins its journey as a purely shape estimator. Its seed or reference envelope, the nonderivative FPT is shown on panel (a) in Fig. 5. Herein, as in Figs. 1-4, a metabolite map gives a number of resonances assigned to the known tissue molecules within 0.35-4.535 ppm. Also marked are some macromolecules (e.g. For reliably aiding clinical in vivo MRS throughout diagnostics, including cancer medicine, it is vital to have an adequate procedure for processing all macromolecules. This is especially important for those macromolecules (lipids, proteins, … ) that invariably block detection of cancer biomarkers (Lac, Cho, PC). The dFPT possesses this sought procedure which is called sequential quantification. It isolates and quantifies the sharp resonances first because their growth pace with the increased derivative order m is fast. Macromolecules as heavy molecules are associated with broader and shorter resonances that are, because of their slower growth rate, forced by the derivative operator D m to merge into the background baseline.
Once the well-isolated sharp resonances have been quantified, the formerly broader macromolecular resonances would emerge from their hideouts in the background baseline. In such a way, macromolecular peaks can also be visualized, isolated and quantified. This visualization can be made merely by some suitable plotting devices. For example, there is a possibility to use two different scales on the ordinates, one for taller and the other for shorter resonances so as to have both kinds of peaks visualized simultaneously. Alternatively, a simple reduction of the dynamic range on the ordinate for the spectral intensities can be made. This would cut-off the top parts of tall sharp resonances, thus making in evidence the presence of the formerly obscured macromolecular peaks.
The sequential quantification is not limited to macromolecular resonances alone. Quite the contrary, it can be done for any shorter resonances that the derivative operator D m makes insufficiently visible after isolation of thinner peaks. Such is the situation seen clearly on panel (j) in Figs. 2-4 for |D 4 FPT|. It is also magnified in Fig. 5e. The smaller peaks seen therein could further been narrowed and sharpened using |D m FPT| for m ≥ 5 (not shown here as Fig. 5 is busy enough).
This sequential quantification constitutes a novel kind of quantification of envelopes, a visualized quantification. To extract the peak areas e.g. from |D 3,4 FPT| (d,e) in Fig. 5 to deduce the metabolite concentrations becomes an indeed straightforward exercise by any convenient method. For example, some of the classical quadrature rules (Gauss-Legendre, … ) can confidently be used in the dFPT which allows unequivocal determination of the integration limits.
As to the metabolite map, the primary focus on panels (b-e) is on the three recognized cancer biomarkers (Lac, Cho, PC) as indicated already in the first derivative |D 1 FPT| (b). However, this does not mean that the other metabolites are not diagnostically relevant. They are and, in fact, even some of the macromolecular resonances could add to the diagnostic value of in vivo MRS. This problem deserves attention and needs to be separately studied as we shall do in the nearest future.
Overlapped resonances coupled with relatively poor SNR represent one of the perennial stumbling blocks of data analysis and interpretation using time signals encoded from patients by in vivo MRS. They are abundantly present also in the total shape spectrum, the nonderivative envelope |FPT| (a). This is seen therein in several main compound peaks, e.g. the lactate doublet ( ∼ 1.30-1.40 ppm), GABA + NAA + NAAG ( ∼ 1.95-2.10 ppm), total creatine Cr + PCr ( ∼ 3.0-3.1 ppm), total choline Cho+PC ( ∼ 3.20-3.22 ppm), etc. (GABA = gamma amino butyric acid, NAAG = nitrogen acetyl aspartyl glutamic acid).
Such a difficult twofold problem is successfully and simultaneously solved by the dFPT through narrowing linewidths and flattening the background baseline. This process begins with the first derivative |D 1 FPT| (b). It continues with the second |D 2 FPT| (c) and the third |D 3 FPT| (d) derivatives. The culmination of such a steadily progressing trend is with the fourth derivative |D 4 FPT| (e). As in Figs. 2, 3 and 4, some ten most tightly packed resonances are marked by the red star symbols on panels (c) and (d) for |D 2 FPT| and |D 3 FPT|, respectively. They are all overlapped in |D 2 FPT| (c). However, six of them are separated in |D 3 FPT| (d).
Finally, the remaining four overlapped resonances around 4 ppm are split apart in |D 4 FPT| (e). As such, the envelopes |D 3 FPT| (d) and |D 4 FPT| (e) succeed in separating 50 and 54 resonances, respectively. In the examined, relatively small frequency band 0.35-4.535 ppm, the total number of the physical resonances, quantified in one go, is 54. The other physical resonances, temporarily immersed into the background baseline, coincident with the chemical shift axis, can also be faithfully reconstructed in the second go in the manner of sequential quantification.

Implications of the present findings for oncology: focus on brain tumor diagnostics
Primary cancerous brain tumors are quite rare in adults compared to other malignancies. However, they engender much attention related to the fear associated with the location, the young age at which these can appear, as well as their frequently poor prognosis [2]. Among children, primary brain tumors are the leading cause of solid tumor-related morbidity and mortality [39,40]. Molecular imaging through MRS is of potentially seminal importance for nearly all aspects of brain tumor diagnostics and management. This encompasses initial detection, tumor characterization and grading, treatment planning in radiation therapy, surgical guidance, as well as post-therapeutic follow-up, including identification of prognostic indicators [41][42][43][44][45][46]. It has been recently stated that a "new era" is foreseen, in which "spectrobiopsy will offer a vital new tool for the non-destructive differential diagnosis and characterisation of brain tumors" [47].
This prospect relies critically on the proper analysis and interpretation of MRS data. Encoded MRS times signal data are uninterpretable directly. Therefore, certain suitable mathematical transforms (Fourier, Padé, Laplace, … ) are needed for passing to the more manageable frequency domain. In the latter domain, peaks appear that in some circumstances could be associated with metabolites present in the excited slice of the tissue scanned by MRS. Within this framework, the limitations of the standard Fourier-based signal processing and fitting become all the more striking. Much of neurodiagnostics through FFT-based MRS has been reduced to a rough, semi-quantitative approach, based on a very small number of metabolites and their often ambiguously assessed concentration ratios. Numerous diagnostic dilemmas remain. phosphocholine, PC) clearly identified and quantified, but, for the first time, the abundance of PC can be easily appreciated on visual inspection of the envelope spectrum. Moreover, the vitally important PC to Cho concentration ratio can be straightforwardly and unequivocally assessed via nonparametric dFPT processing of in vivo encoded MRS time signals generated on a standard clinical MR scanner of magnetic field strength 1.5T.

Conclusion
This study is on derivative signal processing applied to in vivo magnetic resonance spectroscopy, MRS. The employed short time signals or FIDs (512 data points) have been encoded at a 1.5T clinical scanner from a brainstem glioma in a patient. The principal data analysis is carried out using the nonderivative and derivative fast Padé transforms, FPT and dFPT, respectively. The main emphasis is on the nonparametric dFPT. The parametric FPT and dFPT are also used, but only to cross-validate the corresponding Padé nonparametric reconstructions. Moreover, the nonparametric FPT and dFPT are compared with nonderivative and derivative fast Fourier transforms, FFT and dFFT, respectively.
Very recently [21][22][23], in derivative estimations aimed at proton MRS, we set up a comprehensive and robust validation procedure for the dFPT and dFFT, using the FIDs encoded from a test sample, the well-known Philips polyethylene phantom [32,33]. In computations, this procedure systematically involves all the variants of the FPT and dFPT (nonparametric, parametric, nonderivative, derivative) for the retrieved component spectra and envelopes. Further, the nonparametric or shape estimations by the FPT and dFPT have been compared with their Fourier counterparts, FFT and dFFT.
The FFT and dFFT can do only shape estimations. It was found in Refs. [21][22][23] that the dFFT breaks down for derivative estimations because it modifies the input FID in such way that mainly noise is sampled through the later acquired data points.
In sharp contrast, the overall performance of the nonparametric dFPT was outstanding for the polyethylene phantom [21][22][23]. In particular, the exemplary accomplishment of this latter variant of the dFPT was in terms of comparisons between the nonparametric and parametric estimations of various Padé spectra (components, envelopes) in the most useful, magnitude mode.
The same mentioned validation procedure [21][22][23] has also been applied to the present in vivo encoded FIDs. Moreover, all the acquisition parameters were identical in the present case of in vivo MRS and in the polyethylene phantom [21][22][23].
Therefore, the next critical step would be to see how much of this achievement can also be repeated for the present FIDs encoded by in vivo MRS. This type of inquiry is very important in view of the fact that the FIDs from the polyethylene phantom and from a patient are fundamentally different. The former has high concentrations of the constituent molecules and, hence, very good signal-to-noise ratio, SNR. Moreover, the content in the polyethylene phantom is known and the concentrations of the molecules therein are given. The present in vivo MRS time signals are of a relatively low SNR and, of course, the number and abundance of metabolites in the excited slice of the scanned human tissue are totally unknown prior to data analysis.
Despite such fundamental differences between the mentioned two kinds of the FIDs, it is gratifying to report that the net outcome of our current analysis completely matches every single facet of the expectation stemming from the conclusions of the polyethylene phantom. In particular, there is full agreement between the nonparametric and parametric Padé estimations of envelopes in all the nonderivative and derivative reconstructions. The most stringent of the effected cross-validations (all within Padé processings), was to compare the nonparametric envelopes with the components generated by parametric estimations. The obtained extraordinary result is best manifested in coincidence of some 50 and 54 resonances in the third and fourth derivative spectra, respectively, given by the lineshapes for the parametric components and nonparametric envelopes.
Abundant overlapped resonances superimposed on a noisy and elevated background baseline, as a common feature of all encoded in vivo MRS time signals, are also present in the standard nonderivative envelopes (be they from the FPT or FFT). The usual procedures of fitting such envelopes is fruitless for the welldocumented reasons. The dFFT even exacerbates this situation and generally may produce the envelopes that are worse than those in the FFT. By contrast, the nonparametric dFPT steadily improves the FPT with every derivative order. The improvement is in both resolving power and SNR. This is achieved in the dFPT, by lineshape narrowing and background flattening. The result is translated into the increased peak heights and the diminished baseline. Hence increased SNR.
Such a success of the dFPT in its role as a shape estimator cannot be understated given that precisely this type of performance is sought in medical diagnostic applications of in vivo MRS whose implementation in the clinic is basically precluded by overlapped, uninterpretable resonances. Not only that the nonparametric dFPT visualizes the numerous component resonances, hidden within some overlapped compound peaks, but it does so also for the diagnostically most relevant metabolites, the recognized cancer biomarkers (lactate, choline, phospocholine, … ).
This detection does not stop at mere visualizations. It is followed by providing the quantified in vivo MRS data, the concentrations of the metabolites assigned to the recovered resonance peaks. It is a veritable quantum leap by which an initial limitation to a purely qualitative estimation of the lineshape forms or profiles (nonparametric dFPT) subsequently turns into a method which secures clinically quantitative determination of the diagnostically most useful information (metabolite concentrations, chemical shifts, spin-spin relaxation times, … ). We then have here offered a theoretical framework of the proven validity, the nonparametric derivative fast Padé transform, dFPT, as the method of choice at service to dNMR spectroscopy for applications in analytical chemistry, cancer medicine and many other fields relying upon signal processing.
Funding Open access funding provided by Karolinska Institute.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.