1 Introduction

Magnetic resonance spectroscopy (MRS) is one of the most versatile methods for uncovering the internal structure of matter [1,2,3]. For example, MRS can find out how many molecules, and of which specific kind, are present in vastly different specimens/substances irrespective of their structural complexity. Such quantitative information is of utmost importance because the detected larger or smaller concentrations (abundances) of the most informative molecules, when compared with the corresponding standards regarding the matter under investigation, can advantageously be used for many purposes not only in basic and applied sciences but also in technology and industry.

There are two domains associated with MRS. These are the time domain and the frequency domain. The results from MRS can be given in either form of this dual representation. They are equivalent since they contain the same information. The passage from one to the other representations is achieved by means of the direct and inverse transforms (Fourier, Padé).

Being primarily a measuring/detecting device, MRS encodes time signals (in the time domain) stemming from the intrinsic molecules of the examined substance. These time signals, or equivalently, free induction decay (FID) curves come after an external perturbation of a selected slice of the specimen by the triple magnetic perturbations (static as well as gradient magnetic fields and radiofrequency pulses). With the act of encoding, MRS completes its measuring task. However, MRS does not itself provide the key information which is the sought assembly of molecular concentrations in the specimen.

Theory (by way of mathematical methods) does so by data processing of time signals acquired by MRS [2, 3]. Data analysis by mathematical methods can be carried out in either the time or frequency domains. This is opposed to encoding which is only in the time domain. The literature in MRS is abundant with the terms “acquired spectra”. This is wrong, since spectra are never encoded by MRS. Spectra are the results of the theory by which various mathematical transforms are used to pass from the acquired time signals to the computed spectra.

One of such mathematical transforms in MRS is the fast Padé transform (FPT) [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. It can perform both shape and parameter estimations. The customary form of either of these two estimations does not take derivatives of the computed spectra and, thus, it is referred to as the non-derivative (or zero-order derivative) FPT.

The most recent upgrade of the FPT is the derivative fast Padé transform (dFPT) [6, 15,16,17]. This processor can provide quantification by exclusive reliance upon shape estimation. In other words, the non-parametric dFPT deals with spectral envelopes alone, and yet its high-order derivatives can reconstruct the precise numerical values of the number of the physical peaks, their positions, widths, heights, areas (and, thus, molecular concentrations) without recourse to any fitting.

First introduced and benchmarked on synthesized MRS time signals, the dFPT showed its outstanding performance by simultaneously eliminating noise and separating tightly overlapping resonances. For instance, despite being merely 0.001 ppm apart, metabolites phosphocholine (PC) and phosphoethanolamine (PE), have been separated and exactly quantified by the non-parametric dFPT [15,16,17] (hereafter, ppm denotes parts per million). In the reconstructed total shape spectrum or envelope, these two resonances were overlapped to such an extent that they appeared as a single perfectly symmetrical Lorentzian. This envelope gave no hint whatsoever that underneath of a large PE peak resided a small PC peak. The concentrations of the PC and PE molecules were 0.3 and \(2.25~\upmu \hbox {M/g}\), respectively. This was a challenging quantification problem which the non-parametric dFPT solved readily and exactly by computation of higher derivatives of the pertinent total shape spectrum.

The gist of the matter is that the derivative operator of a sufficiently high order is, in fact, a converter which maps the magnitude envelope from the dFPT into a set of well-isolated (i.e. non-overlapping) stable peaks. Most importantly, when appropriately scaled, such lone peaks turned out to be the corresponding component spectra for each physical resonance contained in the envelope from the non-derivative parametric FPT.

These isolated resonances are physical (genuine) because all the unphysical (spurious, noisy) resonances are flattened down to their zero-valued heights (i.e. annulled) since derivatives effectively suppress noise in the dFPT. Stated equivalently, the derivative operator (when teamed up with the Padé-based estimation) is an efficient noise filter.

The reason for which the scaled isolated peaks from shape estimation by the dFPT (multiple derivatives) are the component spectra from parameter estimation by the FPT (no derivatives) is twofold:

  1. (i)

    The same high-order derivatives applied to both the envelope in the non-parametric dFPT and the component spectra in the parametric dFPT coincide.

  2. (ii)

    The peak parameters in the component spectra from the two parametric estimations, dFPT (multiple derivatives) and FPT (no derivatives), are uniquely related by the known analytical expressions.

Therefore, it follows from (i) and (ii) that the peak parameters in the isolated resonances from the given envelope in the non-parametric dFPT are uniquely related to the peak parameters in the component spectra from the parametric FPT.

The reconstruction linelist of the applied non-parametric dFPT, as the output of the computer program, contains the number of the peaks, their positions, widths and heights for each physical resonance in the magnitude mode of the high-order derivative envelopes. All that is needed is to compute these envelopes at a very dense grid of the sweep frequencies. Then, in the program, the command “max” for the values (sticks) of the magnitude envelope in the chosen isolated Lorentzian would simultaneously give the peak height and the chemical shift (resonance frequency). Moreover, within the same magnitude envelope of this Lorentzian, the output of the program also contains the two numerical values (alongside their positions on the chemical shift axis) of the halved maximum stick. The difference between these two latter chemical shifts is the full-width at half-maximum (FWHM). The peak width and the peak height yield the peak area which is proportional to the concentration of the molecule (metabolite) to which the given resonance is assigned. No peak phase is required to be reconstructed since the magnitude mode is phase-insensitive. The same procedure is performed automatically for all the physical (genuine) peaks in the magnitude spectrum. The unphysical (spurious) resonances are washed out.

This expounded procedure amounts to providing quantification without actually solving the spectral analysis problem (and, of course, without fitting, either). In Refs. [15,16,17], this goal was accomplished, with both noiseless and noisy simulated MRS time signals. After a thorough benchmarking of the dFPT for synthesized input data, it would be important to extend the analysis to encoded MRS time signals. The present study is the first step in this direction. Our goal here is to see whether already low derivative orders (from the 1st to the 4th) in the non-parametric dFPT could simultaneously improve frequency resolution and signal-to-noise ratio (SNR). The subsequent studies in our plan will use higher-order derivatives in the non-parametric dFPT to test explicit quantification of encoded MRS time signals by performing only shape estimations.

Very recently [19] it has been demonstrated that, regarding SNR, the FPT can dramatically improve acquisition of MRS time signals. The two main routine pulse sequence encodings used in acquiring time signals in MRS are the stimulated echo acquisition mode (STEAM) and the point-resolved spectroscopy sequence (PRESS). Overwhelmingly, PRESS has been favored in practice because of its better SNR relative to STEAM. This has occurred despite several important advantages of STEAM compared with PRESS. The question asked in Ref. [19] was: is it possible to exploit the known advantages of STEAM by performing FPT-guided encoding instead of the conventional fast Fourier transform (FFT)-assisted acquisition of FIDs from MRS? The favorable SNR capabilities of the FPT motivated this question. The answer from Ref. [19] was in the affirmative: noise root-mean-square (RMS) errors in the spectra were decreased by orders of magnitudes (between 3619 and 14253%) relative to the FFT. Moreover, the SNR of creatine (Cr) and choline (Cho) peaks reconstructed by the FPT were increased more than 95%. Since noise is suppressed even more by the derivative fast Padé transform than by the FPT, it is foreseen that further gains in SNR in conjunction with STEAM can be achieved by implementing the dFPT-based acquisition of time signals encoded by MRS.

2 Recapitulation of benchmarking the dFPT on synthesized MRS time signals

The dFPT is the \(m\hbox {th}\) order \((m=1,2,3,\ldots )\) derivative of the conventional FPT \((m=0)\). The derivative operator \(D_\nu \) is taken with respect to the sweep linear frequency \(\nu \) as \(D_\nu =\text {d}/\text {d}\nu \). It is applied to the customary non-derivative complex envelope \(P_K/Q_K\) in the FPT, where \(P_K\)and \(Q_K\) are polynomials. Here, K is the model order which also represents the number or resonances, i.e. the number of molecules in the substance under study. Both the FPT and dFPT can be processed non-parametrically, i.e. by shape estimations.

The introduction of the dFPT took place in Refs. [15,16,17] for synthesized MRS time signals with zero-valued phases. Therein, the primary focus was on the critical challenge of identifying and quantifying phosphocholine, PC, in a manner which would be the most amenable to unambiguous interpretation. Noiseless data were tested first. Afterward, we proceeded to noise-corrupted time signals. These simulated data were generated to resemble typical encoded MRS time signals that lead to both isolated and tightly overlapped resonances. In MRS, reliable identification and quantification of PC is a exceedingly demanding task for non-parametric signal processors.

The PC resonance is completely invisible in customary non-derivative envelopes. This is due to the dominant adjacent PE resonance under which PC is entirely buried. As stated, the PC-to-PE chemical shift separation is only 0.001 ppm with PC and PE resonating at 3.220 and 3.221 ppm, respectively. In general, the PC and PE metabolites sharply differ in abundance. In such a case, a practical goal with derivatives is to unravel this overlap, i.e. the underlying structure, by progressively narrowing the resonances. Thereby, eventually each resonance becomes clearly visible and its abundance unequivocally displayed at higher derivative orders.

Within the FPT, in the non-derivative absorption envelope \({\mathrm{Re}}(P_K/Q_K)\), the compound resonance near 3.220 and 3.221 ppm is a completely symmetrical single peak with no hint whatsoever that there could be any underlying structure. However, as seen in Refs. [15,16,17], the situation changed subtly in the magnitude mode with application of the 1st derivative of the Padé complex envelope, \(|(\text {d}/\text {d}\nu )P_K/Q_K|\). Therein, the peak at 3.220 ppm was no longer symmetrical, since a slight bulge emerged to the right. Further, already with the 2nd derivative in the magnitude mode \(|(\text {d}/\text {d}\nu )^2P_K/Q_K|\), this bulge became more distinct, suggesting the possibility of another peak besides the dominant PE resonance. The situation became even clearer for the 3rd derivative \(|(\text {d}/\text {d}\nu )^3P_K/Q_K|\) where PC was plainly visible. Additionally, the baseline was flattened throughout the entire displayed chemical shift region from 3.21 to 3.29 ppm.

In Refs. [15,16,17], the detailed examination of increasing derivative orders has been presented in full. Therein, we proceeded to the high order derivatives where complete stabilization was observed. Namely, in Ref. [15], we computed the 42nd, 46th and 50th derivative orders via \(|(\text {d}/\text {d}\nu )^{42}P_K/Q_K|\), \(|(\text {d}/\text {d}\nu )^{46}P_K/Q_K|\) and \(|(\text {d}/\text {d}\nu )^{50}P_K/Q_K|\), respectively. Overall, the full converging effect of the derivative was seen. When the saturation level was reached, the peak widths became systematically narrowed by the pertinent identifiable scaling factors that made it possible to exactly reconstruct the peak parameters. This includes the critical region 3.220–3.221 ppm containing the PC and PE peaks. The two latter overlapping resonances were extremely well delineated, with both of their baselines descending nearly or completely to the zero value of the ordinates (i.e. merged into the abscissae).

The tails of all the other derivative lineshapes were entirely embedded in the chemical shift axis itself. This shows that the differential operator in \(|(\text {d}/\text {d}\nu )^mP_K/Q_K|\), no matter how high the derivative order m is taken, introduces no noise, random or any other kind for that matter (from e.g. round-off errors in numerical computations with finite arithmetic). Such findings suggested that the dFPT could also perform well for noisy time signals reaching the full converging effect of the derivative transform in the so-called “normalized” representation.

As shown in Ref. [15], normalization in the dFPT has been performed, so that the increased peak heights of derivative spectra \((m>0)\) could still be plotted on the same graphs with the non-derivative \((m=0)\) envelopes. Crucially, the cross-checking relationships of the lineshapes for the magnitude modes with the real absorptive and imaginary dispersive parts of the given complex-valued envelopes have been established [16].

In non-derivative spectra \((m=0)\), the absolute value (modulus) of the amplitudes \(d_k\) from an MRS time signal, as well as the resonance peak positions, are the same in the magnitude and absorption modes. On the one hand, the peak widths are broadened by a factor of \(\sqrt{3}\) in the magnitude mode relative to the absorption mode \((m=0)\). In Ref. [16], we provided an in-depth comparison of spectra of increasing derivative orders for magnitude, absorption and dispersion modes. Therein, at higher derivative orders m, the resonances appear as thin, distinct peaks in the magnitude mode, whereas in the absorption or dispersive modes, the resonances have several side lobes. Thus, the magnitude mode of higher order derivative spectra is more conducive to interpretation than the absorption mode of derivative spectra.

As such, the extremely narrow symmetric Lorentzian peaks (from their tips to bottoms) in the magnitude modes for higher derivative order solve a practical problem encountered with conventional non-derivative \((m=0)\) lineshapes. Namely, extended Lorentzian tails always mask adjacent and even distant spectral structures of potentially informative content. This is recognized as the “leakage problem” in MRS. The high-order derivative envelopes circumvent this problem by nullifying the tails and exactly determining the peak areas, thus yielding the sought correct concentrations of metabolites in the scanned substance.

After validation for the noiseless case, we proceeded to examine the performance of the dFPT for time signals with added noise, and compared this with the derivative fast Fourier transform (dFFT), using the absorption mode. This was done in Refs. [6, 15]. Both processors performed envelope estimations using the same noisy time signal points \(\{c_n+r_n\}\) where \(\{c_n\}\) are the noiseless time signal points. The added corruptions \(\{r_n\}\) were a set of random numbers with zero-mean Gaussian distributions (white noise) of standard deviation \(\sigma = 0.0289\) RMS. These simulated time signals were based upon those encoded in Ref. [20]. The total signal length of these encoded time signals was very long (62 MB or 65536 data points). Such a long signal length was needed in Ref. [20] in order to obtain reasonable appearing envelopes in the conventional FFT \((m=0)\).

However, for the simulated noiseless time signals, there was no reason to use all the 65536 data points in the Padé-based reconstructions. The same rationale also applies to the noisy synthesized time signals for which we used a much shorter signal length in the FPT and dFPT. As illustrations, comparisons were made employing the results of the reconstructions by the dFPT and the dFFT for the derivative orders \(m=3-6\). The dFPT yielded results that were the same for the noiseless and noisy time signals. This implies noise suppression by the 3rd to 6th order derivative transform in the dFPT. On the other hand, the same added noise severely distorted the 3rd to 6th order derivative spectral envelopes in the dFFT.

Even though the noise corruption level of the input time signal was of a rather small standard deviation (\(\sigma = 0.0289\) RMS), it is remarkable that this noise was hugely amplified in the dFFT by proceeding from the third to the sixth-order derivatives. For higher-orders m of derivatives \((m\ge 7)\), the dFFT completely fails, as it processes merely the enhanced noise. By contrast, the results from the dFPT remain the same for noise-free and noise-corrupted time signals in the cases with \(m\ge 7\), as was also for \(0\le m\le 6\). This demonstrates the noise suppression capability of the \(m\hbox {th}\) derivative operator \(D^m_\nu \) when applied to the envelopes from the non-parametric FPT.

The dFFT greatly amplifies noise even for relatively low derivative orders m, such that all genuine information is lost. This occurs because the dFFT processes the product of the power function \(t^m\) and the given time signal c(t), both in their digitized forms. Here, the \(m\hbox {th}\) power of time t puts higher weight on the later time signal points, dominated by noise in encoded data.

In Ref. [17], we made one critical step further, by setting up the goal of validating the non-parametric dFPT with the parametric dFPT. This was done by comparing the full lineshapes of derivative envelopes from the non-parametric dFPT with the corresponding derivative component spectra from the parametric dFPT.

As a shape estimator, the non-parametric dFPT never even addresses, let alone solves the quantification problem (i.e. it never consider eigenvalue problems, rooting characteristic/secular equations, ...). The parametric dFPT first solves the quantification problem from which the lineshapes of components and envelopes are plotted. Thus, if the derivative component spectra from the parametric dFPT could be fully reconstructed by the derivative envelopes from the nonparametric dFPT, the goal of achieving quantification would be done by derivative lineshape processing alone. This was indeed possible, as has been shown in Ref. [17].

3 Theory of derivative signal processing

3.1 Derivative fast Fourier transform

For a continuous time signal c(t), the standard finite Fourier integral denoted by \(F(\nu )\) over time t, which runs from zero to T, is given by:

$$\begin{aligned} F(\nu )=\frac{1}{T}\int \limits _0^T\text {d}t\mathrm{e}^{-2\pi i\nu t}c(t). \end{aligned}$$
(3.1)

Here, T is the total duration of the time signal c(t). If time t is equidistantly discretized (digitized) as \(t=n\tau \, (n=0,1,2,\ldots ,N-1)\), the Fourier integral \(F(\nu )\) from (3.1) would become the Riemann sum:

$$\begin{aligned} F(\nu )=\frac{1}{N}\sum \limits _{n=0}^{N-1}\mathrm{e}^{-2\pi i\nu n\tau }c_n, \end{aligned}$$
(3.2)

where \(c_n\equiv c(n\tau )\) and \(T=N\tau \). Here, \(\tau \) is the sampling time. In an MRS encoded time signal, T is the total acquisition time of all the N data points in a single FID set \(\{c_n\}\, (0\le n\le N-1)\). Frequency \(\nu \) is a continuous (analog) variable in (3.2). If \(\nu \) is also digitized by sampling it at the N equidistant points (the so-called Fourier grid), the Riemann sum on the rhs of Eq. (3.2) would be reduced to the discrete Fourier transform (DFT) denoted by \(F_k\) as:

$$\begin{aligned} F_k=\frac{1}{N}\sum \limits _{n=0}^{N-1}\mathrm{e}^{-2\pi ink/N}c_n. \end{aligned}$$
(3.3)

For \(N=2^\ell \,(\ell =1,2,3,\ldots )\), computations in the DFT with the \(N^2\) multiplications is accelerated to only \(N\log _2N\) multiplication using the Tukey–Cooley fast algorithm and this leads to the FFT [2].

If a derivative operator \(D^m_\nu \) of order m with respect to \(\nu \) via \(D^m_\nu =(\text {d}/\text {d}\nu )^m\) is applied to \(F(\nu )\) from Eq. (3.1), it would follow:

$$\begin{aligned} F^{(m)}(\nu )=\frac{1}{T}\int \limits _0^T\text {d}t\mathrm{e}^{-2\pi i\nu t}c_m(t),\quad c_m(t)=(-2\pi it)^mc(t), \end{aligned}$$
(3.4)

where \(F^{(m)}(\nu )\equiv D^m_\nu F(\nu )\) and,

$$\begin{aligned} D^m_\nu =\left( \frac{\text {d}}{\text {d}\nu }\right) ^m\quad (m=1,2,3,\ldots ). \end{aligned}$$
(3.5)

Similarly, applying \(D^m_\nu \) to Eq. (3.2) gives the derivative discrete Fourier transform (dDFT):

$$\begin{aligned} F^{(m)}_k=\frac{1}{N}\sum \limits _{n=0}^{N-1}\mathrm{e}^{-2\pi ink/N}c_{n,m},\quad c_{n,m}=(-2\pi in\tau )^mc_n. \end{aligned}$$
(3.6)

We then see that the derivative operator \(D^m_\nu \) in Fourier-based MRS is equivalent to introducing the apodization function \((-\,2\pi it)^m\) or \((-\,2\pi in\tau )^m\) multiplying the analog or digital time signal c(t) or \(c_n\), respectively. This apodization puts emphasis on the noise-dominated tails of the encoded times signals. The derivative fast Fourier transform, dFFT, is obtained from the total shape spectrum \(F^{(m)}_k\) by using the Tukey-Cooley algorithm for the dDFT in Eq. (3.6). This amounts to simply multiplying the input time signal \(\{c_n\}\, (0\le n\le N-1)\) by the apodization power function \((-\,2\pi in\tau )^m\) and applying the standard non-derivative FFT \((m=0)\) to the product \((-\,2\pi in\tau )^mc_n\) for a fixed derivative order \(m>0\). Simple, it undoubtedly is. But that is about all. The price for this simplicity is, however, very high: even for relatively low values of m, the power function \(\sim (n\tau )^m\) multiplying \(c_n\) takes over to become a dominant factor at longer FID encoding. As a result, the dFFT ends up by processing merely noise from encoded FIDs with all the physical information irretrievably lost.

3.2 Derivative fast Padé transform

The dFPT consists of applying the multiple-derivative operator \(D^m_\nu =(\text {d}/\text {d}\nu )^m\) of order \(m (m=1,2,3,\ldots )\) to the given total shape spectrum from the non-derivative conventional FPT \((m=0)\). For the dFPT, the input envelopes can be from the non-derivative (parametric or non-parametric) FPT. However, the non-parametric dFPT is of the main interest as the goal is to show that quantification can be achieved by exclusive reliance upon shape estimation. The parametric FPT accomplishes quantification already in its non-derivative variant \((m=0)\). Therefore, the parametric dFPT was used only at the benchmarking stage to validate the non-parametric dFPT.

As is well-known, the complex-valued envelopes, i.e the Green response functions \(G^\pm _K(z^{\pm \,1})\), in the two equivalent variants of e.g. the diagonal non-derivative fast Padé transform, denoted by \({\mathrm{FPT}}^{(\pm )}\), are defined by the rational functions:

$$\begin{aligned} {\mathrm{FPT}}^{(\pm )}{:}\quad G^\pm _K(z^{\pm \,1})=\frac{P^\pm _K(z^{\pm \,1})}{Q^\pm _K(z^{\pm \,1})},\quad z=\mathrm{e}^{\pm 2\pi i\nu \tau }, \end{aligned}$$
(3.7)

where \(z^{\pm \,1}\) are the complex harmonics, \(\nu \) is an arbitrary linear frequency (real or complex) and, as before, \(\tau \) is the constant sampling rate for the input time signal \(\{c_n\}\) of full length \(N\, (0\le n\le N-1)\). Here, \(P^\pm _K(z^{\pm \,1})\) and \(Q^\pm _K(z^{\pm \,1})\) are the characteristic polynomials (with a common degree K), given in terms of their expansion coefficients \(\{p^\pm _r,q^\pm _s\}\), respectively:

$$\begin{aligned} P^\pm _K(z^{\pm \,1})=\sum \limits _{r=r^\pm }^Kp^\pm _rz^{\pm r},\quad Q^\pm _K(z^{\pm \,1})=\sum \limits _{s=0}^Kq^\pm _sz^{\pm s}, \end{aligned}$$
(3.8)

where \(r^+=1\) and \(r^-=0\). The superscripts ± are employed to indicate the use of \(z^1\equiv z\) and \(z^{-1}\) as the independent variables in the spectral functions \(G^\pm _K(z^{\pm \,1})\). For \(\nu \) complex, to have damped harmonics, we must have \(|z^{\pm \,1}|<1\). Thus, in the case of complex \(\nu \) with \({\mathrm{Im}}(\nu )>0\), the \({\mathrm{FPT}}^{(+)}\) and \({\mathrm{FPT}}^{(-)}\) converge for \(|z|<1\) and \(|z|>1\), i.e. inside and outside the unit circle, respectively. However, by way of the Cauchy analytical continuation, the \({\mathrm{FPT}}^{(+)}\) and \({\mathrm{FPT}}^{(-)}\) also converge outside and inside the unit circle \(|z|>\) and \(|z|<1\), respectively. In other words, both the \({\mathrm{FPT}}^{(+)}\) and \({\mathrm{FPT}}^{(-)}\) converge everywhere in the complex plane of the harmonic variables z and \(z^{-1}\), respectively, except at the poles \(z^{\pm \,1}_k\) of \(G^\pm _K(z^{\pm \,1})\). The poles of \(G^\pm _K(z^{\pm \,1})\) are the K roots of the denominator polynomials \(Q^\pm _K(z^{\pm \,1})\), i.e. the unique solutions \(\{z_k^{\pm \,1}\}\, (1\le k\le K)\) of the characteristic equations \(Q^\pm _K(z^{\pm \,1})=0\).

As mentioned, we obtain the \({\mathrm{dFPT}}^{(\pm )}\) of order m by applying the operator \(D^m_\nu \) onto the \({\mathrm{FPT}}^{(\pm )}\). Introducing the notation \(G^{(m)\pm }_K(z^{\pm \,1})\) for the \(m\hbox {th}\) derivative of \(G^{\pm }_K(z^{\pm \,1})\) with respect to \(\nu \) as:

$$\begin{aligned} G^{(m)\pm }_K(z^{\pm \,1})=D^m_\nu G_K(z^{\pm \,1}), \end{aligned}$$
(3.9)

the \({\mathrm{dFPT}}^{(\pm )}\) of order m becomes

$$\begin{aligned} {\mathrm{dFPT}}^{(\pm )}{:} \quad G^{(m)\pm }_K(z^{\pm \,1})=D^m_\nu \left( \frac{P^\pm _K(z^{\pm \,1})}{Q^\pm _K(z^{\pm \,1})}\right) . \end{aligned}$$
(3.10)

We emphasize that in Eq. (3.10), the input envelopes \(G^\pm _K(z^{\pm \,1})\), i.e. the polynomial ratios \({P^\pm _K(z^{\pm \,1})}/{Q^\pm _K(z^{\pm \,1})}\) from the customary non-derivative \({\mathrm{FPT}}^{(\pm )}\) are computed non-parametrically. Advantageously, the non-parametric derivative envelopes \(D^m_\nu (P^\pm _K(z^{\pm \,1})/Q^\pm _K(z^{\pm \,1}))\) from Eq. (3.10) provide the exact number of peaks and their parameters (positions, widths, heights and phases) of all the physical resonances. This has explicitly been shown in Ref. [17] for large m by the complete agreement between the lineshapes of the non-parametric derivative envelopes and the derivative component spectra in the \({\mathrm{dFPT}}^{(-)}\) of order m. The latter spectra refer to the lineshapes obtained by applying the derivative operator \(D^m_\nu \) to the component spectra constructed after solving the quantification problem in the parametric \({\mathrm{FPT}}^{(-)}\). Although not reported in Ref. [17], the same findings and conclusions have been verified to hold true for the high-orders m in the \({\mathrm{dFPT}}^{(+)}\), as well.

Relationships exist between the two sets of the peak parameters, one for the \({\mathrm{dFPT}}^{(\pm )}\) with \(m>0\) and the other for the \({\mathrm{FPT}}^{(\pm )}\) with \(m=0\) [17]. Thus, reliance exclusively upon the non-parametric derivative envelopes, permits reconstruction of the exact peak positions, widths, heights and phases of every physical resonance that would have been available from a direct application of the non-derivative parametric \({\mathrm{FPT}}^{(\pm )}\).

The spectra in the non-parametric and parametric \({\mathrm{dFPT}}^{(\pm )}\) are given by the analytical expressions for any arbitrary positive integer m in terms of the Bell and Eulerian polynomials, respectively. With these analytically available derivative spectra, the operator \(D^m_\nu \) does not introduce any additional noise when the input time signal \(\{c_n\}\) is processed by the \({\mathrm{dFPT}}^{(\pm )}\). Note also that analytical (as opposed to numerical) differentiation is itself free from errors. Moreover, in the \({\mathrm{dFPT}}^{(\pm )}\), the derivative operator \(D^m_\nu \) never acts on the \(\nu \)-dependent harmonic variable \(\exp {(-\,2\pi i\nu t)}\) which multiplies the time signal c(t), as opposed to the passage from (3.2) to (3.4) when going from the DFT to dDFT (and afterward from the FFT to dFFT), respectively.

This means that in the \({\mathrm{dFPT}}^{(\pm )}\) there is no dFFT-type apodization by which the power function \(\sim (n\tau )^m\) multiplies the time signal \(c_n\) (nor any other apodization for that matter). As noted, when the time signal number n (in units of the sampling rate \(\tau \)) is augmented, as encoding proceeds toward the completion of the total acquisition time T, more and more noise is encountered. In other words, this power function massively counters the usual exponential decay of \(c_n\). As a result, especially for higher derivative orders m, mainly noise is processed when the FFT is applied to the product \(\sim (n\tau )^mc_n\) in the course of obtaining the envelopes in the dFFT.

According to the theory of derivative resonances in derivative envelopes computed by the \({\mathrm{dFPT}}^{(\pm )}\), when the derivative order m is augmented, the peak positions do not change, whereas the peak heights increase and simultaneously the peak widths decrease. These latter two trends for \(m>0\) are concomitant, such that their product in a normalized derivative resonance can reconstruct the peak area of the absorption mode of the non-derivative \((m=0)\) component spectrum from the parametric \({\mathrm{FPT}}^{(\pm )}\).

Importantly, and this is the main advantage of derivative estimation by the non-parametric \({\mathrm{dFPT}}^{(\pm )}\), for higher derivative orders m, the normalized derivative envelopes coincide with the Padé non-derivative component spectra. This key feature amounts to performing quantification by shape estimation alone. Therefore, particularly for the magnitude mode of the normalized high-order derivative envelopes in the \({\mathrm{dFPT}}^{(\pm )}\), the stabilized (relative to both K and m) estimation permits accurate determination of the positions, widths and heights of the non-derivative absorptive peaks. The resulting products of the peak heights with the peak widths are proportional to the concentrations of molecules in the substance examined by MRS. As mentioned, the magnitude mode is convenient because it is phase-insensitive, i.e. it does not necessitate phase corrections of the encoded time signals.

Moreover, peaks in the magnitude mode are straightforward for interpretation of the peak areas that are the essential ingredient to the sought molecular concentrations. These latter quantities are the ultimate goal of MRS. In non-derivative estimation, the magnitude mode is not used because, as noted, it widens the peak widths by \(\sqrt{3}\). However, in derivative estimation, especially for higher orders m, peak widths are significantly narrower than those in the non-derivative absorptive resonances. Thus, the magnitude mode is the spectral mode of choice for displaying the normalized derivative envelopes. On the other hand, the product of the peak heights with the peak widths remains constant once convergence has been attained. In order to facilitate monitoring of convergence of resonances from derivative envelopes, by comparing them with their non-derivative counterparts, it is necessary that the former envelopes are normalized relative to e.g. the tallest peak in the non-derivative envelopes from the \({\mathrm{FPT}}^{(\pm )}\). With such normalization procedures, the increased peak heights in derivative envelopes for \(m>0\) can be plotted on the same graphs with the non-derivative envelopes \((m=0)\). In the present work, all the Fourier and Padé magnitude derivative envelops are normalized to the NAA peak height (at \({\sim }2\) ppm) in the corresponding magnitude non-derivative envelopes.

4 Results and analysis of derivative reconstructions

4.1 Characteristics of the encoded time signals and the related reconstructed spectra

As described earlier in Ref. [18], the time signals were encoded at a General Electric magnetic resonance (MR) clinical scanner with the static magnetic field of relatively low strength, \(B_0=\) 1.5 T, associated with the Larmor frequency \(\nu _\mathrm{L}=63.87\,\mathrm{MHz}\). The encoding was performed at the Astrid Lindgren’s Children Hospital (Stockholm) from the parietal temporal brain region from a child who suffered cerebral asphyxia. Single-voxel proton MRS was used together with the PRESS sequencing. The number of excitations (NEX) was 128. These 128 encoded FIDs are averaged (arithmetic average) to improve SNR of the data. It is such an averaged time signal that is processed.

The averaged FID is phase corrected through multiplication by \(\exp {(i\phi _0)}\), where \(\phi _0=-1.1831\,{\mathrm{rad}}\). This is inconsequential for spectra in the magnitude mode for which, evidently, \(|\exp {(i\phi _0)}|=1\). Phasing is often used in MRS in an attempt to roughly bring the real part of the given complex spectrum closer to a positive-definite lineshape. As stated, phasing is not necessary for the magnitude mode of the envelope. The bandwidth (BW) and the sampling time \(\tau \) were 1000 Hz and 1 ms (\(\tau =\) 1/BW), respectively. The echo time (TE) and the repetition time (TR) were 24 ms and 2000 ms, respectively. The total length of each of the encoded FIDs was 512, implying that the total acquisition time \((T=N\tau )\) was 512 ms. Afterward, the averaged FID was zero-filled once to have the full signal length \(N=1024\) for both Fourier and Padé estimations.

Zero-filling in the given input FID can occasionally be beneficial esthetically to Fourier envelopes. It is usually not done in a Padé-based estimation because spectra in this processor can be computed at any desired number of frequencies. Moreover, for Padé signal processing, the FID can meaningfully be extended by adding an arbitrary number of realistic time signal points generated from the encoded data points themselves. Such a gain is achieved automatically by convolution of the encoded FID data points with the expansion coefficients of the denominator polynomial in the fast Padé transform. This gives a predictive power to Padé processing. Thus, the fast Padé transform can predict the FID data points that would have been measured had the encoding continued beyond the original total acquisition time (\(T=512\) ms in the case under study). Presently, we did not do this Padé-conceived extrapolation of the encoded 512 FID data points. Instead, for a direct comparison, we consider the same input time signal (the encoded 512 FID points plus 512 zeros) in both Fourier and Padé reconstructions.

Very recently [18], the same above-described FID has also been used in the non-derivative fast Padé transform, \({\mathrm{FPT}}^{(+)}\), within both parametric and non-parametric processing. In Ref. [18], the concept of iterative spectra averaging from Refs. [4, 6] was employed. The total number of iterations in averaging was 12 and 3 of them have been explicitly shown and analyzed in detail. Convergence of iterations was remarkably good such that even the 1st iteration could suffice for practical purposes. In the current study dealing with the first application of the \({\mathrm{dFPT}}^{(+)}\) to encoded FIDs, we will process the raw time signal alone (zero filled once). However, in the nearest future, the derivative fast Padé transform will be employed together with iterative spectra averaging (frequency domain) and extrapolation (time domain).

Presently, only non-parametric estimations will be applied by means of the \({\mathrm{FPT}}^{(+)}\) and \({\mathrm{dFPT}}^{(+)}\). For these two processors, the full signal length will be used \((N=1024)\) whose first 512 data points were encoded and the rest is a set of 512 added zeros (one zero filling). With such a choice, in this study, the model order K in the Padé spectrum \(P_K/Q_K\) is chosen to be 512 \((K=N/2=1024/2)\). As to our Fourier-based computations, besides \(N=1024\) (one zero filling), longer signal lengths will also be employed with repeated zero-padding to illuminate its net effect on the resolution in the dFFT.

In panel (a) of Fig. 1, the real part of complex spectrum, via the non-derivative \((m=0)\), non-parametric envelope \({\mathrm{Re}}(P^+_K/Q^+_K)\), is displayed. Therein, we see a substantial number of metabolite peaks identified according to their resonant frequencies, in the frequency range of interest (ROI) from 0.75 to 4.4 ppm. Notably, a number of peaks are split apart due to destructive interference. Thus, for example, N-acetyl aspartate (NAA) and N-acetyl aspartyl glutamic acid (NAAG) can be identified as two ridges on the top of the peak centered at about 2.0 ppm. Also noteworthy is that especially from about 1.3 to 3.8 ppm, most of the baseline is lifted above the zero level of the ordinate which is given in arbitrary units (au).

Panel (b) displays the non-derivative \((m=0)\) non-parametric envelope generated by the \({\mathrm{FPT}}^{(+)}\) in the magnitude mode \(|P^+_K/Q^+_K|\). Although there are substantial similarities between the two modes, there are also some noteworthy differences. For example, the peak centered at about 2.0 ppm on panel (b) is now nearly twice taller compared to that on panel (a), but without the clear ridges distinguishing NAA and NAAG. This peak at 2.0 ppm is amplified due to constructive interference in the magnitude mode, as opposed to destructive interference (around NAA and NAAG on panel a) yielding a dip in the real part of the complex envelope.

Moreover, at the outermost parts of the ROI, the envelope is markedly lifted. Thus, from about 0.9 to 0.75 ppm, the spectrum rises to about 13 au in the magnitude mode, whereas the real part of the complex spectrum descends to below 0 au. At the other extreme of the ROI, from about 4.2 to 4.4 ppm, the magnitude envelope ascends to nearly 15 au, while the real part of the complex envelope lowers to 4–5 au. The reason for these differences in the extrema of the ROI is related to the background, which is strongly influenced by the residual tail from water (presently set to resonate at 4.61 ppm) and lipids. The residual water peak (not shown) is still about 7 times taller than the NAA peak on panel (b) of Fig. 1.

The great majority of spectral peaks on panels (a) and (b) of Fig. 1 are seen not to be pure Lorentzians (bell-shaped symmetric lineshapes). The presence of non-zero phases of oscillatory waveforms in any encoded FID would lead to mixing of absorptive and dispersive profiles even for well-isolated resonances. Non-zero phases of the present FID are the reason for which, in Fig. 1a, the real part \({\mathrm{Re}}(P^+_K/Q^+_K)\) of the complex spectrum \(P^+_K/Q^+_K\) contains mainly non-Lorentzians. On the other hand, the magnitude mode is phase insensitive, in which case all isolated resonances should be pure Lorentzians. Nevertheless, this is still not the case for most peaks in the magnitude envelope \(|P^+_K/Q^+_K|\) from Fig. 1b because therein nearly all resonances are overlapped, not isolated.

Fig. 1
figure 1

Two spectral modes (real part vs. magnitude) in the \({\mathrm{FPT}}^{(+)}\) at the signal length \(N=1024\, (K=512)\). Panels (a, b): The real part and the magnitude of the complex envelope, respectively. The two spectral modes are similar except at the edges and around 2 ppm (near NAA and NAAG). On panel (a), a destructive interference yields a dip between the NAA and NAAG peaks close to 2 ppm. Near 2 ppm, a constructive interference on panel (b) enhances (nearly by a factor of 2) the peak height of NAA, whereas merely a slight shoulder is seen for NAAG (Color figure online)

Full size image

Further comparison is provided in Fig. 2 between these two modes of the non-derivative \((m=0)\) non-parametric envelopes from the \({\mathrm{FPT}}^{(+)}\). As reference, the real part \({\mathrm{Re}}(P^+_K/Q^+_K)\) and magnitude \(|P^+_K/Q^+_K|\) borrowed from Fig. 1 are re-displayed on panels (a) and (b) of Fig. 2. Herein, to facilitate another comparison, we used color coding for the two curves: blue for \({\mathrm{Re}}(P^+_K/Q^+_K)\) and red for \(|P^+_K/Q^+_K|\). On panel (c) these two curves are superimposed. Thereby, the similarities and differences, as noted in Fig. 1, are even clearer. In particular, the extent by which a large part of the magnitude envelope is lifted relative to the real part of the complex spectrum becomes more evident than in Fig. 1. This is the case especially in the region between 1.5 and 2.0 ppm, as well as from about 2.7 to 3.7 ppm. Moreover, a marked baseline lifting also occurs above 4.2 ppm due to constructive interference of metabolite resonances with the remnants from the still strong water residual tail. Also, the peak at about 2.0 ppm is substantially widened in the magnitude mode, as is the peak centered at about 3.5 ppm corresponding to a part of the myo-inositol (m-Ins) triplet. Also noteworthy is that the valley in between creatine, Cr, and phosphocreatine (PCr) at 4.0 ppm is much deeper in the real part of the complex envelope compared to the magnitude mode.

Fig. 2
figure 2

Two spectral modes (real part vs. magnitude) in the \({\mathrm{FPT}}^{(+)}\) at the signal length \(N=1024\, (K=512)\). Panels (a, b) borrowed from Fig. 1: the real part and the magnitude of the complex envelope, respectively. Panel (c): Superimposed envelopes from the real (blue) and magnitude (red) modes to expose more clearly the extent of broadening of the peak widths in the latter versus the former mode (Color figure online)

Full size image

Next we pass to derivative signal processing. With this goal, we begin with the Fourier technique in Fig. 3. Panels (a) and (b) of this figure display the magnitude envelopes in the FFT and dFFT, respectively. In the case of the dFFT, the 4th order derivative envelope \(|(\text {d}/\text {d}\nu )^4F_k|\) is depicted on panel (b). In the dFFT (panel b: \(m=4)\), the resonance widths are generally broadened compared to the FFT (panel a: \(m=0)\). Panel (c) superimposes the two envelopes, one from the FFT \((m=0)\) and the other from the dFFT \((m=4)\). Therein, much of the finer structure is lost with the derivative operator \((\text {d}/\text {d}\nu )^4\) in the dFFT. Relative to the FFT \((m=0)\), not only that most of the peaks are broadened, but they are also much coarser in the dFFT \((m=4)\).

Fig. 3
figure 3

Comparison of the magnitude modes of the complex envelopes in the FFT (a: \(m=0\)) and normalized dFFT (b: \(m=4\)) at the signal length \(N=1024\) (\(K = 512\)). Panel (c): Superimposed envelopes from the FFT (green) and dFFT (red). Therein, resonance widths are broadened in the dFFT \((m=4)\) compared to the FFT \((m=0)\). This tendency worsens with augmentation of the derivative order in the dFFT (not shown) (Color figure online)

Full size image

Overall, Fig. 3 amply illustrates the effect of the apodizing power function \((-2\pi in\tau )^4\) in the dFFT. As per the Theory section, this apodization is generated by applying the derivative operator \((\text {d}/\text {d}\nu )^4\) to the complex envelope in the FFT. At larger times \(n\tau \), encoding is corrupted by noise more than at the earlier instances of data acquisition. Noise is further amplified by the power function \((n\tau )^4\). Such a circumstance exacerbates the situation with the dFFT \((m=4)\), which processes the noise amplified FID data set, \(\{(n\tau )^4c_n\}\). This broadens the peaks from the FFT \((m=0)\), which processes directly the FID data points, \(\{c_n\}\).

Figure 4 deals with non-derivative and derivative magnitude envelopes from the \({\mathrm{FPT}}^{(+)}\) and \({\mathrm{dFPT}}^{(+)}\), respectively. Therein, the non-derivative \((m=0)\) envelope \(|P^+_K/Q^+_K|\) (panel a) from the \({\mathrm{FPT}}^{(+)}\) is compared with the 4th derivative envelope \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\) (panel b) from the \({\mathrm{dFPT}}^{(+)}\). The former is marked in green, whereas the latter is in blue. Strikingly, in panel (b) for \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\), a substantial portion of the spectrum is significantly lowered almost to the baseline. Most notably between 2.5 and 3.3 ppm, a large number of distinct, narrow resonances have emerged. Also on panel (b), centered at 2.0 ppm, there are now two clearly distinguished peaks with a deep bifurcation between them.

Fig. 4
figure 4

Comparison of the magnitude modes of the complex envelopes in the \({\mathrm{FPT}}^{(+)}\) (a: \(m=0\)) and \({\mathrm{dFPT}}^{(+)}\) (b: \(m=4\)) at the signal length \(N=1024 (K=512)\). Panel (c): Superimposed envelopes for \(m=0\) (red) and \(m=4\) (blue). Herein, relative to the \({\mathrm{FPT}}^{(+)}\)\((m=0)\), the \({\mathrm{dFPT}}^{(+)}\)\((m=4)\) shows resolution improvement (splitting overlapped peaks) and enhanced SNR (lowering the background baseline) (Color figure online)

Full size image

In Fig. 4, the two envelopes from panels (a) and (b) are superimposed in panel (c). On this plot, numerous thin resonances and the substantial lowering of the baseline in \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\) are clearly exhibited in direct comparison with the non-derivative magnitude mode \(|P^+_K/Q^+_K|\). It is fascinating that the relatively low 4th order derivative in the \({\mathrm{dFPT}}^{(+)}\) is capable of converting many complexes from \(|P^+_K/Q^+_K|\) into clear and well-delineated resonances, a good number of which are nearly isolated despite being tightly packed. This shows that even a low-order derivative envelope in the \({\mathrm{dFPT}}^{(+)}\) can split apart spectrally crowded regions that are largely smeared and/or smoothed out in \(|P^+_K/Q^+_K|\) from the non-derivative \({\mathrm{FPT}}^{(+)}\). Such an observation is evidenced in Fig. 4c throughout the displayed range of frequencies.

Figure 5 is concerned with comparison between the two derivative shape estimators (Padé and Fourier): the \({\mathrm{dFPT}}^{(+)}\) and dFFT. Panels (a) and (b) show the 4th derivative envelopes \(|(\text {d}/\text {d}\nu )^4F_k|\) (Fourier) and \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\) (Padé), respectively. In Fig. 5 too, as in the first four figures, the full signal length \(N=1024\) is employed by both processors, dFFT and \({\mathrm{dFPT}}^{(+)}\). As stated, this includes only one zero filling (original 512 data points appended by 512 zeros). In Fig. 5a for \(|(\text {d}/\text {d}\nu )^4F_k|\) in the dFFT, the derivative Fourier envelope appears to be substantially rougher. It shows a notably lower resolution than even the non-derivative Padé envelope \(|P^+_K/Q^+_K|\) in the \({\mathrm{FPT}}^{(+)}\) from Fig. 4a.

Specifically, all of the peaks on panel (a) in Fig. 5 are wide and rough, with minimal or no indication of any substructures. There is practically no hint of a split in the broad peak at 2.0 ppm, where merely a vague shoulder is noticed. Moreover, none of the abundant fine resonances between 2.5 and 3.3 ppm can be seen. At about 3.7 ppm, a single wide structure is present, instead of the triple serrated peak corresponding to m-Ins.

Panel (b) of Fig. 5 re-depicts the Padé magnitude envelope of the 4th derivative \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\) borrowed from Fig. 4b for comparison with the dFFT. On panel (c), the curves from panels (a) and (b) are superimposed to have a clearer differentiation between the \({\mathrm{dFPT}}^{(+)}\) and dFFT. Therein, almost everywhere, the Fourier envelope (red), with its multiple bumps, appears as some wider umbrellas over sharply separated individual resonances from the Padé envelope (blue). Often, e.g. a single spectral umbrella from the dFFT is seen to cover 2–4 resonances from the \({\mathrm{dFPT}}^{(+)}\). As such, the entire panel (c) of Fig. 5 looks as if the blue Padé envelope were the components of the red Fourier envelope. This analogy helps trace the reason for which the Fourier envelope is nearly everywhere more elevated than the Padé envelope. For example, by not being able to resolve NAA and NAAG through destructive interference, the dFFT performs constructive interference near 2 ppm. This yields a wide shoulder at the location of the dip/valley between NAA and NAAG. So indeed, close to 2 ppm, the Fourier red curve is like a folding (convolving) function of the Padé blue curve. Similar analogies also apply to all the other parts of the spectrum where Fourier complexes are juxtaposed to Padé splitting of compound spectral structures. In other words, Fig. 5c gives the impression that the red Fourier curve is the total shape spectrum of the blue Padé component-type spectra, even though both curves are merely the non-parametric envelopes.

The just evoked analogy, “a spectral umbrella-type envelope” versus “the hidden, resolved peaks” is extended to Fig. 6 by posing the question: could the blue substructures underneath their red cover be revealed by the dFFT itself through drawing the Fourier curve at more frequency points than those used in Fig. 5? A denser sampling grid on the chemical shift axis in the dFFT can easily be achieved by putting more zeros at the tail of the FID, i.e. by using repeated zero-padding of the time domain data. To check for this eventuality, we systematically enhanced zero-filling of the originally encoded 512 FID data points (from 2 to 8 times). The sole outcome of this exercise with 2–8 zero-padding is to minimally smooth the appearance of the derivative Fourier spectrum. For instance, with 8 zero-fillings, the only improvement by the dFFT in Fig. 6 is in smoothing the rough edges of spectral “bumps”. These bumps are still seen as the same type of umbrellas over the well resolved envelope in the \({\mathrm{dFPT}}^{(+)}\). However, the substantive problem inquiring whether any new splitting by the dFFT in Fig. 6 occurs is not solved. Namely, there are no further resonances identified by the dFFT in Fig. 6 with 8 zero-fillings compared to Fig. 5 with only one zero filling. This implies that supplementary zero-padding of the FID does not improve resolution in the dFFT.

Fig. 5
figure 5

Comparison of the magnitude modes of the complex envelopes in the dFFT (a) and \({\mathrm{dFPT}}^{(+)}\) (b) for \(m=4\) at the signal length \(N=1024 (K = 512)\). Panel (c): Superimposed envelopes from the dFFT (red) and \({\mathrm{dFPT}}^{(+)}\) (blue). Superior performance of the \({\mathrm{dFPT}}^{(+)}\) over that in the dFFT is seen. A number of peaks are clearly delineated and resolved in the \({\mathrm{dFPT}}^{(+)}\) as opposed to the dFFT (Color figure online)

Full size image
Fig. 6
figure 6

Comparison of the magnitude modes of the complex envelopes in the dFFT (a: 8 zero-fillings, \(N=4096\)) and \({\mathrm{dFPT}}^{(+)}\) (b: one zero-filling, \(N=1024\)) for \(m=4\). Panel (c): Superimposed envelopes from the dFFT (red) and \({\mathrm{dFPT}}^{(+)}\) (blue). Extensive zero-filling of the encoded FID does not improve resolution in the dFFT. Only smoothing is obtained in the dFFT for extra zero-padding, but the Fourier compound structures remain unresolved (Color figure online)

Full size image

Figures 7 and 8 are of the type of Figs. 5 and 6 with the same idea of comparison between \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\) (Padé) and \(|(\text {d}/\text {d}\nu )^4F_k|\) (Fourier). The only difference is that Figs. 7 and 8 single out two specific narrower frequency intervals. They both zoom into two spectral subregions of particular interest, one from 1.8 to 2.5 ppm (right column, panels e–h), and the other from 2.75 to 3.5 ppm (left column, panels a–d). These are the two innermost regions of the initial wider ROI extending from 0.75 to 4.2 ppm. Figures 7 and 8 begin with the non-derivative reference envelopes in the \({\mathrm{FPT}}^{(+)}\) (Padé) for the two frequency subintervals, \(\nu \in [2.75,3.5]\) ppm (panel a) and \(\nu \in [1.8,2.5]\) ppm (panel e).

Fig. 7
figure 7

Comparison of the magnitude modes of the complex envelopes in the dFFT and \({\mathrm{dFPT}}^{(+)}\) at the same signal length \(N=1024\). The left and right columns are for two innermost bands of the fuller ROI from Figs. 3, 4, 5 and 6. Panels (a, e): The \({\mathrm{FPT}}^{(+)}\)\((m=0)\), (b, f): the dFFT \((m=4)\), (c, g): the \({\mathrm{dFPT}}^{(+)}\)\((m=4)\). Panels (d, h): Superimposed envelopes from the dFFT (red: \(m=4\)) and \({\mathrm{dFPT}}^{(+)}\) (blue: \(m=4\)). Therein, the \({\mathrm{dFPT}}^{(+)}\) narrows the peak widths, as opposed to line broadening in the dFFT (Color figure online)

Full size image

In the 1st row of Fig. 7, on panel (a) for \(|P^+_K/Q^+_K|\) in the \({\mathrm{FPT}}^{(+)}\), the most prominent resonances are identified as Cr at 3.0 ppm and Cho at 3.2 ppm. Further, N-acetyl aspartate, NAA centered at 2.0 ppm is the largest resonance on panel (e) of Fig. 7, which is also for \(|P^+_K/Q^+_K|\). Several other peaks are suggested within the highly elevated rolling baseline. Long-extending tails of each individual resonance contained in the magnitude envelope add up and this lifts the entire spectrum considerably above the abscissae on panels (a) and (e) of Fig. 7.

The 2nd row (panels b and f) of Fig. 7 presents the envelopes \(|(\text {d}/\text {d}\nu )^4F_k|\) (Fourier) generated by the dFFT. Comparing the envelopes on panels (a) and (b), the latter spectrum is overall much rougher with only a suggestion of an extended and bumpy Cho peak around 3.2 ppm. On panel (f) for the dFFT, it is seen that the NAA peak is markedly broadened. Also, most of the other structures from the reference envelope on panel (e) are obscured on panel (f).

Panels (c) and (g) in the 3rd row of Fig. 7 show the \({\mathrm{dFPT}}^{(+)}\) as the 4th derivative envelopes \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\). Therein, notably, on panel (c), some sixteen clear peaks have emerged, and they are resolved to such an extent that some of the valleys among them descend almost or completely to the background baseline. For example, in the narrow region between 3.13 and 3.3 ppm, five distinct resonances are predicted by the \({\mathrm{dFPT}}^{(+)}\). Of utmost importance is that the choline compound in the middle of panel (c) is fully resolved into all its components even though shape estimation alone is used in the \({\mathrm{dFPT}}^{(+)}\). Similarly, between 2.75 and 2.95 ppm, another five peaks are clearly delineated. In panel (g) there is a marked split between NAA and NAAG at 2.0 ppm, and several other fine structures have begun to emerge, as well.

In the 4th row of Fig. 7, on the bottom panels (d) and (h), shown together are the two envelopes: \(|(\text {d}/\text {d}\nu )^4F_k|\) (Fourier, red) and \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\) (Padé, blue). Herein, it becomes patently clear what we meant earlier on Figs. 5 and 6 by the notion of “the Fourier envelope-umbrella” extending over the Padé envelope. Some eight peaks from the Padé envelope \(|(\text {d}/\text {d}\nu )^4P^+_K/Q^+_K|\) are seen underneath 2–3 bumps in the Fourier envelope \(|(\text {d}/\text {d}\nu )^4F_k|\) even in the two very narrow frequency bands \(\nu \in [3.05,3.30]\) (panel d) and \(\nu \in [2.10,2.35]\) ppm (panel h).

Fig. 8
figure 8

Comparison of the magnitude modes of the complex envelopes in the dFFT (8 zero-fillings, \(N=4096\)) and \({\mathrm{dFPT}}^{(+)}\) (one zero-filling, \(N=1024\)) for \(m=4\). The left and right columns are for two innermost bands of the fuller ROI from Figs. 3, 4, 5 and 6. Panels (a, e): The \({\mathrm{FPT}}^{(+)}\)\((m=0)\), (b, f): the dFFT \((m=4)\), (c, g): the \({\mathrm{dFPT}}^{(+)}\)\((m=4)\). Panels (d, h): Superimposed envelopes from the dFFT (red: \(m=4\)) and \({\mathrm{dFPT}}^{(+)}\) (blue: \(m=4\)). Herein, no resolution improvement in the dFFT follows from more extensive zero-filling of the encoded FID. Merely smoothing without resolving the Fourier amalgamated structures is observed in the dFFT (Color figure online)

Full size image
Fig. 9
figure 9

Comparison of the magnitude modes of varying derivative orders \((m=1-4)\) in the \({\mathrm{dFPT}}^{(+)}\). The left and right columns are for the two outermost bands of the fuller ROI from Figs. 3, 4, 5 and 6. These bands are complementary to those from Figs. 7 and 8. Panels (a, f): The \({\mathrm{FPT}}^{(+)}\)\((m=0)\). Panels (b–j) are all for the \({\mathrm{dFPT}}^{(+)}\) as: \(m=1\) (b, g), \(m=2\) (c, h), \(m=3\) (d, i) and \(m=4\) (e, j). Relative to \(m=0\) (a), better resolution is seen for \(m=1-4\) (b–e) (Color figure online)

Full size image

Figure 8 differs from Fig. 7, only in the number of zero fillings of the FID for the dFFT, eight and one, respectively. Seeing the sharp edges on panels (d) and (h) of Fig. 7, at the Fourier grid frequencies in the dFFT with the signal length \(N=1024\) (one zero filling), it is tempting to think that some splitting of the Fourier envelope-umbrella could potentially be achieved by further zero-filling of the FID. This would lead to the additional Fourier grid frequencies (exceeding the number 1024) on the chemical shift axis. Such hopes proved futile on Fig. 6c for the wider ROI within \(\nu \in [1.1,4.2]\) ppm and definitely remain so in the two narrower bands \(\nu \in [2.75, 3.50]\) (panel d) as well as for \(\nu \in [1.8, 2.5]\) ppm (panel h) from Fig. 8. Zooming in Fig. 8 relative to Fig. 6 only puts in more transparent evidence the fact that extensive zero-padding for estimation by the dFFT may smooth out the sharp edges, but cannot improve resolution.

In Fig. 9, we focus exclusively upon the performance of the \({\mathrm{dFPT}}^{(+)}\) for \(|(\text {d}/\text {d}\nu )^mP^+_K/Q^+_K|\) with four consecutive derivative orders \((m=1-4)\). As reference spectral data, the non-derivative \((m=0)\) variant \(|P^+_K/Q^+_K|\) in the \({\mathrm{FPT}}^{(+)}\) is given on panels (a) and (f). To complement Figs. 7 and 8, here in Fig. 9, we consider the two outermost frequency intervals, one from 3.4 to 4.0 ppm (left column, panel a–e) and the other from 1.0 to 1.8 ppm (right column, panels f–j). The first range \(\nu \in [3.4,4.0]\) ppm on the left column of Fig. 9 contains two groups of different molecules, one with Glx in a narrow band (3.82–3.91 ppm), and the other with m-Ins, which is spread out in a wider band (3.50–3.82 ppm). Here, Glx is the common abbreviation for Glutamine (Gln) and Glutamate (Glu). The right column in the interval \(\nu \in [1.0,1.8]\) ppm is a bit more diverse as it contains several molecules, such as lipid (Lip), lactate (Lac), alanine (Ala) and acetate (Ace).

On the left column of Fig. 9, relative to the case with \(m=0\) (panel a), the main m-Ins resonance near 3.77 ppm remains unresolved and appears to be wider on panels (b) and (c) for the 1st \((m=1)\) and second \((m=2)\) derivative, respectively. As to Glx, these two lowest derivative transforms are seen to achieve some initial splitting, albeit rough. Similarly, width broadening for resonances assigned to Lip, Lac, Ala and Ace also occurs with the 1st and 2nd derivatives on the respective panels (g) and (h) from the right column of Fig. 9.

On the right column of Fig. 9, relative to the 2nd and 3rd derivative envelopes on panel (h) and (i), respectively, the 4th derivative envelope on panel (j) brings a clear progress toward better resolution. For example, the trend of splitting in Ala as well as in one of the Lac doublets (upstream, near 1.3 ppm) begins to appear on panel (j). A closer look shows that Lip and Lac are, in fact, less resolved on panel (h) for \(m=2\) than on panel (i) for \(m=3\). Eventually, better resolution within these latter two resonance complexes is achieved for \(m=4\), as seen on panel (j). Recall, by reference to study [15], some oscillatory patterns in derivative lineshape profiles at smaller values of m are not unexpected. Namely, the same occurrence has also been detected with synthesized time signals (both noise-free and noise-contaminated) for lower-order derivative envelopes from the \({\mathrm{dFPT}}^{(-)}\) [15,16,17].

However, in closing the analysis of Fig. 9 by returning to its left column, it can be seen that the 3rd derivative envelope \((m=3)\) on panel (d) is actually better resolved than that for the 2nd derivative envelope the \((m=2)\) on panel (c). On panel (d), the triplet in Glx begins to emerge. Moreover, the triplet in m-Ins, which appeared as a single symmetrical Lorentzian in the 1st derivative envelope \((m=1)\) on panel (b), now acquires two slight shoulders on panel (d) for \(m=3\) around the central peak (\(\sim \) 3.77 ppm). Finally, on panel (e) for the 4th derivative envelope \((m=4)\), the shoulders in m-Ins clearly indicate that there are indeed two side peaks surrounding the central resonance at 3.77 ppm. Also, the two other m-Ins triplets are nicely resoloved in the frequency interval \(\nu \in \) [3.5, 3.72] ppm. Moreover, a deeper splitting within Glx permits the emergence of the constituent, distinctly delineated peaks.

4.2 Clinical relevance of the studied MRS problem

For early tumor diagnostics, molecular imaging is vital. This is the case because malignant transformation begins to develop at the molecular level. Diagnosing such pathologic alterations at their early stages can greatly improve the chances for more successful tumor control, potential cure, and impact positively on the patient’s survival and quality of post-therapeutic life.

We have herein applied the dFPT to process MRS time signals encoded in vivo on a 1.5 T MR clinical scanner from a pediatric patient with cerebral asphyxia [4, 6, 18]. This problem was selected, in part, because MR spectra associated with ischemia or asphyxia are very dense and, thus, have represented a great challenge for resolution and interpretation. This clinical problem is, furthermore, of major importance within neurodiagnostics for children, as well as for the adult population. Specifically, a critical differential diagnostic dilemma in neuro-oncology is to distinguish non-tumorous cerebral hypoxia/ischemia from brain tumors [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39].

Heretofore, in vivo MRS within neurodiagnostics has mainly relied upon a few metabolites and their concentration ratios. Among these are N-acetyl aspartate, NAA, resonating at about 2.0 ppm, which reflects viability and abundance of neurons. The concentration of NAA is generally low with cerebral ischemia or asphyxia [21], since cerebral neurons are exceedingly vulnerable to hypoxia. However, in brain tumors, as well as in most other brain injury, diminished NAA is also generally reported [1, 3]. A lactate, Lac, doublet, resonating at approximately 1.3 ppm, is often observed in MR spectra with cerebral ischemia/hypoxia, due to the predominance of anaerobic glycolysis. In brain tumors, Lac is frequently seen, as well [21,22,23]. In approximately the same chemical shift region, i.e. at about 1.3 ppm, the appearance of lipid, Lip, has also been reported, both in association with reperfusion after hypoxia and in brain tumors.

There are additional challenges for neurodiagnostics through MRS in the pediatric population, since account must be taken of the age of the child. In newborns, myo-inositol, m-Ins, is the dominant brain metabolite. Total choline, Cho, resonating at about 3.2 ppm is reportedly the largest peak in older infants [24]. However, Cho as an indicator of membrane damage, cellular proliferation and cell density is also very often reported to be elevated in brain tumors, as well as in many non-tumorous brain pathologies [1, 3]. Cerebral energy metabolism is reflected by creatine, Cr, which resonates at about 3.0 ppm and whose concentration in the brain is usually stable after the first year of life [23]. Concentrations of NAA increase as the child’s brain matures. Thus, concentration ratios of Cho to NAA and Cho to Cr reportedly fall with the normal development of the child’s brain [24].

One rationale for the clinical focus of this work is to improve neurodiagnostics through MRS, striving toward greater accuracy in identifying cerebral hypoxia/ischemia versus brain tumors, as well as other pathology that present differential diagnostic dilemmas. This aim is highly relevant for pediatric patients, in addition to addressing the pressing needs of the adult population.

These investigations can have wider implications for oncology. Notably, hypoxic areas are frequently found not only within brain tumors, but also within many other malignancies. The hypoxic regions are generally resistant to treatment, both by radiation or chemotherapy. In addition, hypoxia generates genomic instability and promotes invasiveness [40]. Thus, identifying hypoxic regions via MRS could help improve cancer treatment planning. Of particular note is that phosphocholine, PC, is an identified cancer biomarker which also reflects hypoxia [41].

The present feasibility study is of high clinical relevance for in vivo MRS, and motivates further applications of the derivative fast Padé transform to encoded data associated with various tumors from e.g. brain, breast, prostate, ovary and other organs. This is what we are set to do in the nearest future.

5 Discussion and conclusion

The most recent advance in magnetic resonance spectroscopy, MRS, is the introduction of the derivative fast Padé transform, dFPT. This processor complements the standard non-derivative fast Padé transform, FPT. The dFPT was benchmarked on synthesized MRS time signals or FIDs that were both noise-free and noise-contaminated. It has been shown that the dFPT for high-derivative orders is capable of exactly reconstructing all the peak parameters for every physical resonance by carrying out solely estimation of total shape spectra or envelopes.

In the dFPT, the linewidths are narrowed and the peak heights concomitantly enhanced. Simultaneously, the dFPT suppresses noise by lowering the tails of lineshape profiles. This leads to clearly visualized separation of overlapping peaks, implying resolution improvement and noise reduction. All of the mentioned features of the dFPT have been established with a primary focus on high-order derivatives. The reason is that in such cases the derivative envelopes reduce to the underlying derivative components that, in turn, permit unequivocal quantification.

The topic of the present paper is extension of the dFPT to encoded FIDs. However, here, the focus is exclusively on low-order derivatives (from the 1st to the 4th). We applied the dFPT to MRS time signals encoded at short echo time (24 ms) with static magnetic field of strength 1.5 T. The total signal length was only 512 data points. Both the FPT (no derivatives) and dFPT (multiple derivatives) are implemented in their non-parametric forms (conventionally used in shape estimation).

The dFPT for low-order derivatives based on encoded FIDs is shown here to uphold the two main features of this processor by improving frequency resolution and signal-to-noise ratio, SNR, relative to the customary, non-derivative FPT. Low-order derivatives of envelopes are chosen because this is a feasibility study aimed at visualizing the improved resolution and SNR. For testing the precision of estimation in the present applications of the non-parametric dFPT, we plan to use the parametric dFPT, similarly to the previous studies on synthesized FIDs [15,16,17].

The magnitude mode of spectral envelopes is employed because it is phase-insensitive. This means that each isolated peak is a symmetric Lorentzian for any derivative order. Such a lineshape is most convenient for interpretation as if it were a non-derivative pure absorption mode with no phase. In the usual quantification with non-derivative estimators, the magnitude mode is not preferred because the peak width is broadened by a factor of \(\sqrt{3}\) relative to a symmetric absorptive Lorentzian. Some peak width broadening also exists in the present magnitude mode of lower-order derivatives of complex-valued amplitudes. However, this has nothing to do with the fact that we are dealing here with encoded MRS data. Similar oscillatory patterns of resonances toward their stabilization have also been reported with simulated MRS time signals processed by low-order dFPT.

There is a clear trend of resolution enhancement when going from the 1st to the 4th derivative envelope in the dFPT. In this regard, especially the 4th order dFPT is performing excellently since it brings significant improvement over the 2nd order dFPT. As expected, derivative resonances cannot fully stabilize at low-derivative orders. Nevertheless, this trend is important as it sets the stage for convergence of normalized envelopes at higher derivative orders to be investigated soon using encoded MRS time signals.

Notably, even though low-order derivatives are used, the dFPT is still found to perform significantly better than the non-derivative FPT. This fact itself justifies the use of low-order derivatives in the dFPT as the first attempt in improving the performance of non-derivative estimations.

The derivative Fourier transform, dFFT, is observed to widen the resonance widths relative to the FFT. Furthermore, with increasing derivative order, the performance of the dFFT is worsened with respect to the FFT. This is just the opposite to what is experienced with the dFPT versus FPT (presently, for encoded FIDs, and in our previous reports for synthesized MRS time signals). Noise has a much more detrimental effect on the dFFT than on the FFT. The reason is in the fact that the derivative operator in the dFFT apodizes the input FID with a power function of time raised to the given derivative order. This apodization amplifies noise since the multiplicative time power function puts the emphasis on the FID tail which is itself more corrupted with noise than the earlier encoded FID data points.

By contrast, the same derivative operator acts entirely differently in the dFPT: it processes the intact FID with no apodization whatsoever. The operation of derivation is applied to a complex spectrum in the FPT and this leads to no time power function (nor to any other apodizing function for that matter). The envelope in the FPT is given by an analytical expression (a polynomial quotient) and so is the general derivative envelope in the dFPT (a closed formula in terms of the Bell polynomials). A key consequence is that the derivative operator itself does not generate any additional noise in the dFPT. Quite the contrary, the derivative operator in the dFPT suppresses noise because it flattens the rolling background by narrowing the peak widths and simultaneously elongating the peak heights compared to the non-derivative FPT.

Overall, derivative envelopes from the dFPT have a better SNR than those in the FPT. In particular, the up-lifted rolling background from the FPT is considerably reduced in the dFPT. Resonances that were packed and superimposed on high background in the FPT are brought down by the dFPT close to the zero-valued baseline which is almost embedded in the chemical shift axis.

Moreover, derivative envelopes in the dFPT are better resolved than those in the FPT. A number of compound, amalgamated spectral structures in the FPT are disentangled by the dFPT into the constituent isolated peaks. Thus, especially with the 4th order derivative, the dFPT is shown to fully resolve the hidden multiplets of glutamine-glutamate complex, Glx, as well as those associated with myo-inositol, m-Ins and taurine, Tau. The same 4th order derivative in the dFPT completely splits apart the choline-containing compound into its several components. It is the first time that a shape estimator was able to clearly visualize all the constituents of the choline multiple structure in a total shape spectrum. This finding is all the more significant by the fact that merely a low-order derivative in the dFPT is applied to the encoded MRS time signals. The present feasibility study holds promise for further exploration within the ’theme cancer’ in neurodiagnostics, and well beyond.