1 Introduction

The present study is on the estimation problem by parametric and non-parametric analysis in magnetic resonance spectroscopy (MRS) [1, 2]. In MRS, times signals are measured, whereas spectra are computed. Time signals are not directly interpretable due to many coupled components that are attenuated sinusoids and cosinusoids. At first glance, spectral lineshapes have an advantage since they present at least some recognizable structures appearing as peaks or resonances.

However, this is only an apparent advantage because the initial obstacle with tightly packed complex damped harmonics in the input time signal did not magically disappear from the subsequent frequency-domain analysis. Rather, the related trouble exists in spectral lineshapes too where, in fact, coupled time-signal constituents are now disguised through their direct influence on the emergence of closely overlapping resonances.

The art is then in splitting apart the overlapped peaks even if the input time signals were noise-free as in the case with theoretical simulations. Measured time signals dramatically aggravate the situation due to unavoidable contamination of recorded data with noise.

Historically, it was the fast Fourier transform (FFT) [1], which helped research in data analyses across inter-disciplinary fields. Depending on total signal lengths, the use of the FFT could mean huge reduction of computational demands. Nevertheless, computational speed is not all that is needed. It may be necessary but not sufficient. For example, nowadays, the ever pursuing quest for much improved resolution de facto rules out the FFT from useful applications especially in MRS. The reason is low resolution of the FFT coupled with its other well-documented drawbacks, including e.g. linearity which translates into inability to suppress noise from the input time-signal data.

Therefore, to better reconstruct the physical information and improve resolution, it is critical to have a signal processor which can accurately analyze the measured time signals by simultaneously reducing noise and separating overlapping resonances. We show that such a key twofold goal for encoded time signals can reliably be achieved by using the fast Padé transform (FPT) [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55]. This is demonstrated by several representative illustrations of both the shape and parameter estimations with the FPT.

2 Theory and methods

2.1 Mathematics of signal processing

To account for the physical nature of time signals, or equivalently, free induction decay (FID) curves, encountered in MRS, the following decomposition is appropriate [1, 2]:

$$\begin{aligned} c_n=\sum \limits _{k=1}^Kd_k{\mathrm{e}}^{2\pi i\nu _kn\tau }\quad (0\le n\le N-1), \end{aligned}$$
(2.1)

where \(\tau \) is the sampling or dwell time, N is total signal length and i is the imaginary unity, \(i=\sqrt{-1}.\) Here, \(\nu _k\) is the fundamental complex linear frequency, \(\nu _k={\mathrm{Re}}(\nu _{\mathrm{k}})+i{\mathrm{Im}}(\nu _{\mathrm{k}}),\) while \(d_k\) is the associated complex amplitude, \(d_k=|d_k|\exp {(i\varphi _k)},\) whose magnitude and phase are \(|d_k|\) and \(\varphi _k,\) respectively. Quantity \({\mathrm{Re}}(\nu _{\mathrm{k}})\) represents the chemical shift which is a dimensionless frequency at which the \(k\,\)th molecule resonates with the external excitation. The imaginary part \({\mathrm{Im}}(\nu _{\mathrm{k}})\) of the complex frequency \(\nu _k\) is proportional to the full width at half maximum (FWHM) of the \(k\,\)th resonance. The reciprocal of \({\mathrm{Im}}(\nu _{\mathrm{k}})\) is proportional to the lifetime of the \(k\,\)th resonance.

For the estimation problem, the input data in the time domain are the FID data points, \(\{c_n\}.\) These can be either synthesized (i.e. simulated theoretically through computations), or available through experimental data from measurements (encoding). In the case of the frequency domain, the input data for estimation is a single response function in the form of the finite-rank Green function \(G_N(z^{-1})\):

$$\begin{aligned} G_N(z^{-1})=\frac{1}{N}\sum \limits _{n=0}^{N-1}c_nz^{-n},\quad z^{-1}={\mathrm{e}}^{-2\pi i\nu \tau }. \end{aligned}$$
(2.2)

This is a Maclaurin polynomial of degree \(N-1\) with the expansion coefficients given by the input signal points \(\{c_n\}\, (0\le n\le N-1).\) If the time signal were infinitely long (\(N=\infty \)), function \(G_\infty (z^{-1})=\lim _{N\rightarrow \infty }G_N(z^{-1})\) would become a Maclaurin series (with an unlimited number of terms) in which case it is necessary to define the convergence radius. Since such a series is an expansion in powers of \(z^{-1},\) its convergence radius would be outside the unit circle \((|z|>1).\)

2.2 Discrete and fast Fourier transforms

The discrete Fourier transform (DFT) is deduced directly from Eq. (2.2) by taking the sweep frequency \(\nu \) to be at the Fourier grid frequencies \(k/T\, (k\le 0\le N-1),\) where T is the total duration of the FID (i.e. the total acquisition time in encoding), \(T=N\tau \):

$$\begin{aligned} {\mathrm{DFT}}_k=\frac{1}{N}\sum \limits _{n=0}^{N-1}c_n{\mathrm{e}}^{-2\pi ink/N}\quad (0\le k\le N-1). \end{aligned}$$
(2.3)

In other words, to pass from (2.2) to (2.3) all that is needed is to discretize (digitize) the continuous running frequency \(\nu \) according to \(\nu =k/T\, (k=0,1,2,\ldots , N-1).\) Afterward, using the relation \(T=N\tau ,\) it follows that \(\nu \tau =k/N,\) which is the term appearing in the exponential from (2.3). Further, the FFT is defined by means of the DFT from Eq. (2.3) by employing therein the Tukey-Cooley fast algorithm. This algorithm converts the initial \(N^2\) multiplications into \(N{\mathrm{log}}_2N\) multiplications, when N is a composite number of the form \(N=2^\ell \, (\ell =1,2,3,\ldots )\). The starting \(N^2\) multiplications in the DFT are due to multiplications of two vectors of length N: one is the so-called work vector \(W_k\) where \(W_k={\mathrm{e}}^{-2\pi ink/N}\, (0\le k\le N-1)\) and the other is the vector comprised of N time signal points \(\{c_n\}\, (0\le n\le N-1).\) For smaller N,  the DFT can occasionally be even faster than FFT. However, for large N (for which the so-named fast algorithms are used), the FFT introduces an enormous advantage over the DFT.

2.3 Fast Padé transform

The spectra defined by e.g. the diagonal FPT are defined by the two equivalent and complementary forms denoted by \({\mathrm{FPT}}^{(-)}\) and \({\mathrm{FPT}}^{(+)}\) depending whether they use the independent variable \(z^{-1}\) or z,  respectively, to connect to the input data \(G_N(z^{-1}).\) These spectra in the \({\mathrm{FPT}}^{(\pm )}\) are Padé-Green functions denoted by \(G^\pm _K(z^{\pm 1}).\) Their diagonal forms are non-linear response functions given by the polynomial ratios of a common degree K with the complex harmonic as the independent variable, which can be either \(z^{-1}\) or z. Both forms \(G^-_K(z^{-1})\) and \(G^+_K(z)\) of the polynomial quotients \(P^-_K(z^{-1})/Q^-_K(z^{-1})\) and \(P^+_K(z)/Q^+_K(z)\) are unique in their respective independent variables \(z^{-1}\) and z:

$$\begin{aligned} G^-_K(z^{-1})= & {} \frac{P^-_K(z^{-1})}{Q^-_K(z^{-1})},\quad |z|>1, \end{aligned}$$
(2.4)
$$\begin{aligned} G^+_K(z)= & {} \frac{P^+_K(z)}{Q^+_K(z)},\quad |z|<1. \end{aligned}$$
(2.5)

The polynomials \(P^\pm _K(z^{\pm 1})\) and \(Q^\pm _K(z^{\pm 1})\) can be represented as the finite developments:

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

where \(\{p^\pm _r,q^\pm _s\}\) are the expansion coefficients that are yet to be determined by a straightforward procedure which will be given in the continuation of this exposition. In the meantime, notice that there is no free term \(z^0\) in polynomial \(P^+_K(z),\) or stated equivalently, \(p^+_0\equiv 0.\)

By default, functions \(G^+_K(z)\) and \(G^-_K(z^{-1})\) converge inside \((|z|<1)\) and outside \((|z|>1)\) the unit circle, respectively. However, by the Cauchy analytical continuation, both \(G^+_K(z)\) and \(G^-_K(z^{-1})\) converge in the entire plane of the complex frequency, i.e. inside \((|z|<1)\) and outside \((|z|>1)\) the unit circle, with the exception of the K poles that are the solutions of equations \(Q^\pm _K(z^{\pm 1})=0\).

At first glance, since \(G^+_K(z)\) and \(G^-_K(z^{-1})\) are comprised of polynomials (with the finite number K of the expansion terms), one may think that consideration of the convergence region would be unnecessary. However, implicitly both \(G^-_K(z^{-1})\) and \(G^+_K(z)\) contain their series since the polynomial reciprocals \(1/Q^-_K(z^{-1})\) and \(1/Q^+_K(z)\) expanded in powers of \(z^{-1}\) and z,  respectively, have infinitely many terms (i.e. they are series).

Thus for the reason of principle, i.e. irrespective of whether or not the inversions of \(Q^-_K(z^{-1})\) and \(Q^+_K(z)\) are carried out in some concrete computations, the issue of the convergence region has to be settled (note that, generally, the model order K need not be finite). Nevertheless, it turns out that this matter of convergence is of a purely formal nature for the functions \(G^-_K(z^{-1})\) and \(G^+_K(z)\) from Eqs. (2.4) and (2.5) because of a twofold reason: the polynomial quotient forms of these functions and the Cauchy analytical continuation.

The same input Green function \(G_N(z^{-1})\) from Eq. (2.2), as the expansion in powers of \(z^{-1},\) is approximated by two different Padé functions \(G^-_K(z^{-1})\) and \(G^+_K(z).\) The relationships of the spectra \(G^\pm _K(z^{\pm 1})\) in the \({\mathrm{FPT}}^\pm \) with \(G_N(z^{-1})\) are established through the standard definitions of the Padé approximants:

$$\begin{aligned}&G_N(z^{-1})\approx G^-_K(z^{-1})\quad {\mathrm{or}}\quad G_N(z^{-1})\approx \frac{P^-_K(z^{-1})}{Q^-_K(z^{-1})}, \end{aligned}$$
(2.8)
$$\begin{aligned}&G_N(z^{-1})\approx G^+_K(z)\quad {\mathrm{or}}\quad G_N(z^{-1})\approx \frac{P^+_K(z)}{Q^+_K(z)}. \end{aligned}$$
(2.9)

2.4 Error estimates

All meaningful mathematical models should have some reasonable ways to estimate the invoked errors when describing the given function by the employed approximations. To learn about the error estimate within the FPT, let us rewrite, e.g. (2.8) in its equivalent form:

$$\begin{aligned} G_N(z^{-1})-\frac{P^-_K(z^{-1})}{Q^-_K(z^{-1})}={\mathcal {O}}(z^{-2K-1}). \end{aligned}$$
(2.10)

Here, the remainder \({\mathcal {O}}(z^{-2K-1})\) represents a series beginning with the power \(z^{-2K-1}\):

$$\begin{aligned} {\mathcal {O}}(z^{-2K-1})\equiv \sum \limits _{n=2K+1}^\infty a_nz^{-n}, \end{aligned}$$
(2.11)

where all the expansion coefficients \(\{a_n\}\, (2K+1\le n\le \infty )\) are obtainable in terms of the input data \(\{c_n\}.\) Thus, we see that the definition (2.8) of the \({\mathrm{FPT}}^{(-)}\) provides the explicit formula for the error committed when the original Green polynomial \(G_N(z^{-1})\) is described by the Padé-Green rational polynomial \(P^-_K(z^{-1})/Q^-_K(z^{-1}).\) This can be summarized in a manner familiar to estimation problems (Input-Model\(=\)Error):

$$\begin{aligned} \left. \begin{array}{c} {\displaystyle {G_N(z^{-1})-\frac{P^-_K(z^{-1})}{Q^-_K(z^{-1})}={\mathcal {O}}(z^{-2K-1})}}\\ {\displaystyle {{\mathrm{Input}}\, -\, {\mathrm{Model}}\, = \,{\mathrm{Error}}}} \end{array}\right\} . \end{aligned}$$
(2.12)

In other words, by definition, the \({\mathrm{FPT}}^{(-)}\) for the initial function \(G_N(z^{-1})\) from (2.2), as given by (2.10), signifies that the following strict equality (in lieu of the approximate equality in (2.8)):

$$\begin{aligned} G_N(z^{-1})=\frac{P^-_K(z^{-1})}{Q^-_K(z^{-1})}, \end{aligned}$$
(2.13)

is valid through all orders \(z^{-2K}.\) This means that every coefficient \(a_n\) of the higher terms \(z^{-2K-\ell }\) with \(\ell \ge 1\) is automatically taken as zero. A similar reasoning can be extended to cover also the \({\mathrm{FPT}}^{(+)}\) with an analogous interpretation of relationship (2.9) leading to the equality:

$$\begin{aligned} G_N(z^{-1})=\frac{P^+_K(z)}{Q^+_K(z)}. \end{aligned}$$
(2.14)

2.5 Extrapolation and interpolation

Crucially, in e.g. the \({\mathrm{FPT}}^{(-)},\) the definition (2.13) itself can be used to generate all the K expansion coefficients \(\{p^-_r\}\) and \(\{q^-_s\}\) of polynomial \(P^-_K(z^{- 1})\) and \(Q^-_K(z^{- 1}),\) respectively. This is achieved through multiplication of both sides of (2.13) by \(Q^-_K(z^{-1}).\) Then, in the ensuing relation, \(G_N(z^{-1})Q^-_K(z^{-1})=P^-_K(z^{-1}),\) the product \(G_N(z^{-1})Q^-_K(z^{-1})\) is treated as a convolution.

Such a procedure yields a single system of linear equations from which the whole set of the expansion coefficients \(\{q^-_s\} \, (0\le s\le K)\) of the denominator polynomial \(Q^-_K(z^{-1})\) is obtained. The other system of linear equations for the expansion coefficients \(\{p^-_r\} \, (0\le r\le K)\) of the numerator polynomial \(P^-_K(z^{-1})\) is, in fact, an analytical expression given by the convolution of \(\{c_n\}\) with the already obtained set \(\{q^-_s\}.\) This completes the generation of the polynomial expansion coefficients \(\{p^-_r,q^-_s\}\) (and similarly for \(\{p^+_r,q^+_s\}\) starting from Eq. (2.14)).

It is a remarkably robust procedure since solving a system of linear equation is a success story of linear algebra in all of computing. It belongs to a prime category of the most accurately solved numerical problems in an exemplary stable manner [1]. The obtained solutions for \(\{q^{\pm }_s\}\) can be refined by the well-known singular value decomposition (SVD) [55], as regularly done in all our applications of the Padé-based signal processing.

Once the polynomials \(P^-_K(z^{-1})\) and \(Q^-_K(z^{-1})\) become available in the explained way, we can expand their ratio \(P^-_K(z^{-1})/Q^-_K(z^{-1})\) in a Maclaurin series in powers of \(z^{-1}\) and compare it term-by-term with the input development \(G_N(z^{-1}).\) The outcome would be pleasing: the first M terms in the said Maclaurin series for \(P^-_K(z^{-1})/Q^-_K(z^{-1})\) coincide with the 2M terms of the development in the input Maclaurin polynomial \(G_N(z^{-1}).\) This means that, if the input \(G_N(z^{-1})\) is truncated by retaining e.g. its first N / 2 terms alone, then modeling \(G_{N/2}(z^{-1})\) by \(P^-_K(z^{-1})/Q^-_K(z^{-1})\) in the \({\mathrm{FPT}}^{(-)}\) would reconstruct exactly the remaining N / 2 missing terms left out after passing from \(G_N(z^{-1})\) to \(G_{N/2}(z^{-1})\).

In other words, the \({\mathrm{FPT}}^{(-)}\) extrapolates the truncated input data set \(\{c_n\}\, (0\le n\le M-1)\) to \(\{c_n\}\, (0\le n\le 2M-1)\) for any positive integer \(M<N.\) It is in this way that the \({\mathrm{FPT}}^{(-)}\) accelerates convergence of the input \(G_N(z^{-1})\) when the signal length N is systematically augmented. In practice, such a gain translates into the ability of the \({\mathrm{FPT}}^{(-)}\) to predict the time signal points that would have been encoded had the measurements continued beyond the total acquisition time \(T=N\tau .\) The FFT cannot extrapolate as it is limited to only those time signal points that are present in the input FID of length N.

Moreover, the \({\mathrm{FPT}}^{(-)}\) can perform interpolation, as well. This is achieved by computing the given spectrum \(P^-_K(z^{-1})/Q^-_K(z^{-1})\) at any desired set of values of the sweep frequencies \(\nu .\) On the other hand, the Fourier grid frequencies are pre-assigned by the total acquisition time T. This implies that resolution of the FFT is the same for all signals with the same T. Stated equivalently, the FFT does not account for the different nature of different time signals. The number of the grid frequencies is the same as the number N of the input time signal points. Often, zero-padding is used by doubling the input FID data set through appending N zeros to the tail of the original time signal. This may somewhat ameliorate the appearance of the resulting Fourier spectrum, but resolution cannot be improved by zero-filling.

The same advantageous interpolation feature is also shared by the \({\mathrm{FPT}}^{(+)}.\) Although they are equivalent and complementary, there is a clear difference between the \({\mathrm{FPT}}^{(-)}\) and the \({\mathrm{FPT}}^{(+)}.\) The latter transform is in terms of z which means that the initial convergence radius of the \({\mathrm{FPT}}^{(+)}\) is inside the unit circle \((|z|<1)\) where the input \(G_N(z^{-1}),\) which is in terms of \(z^{-1},\) diverges. Hence, inside the unit circle, the \({\mathrm{FPT}}^{(+)}\) performs analytical continuation by inducing convergence into the divergent development input \(G_N(z^{-1})\) for \(N\rightarrow \infty .\)

We always simultaneously apply the \({\mathrm{FPT}}^{(-)}\) and \({\mathrm{FPT}}^{(+)}\) for the reason of cross-validation within the same type of signal processing. In the end of computations, after full convergence has been reached in the \({\mathrm{FPT}}^{(-)}\) and \({\mathrm{FPT}}^{(+)},\) only the common, joint reconstructions (spectral parameters, spectral lineshapes,...) are retained as the final, meaningful results. With this intrinsic checking of the output of the Padé-based estimation, chances are minor (if at all) to retain spurious information which, being random, is unlikely to be shared by both the \({\mathrm{FPT}}^{(-)}\) and \({\mathrm{FPT}}^{(+)}.\) There is no other signal processor with such a systematic, robust and reliable way of preserving true, while simultaneously discarding false and misleading information from data analyses.

2.6 Signal-noise separation and denoising Froissart filter

In the spectrum \(P_K/Q_K\) from the FPT, the numerator polynomial \(P_K\) gives the spectral zeros and suppresses noise by a moving average (MA) process. Likewise, the denominator polynomial \(Q_K\) yields the spectral poles. Moreover, the expansion coefficients of \(Q_K\) coincide with those of an auto-regressive (AR) process. The combination of AR and MA yields the auto-regressive moving average (ARMA) [1]. The zeros and poles of the complex-valued envelope \(P_K/Q_K\) are determined exclusively by \(P_K\) and \(Q_K,\) respectively, because the spectrum \(P_K/Q_K\) is a meromorphic function. A function having poles as its only singularities is called a meromorphic function. By its rational polynomial representation \(P_K/Q_K,\) the FPT is automatically of a polar structure, which naturally yields the peaks as the main signature of spectral envelopes.

Quite the contrary, the FFT, as a single Maclaurin polynomial, has no polar structure whatsoever. As a consequence, the FFT needs too many of its basis set functions to approximately mimic a polar structure, i.e. to describe a peak in a spectrum. According to Eq. (2.3), the basis set function in the FFT is the unattenuated, equidistantly discretized harmonic variable \(\exp {(-2\pi ink/N)},\) with N being the total time signal length, whereas n and k are the running numbers of the data in the time and frequency domain, respectively. The Fourier grid linear frequencies, \(k/T\, (0\le k\le N-1),\) being real-valued, are all located exactly on the circumference (\(|z|=1\)) of the unit circle. On the other hand, the natural accumulation limit of noise poles is the same circumference \(|z|=1\) of the unit circle. As such, the Fourier grid frequencies are completely embedded in noise and, hence, they are likely to often coincide with noise poles. Since the Fourier grid frequencies are real, they cannot reconstruct the sought complex eigen-frequencies of metabolites.

In the expansion methods, exponential convergence, \(\exp {(-M)},\) with respect to the size M of the basis set can be achieved only with basis functions that take into account singularities of the expanded function. This is automatically the case with the FPT due to its polar structure. As stated, poles are the only singularities of rational polynomials \(P_K/Q_K\) from the FPT. Moreover, the signatures of spectra are peaks. These stem directly from the singularities (poles) in any MRS spectrum. The mentioned undamped harmonics, as the basis functions in the FFT, do not describe the singularities of the expanded function \(G_N(z^{-1}).\) Therefore, convergence in the FFT relative to an increase of the basis size N (which is also the total signal length) can follow only the inverse power law (1 / N).

In sharp contrast to Fourier analysis, the parametric \({\mathrm{FPT}}^{(\pm )}\) reconstruct the complex eigen-harmonics \(\{z^{\pm 1}_k\},\) where \(z_k^{\pm 1}=\exp {(\pm 2\pi i\nu _k\tau )},\) whose imaginary parts shift the genuine poles away from the circumference \(|z|=1\) of the unit circle. Such a shifting is inside (\(|z|<1\)) or outside (\(|z|>1\)) the unit circle since \(|z|<1\) or \(|z|>1\) in the \({\mathrm{FPT}}^{(+)}\) or \({\mathrm{FPT}}^{(-)},\) respectively. The benefit is in significantly reducing overlap of genuine with spurious poles and this, in turn, improves signal-to-noise ratio (SNR). The gain is achieved under the provision that there is a way to minimize or possibly eliminate altogether the presence of spurious poles resulting from noise in the input FID and/or from round-off errors in signal processing. It is here that the concept of signal-noise separation (SNS) via identification of Froissart doublets comes to the rescue for the FPT.

Froissart doublets are pole-zero coincident pairs characterized by small amplitudes (like noise background). The word “doublet” signifies that a Froissart structure always appears as a pair (a couple comprised of a pole and a zero). Generally (i.e. irrespective whether we are dealing with spurious or genuine resonances), a pole-zero distance is proportional to the amplitude of the given component in the FID. That is why these amplitudes of time signals are alternatively called residues: they are proportional to whatever is left after subtracting the poles from the zeros.

This explains why Froissart doublets are feeble: a small pole-zero distance in a spurious couple produces a small amplitude. Since they are weak, they can be easily perturbed and, thus, distorted. Even small changes of some of the input parameters (e.g. slight truncations of the full signal length N) can lead to large alteration of Froissart doublets e.g. in their locations in the complex plane and/or strength. This makes them random, like noise. Such a marked instability classifies Froissart resonances as unphysical (spurious) spectral structure.

These spurious resonances are automatically removed from the Padé spectrum through its self-corrected feature: a spurious poles from \(Q_K\) cancels the spurious zero from \(P_K\) in the given Froissart doublet. Such a pole-zero cancellation of spurious Froissart pairs, as the result of the underlying pole-zero coincidence, occurs thanks to the rational form \(P_K/Q_K\) in the spectrum of the FPT. As such, detecting pole-zero coincidences amounts to filtering noise. This gives the name denoising Froissart filter (DFF). The DFF is unique to the FPT because the spectrum of this processor is a polynomial quotient \(P_K/Q_K\) where we have pole-zero cancellations due to the emergence of Froissart doublets as pole-zero coincidences.

On the other hand, there is part of the Padé spectrum \(P_K/Q_K,\) which is resilient to perturbations (such as truncation of the signal length N, etc). Therefore, these stable spectral structures are categorized as physical (genuine). They are comprised of genuine poles as well as genuine zeros that do not coincide. For such structures, there are no pole-zero cancellations in the spectral quotient \(P_K/Q_K.\) Thus, since the coincident, i.e. noisy poles and noisy zeros are gone through their cancellations, all that is left in the Padé spectrum \(P_K/Q_K\) is the true, physical information (genuine poles and genuine zeros). This is how the FPT improves SNR.

The gist of the matter is in the concept of signal-noise separation, SNS, in the FPT: this processor accomplishes its self-correction by reducing non-physical, noisy information (via pole-zero cancellation in Froissart doublets) and, at the same time, by retaining the physical signal (no pole-zero cancellations in genuine resonances) [14]. It is a paradigm shift in signal processing which can be given an alternative, transparent name “signal modulated estimation” (SME).

\(\bullet \) The SME is summarized as the following algorithm: optimize signal analysis through a modulation of estimation by maximizing extraction of true, while simultaneously minimizing the chance for unknowingly retaining the false information.

Such a weeding of the final linelist in data analysis is of paramount importance in all applications of signal processing to encoded FIDs. This is the case because corruption of physical by unphysical information adversely impacts on the interpretation of the measured data.

3 Spectral representations

The initial or default spectral representations in the \({\mathrm{FPT}}^{(\pm )}\) are given by the polynomial quotients \(P^\pm _K/Q^\pm _K.\) These representations are common to both non-parametric and parametric \({\mathrm{FPT}}^{(\pm )}.\) The parametric versions of the \({\mathrm{FPT}}^{(\pm )}\) have two additional representations given by the Heaviside expansions and the canonical forms. The Heaviside partial fraction decompositions read as:

$$\begin{aligned} \frac{P^\pm _K(z^{\pm 1})}{Q^\pm _K(z^{\pm 1})}=d^\pm _K+\sum \limits _{k=1}^K\frac{d^\pm _kz^{\pm 1}}{z^{\pm 1}-z^{\pm }_{k,Q}}. \end{aligned}$$
(3.1)

Here, the functions \(d^\pm _kz^{\pm 1}/(z^{\pm 1}-z^{\pm }_{k,Q})\) represent the component spectra for the \(k\,\)th molecule. On the other hand, the canonical representations are defined by:

$$\begin{aligned} \frac{P^\pm _K(z^{\pm 1})}{Q^\pm _K(z^{\pm 1})} =\frac{p^\pm _K}{q^\pm _K}\prod \limits _{k=1}^K\frac{z^{\pm 1}-z^{\pm }_{k,P}}{z^{\pm 1}-z^{\pm }_{k,Q}}. \end{aligned}$$
(3.2)

The quantities \(z^{\pm }_{k,P}\) and \(z^{\pm }_{k,Q}\) are the roots of the characteristic equations \(P^\pm _K(z^{\pm 1}_k)=0\) and \(Q^\pm _K(z^{\pm 1}_k)=0,\) respectively, with the simplified notations \(z^{\pm 1}_{k,P}\equiv z^{\pm }_{k,P}\) and \(z^{\pm 1}_{k,Q}\equiv z^{\pm }_{k,Q}.\)

Advantageously, the amplitudes \(d^\pm _k\) in the \({\mathrm{FPT}}^{(\pm )}\) are available in the analytical forms as the Cauchy residues of the spectral envelopes \(P^\pm _{K}(z^{\pm 1})/Q^\pm _{K}(z^{\pm 1})\) taken at the roots \(z^{\pm }_{k,Q}.\) They have the following two equivalent expressions:

$$\begin{aligned} d^\pm _k =\left\{ \frac{P^\pm _{K}(z^{\pm 1})}{(\text {d}/\text {d}z^{\pm 1})Q^\pm _{K}(z^{\pm 1})}\right\} _{z^{\pm 1}=z^{\pm }_{k,Q}}, \end{aligned}$$
(3.3)

and

$$\begin{aligned} d^\pm _k=\frac{p^\pm _K}{q^\pm _K}\frac{\prod _{k'=1}^Kz^\pm _{k,Q}-z^\pm _{k,P}}{\prod ^K_{k'=1,k'\ne k} z^\pm _{k,Q}-z^\pm _{k',Q}}. \end{aligned}$$
(3.4)

The peak heights are obtained from the component spectra \(d^\pm _kz^{\pm 1}/(z^{\pm 1}-z^{\pm }_{k,Q})\) for \(\nu \) taken at the chemical shift \({\mathrm{Re}}(\nu ^\pm _{k,Q})\) of the \(k\,\)th resonance and their explicit expressions are shown in the pertinent figures from the Results section.

4 Results

In this section, we illustrate the main features of the expounded theory of signal processing within the FPT. Both non-parametric and parametric estimations will be performed and the ensuing total shape spectra or envelopes will be compared with each other.

4.1 Characteristics of the encoded time signals

We shall use the time signal encoded at a General Electric magnetic resonance (MR) scanner with the static magnetic field strength \(B_0=1.5\)T corresponding to the Larmor frequency \(\nu _{\mathrm{L}}=63.87\) MHz. The encoding was made at the Astrid Lindgren Children’s Hospital (Stockholm) from the parietal temporal brain region from a child who had cerebral asphyxia. The total signal length of this FID is 512,  the bandwidth (BW) is 1000 Hz and the sampling time \(\tau =\) is \(\tau =1/\hbox {BW}=1\) ms. Single-voxel proton MRS was used in encoding this FID with the point-resolved spectroscopy sequence (PRESS). More precisely, as usual in MRS, some 128 such time signals were encoded and subsequently averaged to improve SNR (a single encoded FID would be too noisy to be of any practical use). Such an occurrence is alternatively referred to as 128 NEX (number of excitations). It is this averaged FID which is used in our signal processing. Altogether, three values of echo time (TE) were employed in encoding (\(\hbox {TE}=24\), 136 and 272 ms) alongside the repetition time (TR) of 2000 ms. The present illustrations deal only with \(\hbox {TE}=24\) ms since, in this case, the short lived resonances have not decayed and, thus, are advantageously the subject of reconstructions. They are often very important for overall interpretation. Zero-filling to 1024 FID points is used in the FFT and, for consistent comparisons, the same zero-padding is also employed in the FPT. As such, the total time duration of the zero-padded FID becomes \(T=1024\) ms. During the process of encoding, the giant water resonance was partially suppressed by means the standard inversion recovery.

4.2 Reconstructions in signal processing

Prior to reconstruction, the FID was phase corrected. This was done through multiplication of the complex-valued FID by the phase factor \(\exp {(i\phi _0)}\) where \(\phi _0=-1.1831\) rad. This particular value of \(\phi _0\) represents the minimum of the real parts of the complex FFT spectrum with the originally encoded, raw FID. As always, for cross checking, we apply both the \({\mathrm{FPT}}^{(-)}\) and \({\mathrm{FPT}}^{(+)}.\) However, for brevity, all the results will be reported for the \({\mathrm{FPT}}^{(+)}\) alone. As such, hereafter, it is sufficient to refer to \({\mathrm{FPT}}^{(+)}\) by FPT, for short. Specifically, for parametrically computed Padé spectra, it is useful to consider the so-called “Usual” (U) and “Ersatz” (E) spectra. The U- and E-spectral envelopes are computed from the standard Heaviside partial fraction decomposition with \(d_k\) (complex amplitude) and \(|d_k|\) (absolute value of complex amplitude), respectively. The absorption mode of an “Ersatz” spectrum is convenient because it is always positive-definite. We can say that in an Ersatz spectrum, phase of each reconstructed resonance is corrected thanks to the extraction of the pertinent phase (which may contain some distortions due to various factors in encoding of time signals). By default, the terms “Usual” and “Ersatz” refer only to the parametric FPT. No such distinction can be made in the non-parametric FPT because this signal processor does not encounter \(d_k\) as quantification is not performed. As such, naturally, the present non-parametrically produced total shape spectra will not have the acronyms U (“Usual”) nor E (“Ersatz”).

We first address the effect and significance of the procedure called spectra averaging. As stated, it would not make sense to encode a single time signal. Likewise, it may not be of much practical use either to complete computations for a single model order K. In both cases, the reason is the same: presence of too much noise. For a single FID, intolerably large noise comes from various sources of noise (mostly of uncontrollable, i.e. random nature). Even after yet another partial suppression of the residual water by means of some theoretical procedures, a spectrum computed by any signal processor may also be heavily corrupted with noise (at least at some frequency intervals) which stems from both encoding and reconstructions of spurious resonances. All such noise is random, meaning that it may widely fluctuate with even minor external changes. To mitigate the influence of noise, averaging is performed in the time and frequency domain. As mentioned, in the time domain, many FIDs are encoded and then averaged. Likewise, in the frequency domain, many spectra are computed for different values of K and the results are averaged.

As an example, we first computed non-parametrically some 31 total shape spectra or envelopes for K varying in the interval \([K_{\mathrm{min}},K_{\mathrm{max}}]\) where \(K=K_{\mathrm{min}}+\Delta K\) with \(\Delta K=1,\, K_{\mathrm{min}}=385\) and \(K_{\mathrm{max}}=415\) (i.e. \(K=385,386,\ldots ,415\)). Model order K and the partial signal length \(N_{\mathrm{P}}\) are connected by the relation \(N_{\mathrm{P}}=2K.\) As such, the values of the partial signal length belong to the interval \(N_{\mathrm{P}}\in [770,830].\) These are called partial signal lengths due to the implied truncation \((N_{\mathrm{P}}=770, 772, \ldots , 830)\) of the full signal length, \(N=1024.\)

Fig. 1
figure 1

Spectra averaging in the 1st iteration. Panel (a) shows the real parts of 31 envelopes (total shape spectra) in the non-parametric FPT for model orders K in the interval \(K\in [385,415]\) corresponding to the partial signal lengths \(N_{\mathrm{P}}\in [770,830].\) These partial signal lengths are truncations of the full signal length \(N=1024.\) The first half of N contains the encoded FID. The remaining 512 entries are the added real numbers with zero-valued intensities, i.e. they are zeros (one zero-filling from 512 to 1024). Panel (b) displays the real part of the 1st average complex envelope. We took the arithmetic average of the 31 complex envelopes from the lineshape estimations by the non-parametric FPT with \(K\in [385,415].\) This implies that for real-valued spectra, the curve on panel (b) is the average of the curves from panel (a). It is observed that the spurious, noisy spikes from panel (a) are smoothed out on panel (b) by the process of spectra averaging (Color figure online)

Figure 1 illustrates spectra averaging. Panel (a) of this figure exhibits the differences in 31 envelopes of varying model orders K. Some fluctuations (small or large) are present almost throughout the displayed frequency region of interest (ROI). Most notable changes, appearing as sharp spikes with large amplitudes and thin linewidths, are around 2.55 and 4.3 ppm. All these random fluctuations practically disappear from the average spectrum depicted in panel (b) of the same figure. Simultaneously, due to their robust stability (resilience to noise perturbations), physical or genuine resonances are seen in Fig. 1b to survive averaging.

Mathematically, the spectrum in Fig. 1b corresponds to the arithmetic average of 31 spectra from Fig. 1a. Note the presence of the labels “The 1st iteration” and “The 1st average” on panels (a) and (b), respectively. This is written in anticipation that the process of generating an average spectrum may be repeated to encompass higher iterations (2nd, 3rd,...), as will indeed be the case in the analysis which follows. Observe also that the acronyms written for a number of peaks in Fig. 1b abbreviate the names of the molecules present in the examined substance, which is a slice excited by the applied pulses in the radiofrequency (RF) part of the electromagnetic radiation.

The reason for which the fluctuations in lineshapes for varying K are mitigated in Fig. 1b is simply in the fact that SNR is improved by a factor \(\sqrt{M}\) by averaging where M is the number of individual envelopes (\(M=31\) in the present computations). This is a general feature of improvements of SNR when it comes to random perturbations. However, the present mechanism for reduction of randomness (noise and noise-like features) is transpired by pole-zero cancellations. Recall here that spurious (noisy) resonances appear in pairs (Froissart doublets) that are washed out by pole-zero cancellation in the Padé spectrum. In principle, such cancellations should occur for each value of K. If this was indeed the case, no fluctuation and spurious resonances (spikes) should, in fact, be seen throughout the entire Nyquist range of frequencies, including the ROI from Fig. 1a. This, however, is not observed in Fig. 1a. The reason is twofold. First, unlike noiseless synthesized time signals where pole-zero cancellation is exact, noise in any encoded FID precludes exact pole-zero cancellation. Consequently, the less exact pole-zero cancellation, the more spike-looking the envelope. Second, the lack of pole-zero cancellations for noisy spectral structures may also partially be due to those values of K for which the Padé spectra have not yet fully converged.

Fig. 2
figure 2

Non-parametric versus parametric lineshape estimations. Panel (a) is a copy borrowed from Fig. 1b (the real part of the 1st average complex envelope from the non-parametric FPT). Panels (b) and (c) concern the non-parametric and parametric FPT, respectively, at a single model order \(K=512\, (N_{\mathrm{P}}=1024),\) i.e. at the full signal length. Both processors employ the same FID generated by the IFFT-inverted 1st average complex envelope from the non-parametric FPT. It is seen on panels (b) and (c) that the non-parametric and parametric FPT give the same envelope. Hence self-cross-validation within the FPT (Color figure online)

We saw that Fig. 1 deals exclusively with non-parametric signal processing. To complement and cross-validate this information, Fig. 2 compares the envelopes computed by both the non-parametric and parametric estimations. First, on Fig. 2a, we show again the real part of the 1st average complex envelope (as in Fig. 1b). Second, on Fig. 2b, we show the real part of a new complex total shape spectrum computed non-parametrically at a single value of model order, \(K=512\) (i.e. using the full signal length \(N_{\mathrm{P}}=N=1024\)). The said new complex spectrum is generated using the FID obtained by inversion of the average complex envelope whose real part is in Fig. 1b. Here, inversion is performed by means of the standard inverse fast Fourier transform (IFFT). The new FID is afterward subjected to the Padé-based processing and the results are shown in panels (b) and (c) of Fig. 2 for non-parametric and parametric estimations, respectively. It is observed that these two predictions are in full agreement throughout the spectrum. When plotted together on the same graph (not shown to avoid clutter), the curves from panels (b) and (c) would be indistinguishable from each other. This outcome of the effected comparison, retrospectively validates quantification (reconstruction of peak positions, widths, heights, phases) performed by the parametric FPT.

As elaborated, Fig. 1 illustrates the effect and significance of reduction of a marked sensitivity to changes in model order K. Therein, a substantial reduction of this sensitivity is seen already with the 1st arithmetic average. Computation-wise, such a finding is of a considerable practical importance. Nevertheless, this does not imply that the analysis should stop at the 1st average spectrum from Fig. 1b without any further checking. Quite the contrary, it is prudent to go beyond the 1st average and perform iterations. In other words, for more comprehensive verification purposes, it would be important to have a follow-up iterative process, as anticipated earlier. This would permit seeing what actually remains after repeated averaging as a more robust sign of the stability of the whole reconstructed physical information (due to all the genuine resonances) from the processed FID.

Fig. 3
figure 3

The 2nd average envelope. Panel (a) is a copy borrowed from Fig. 1b (the real part of the 1st average complex envelope from the non-parametric FPT). Using the non-parametric FPT, panels (b) and (c) focus on the 2nd iteration and the associated 2nd average envelope, respectively. The 2nd iterative envelopes with \(K\in [385,415]\) are generated through the FID obtained by the IFFT inversion of the 1st average complex envelope. The resulting 31 real-valued envelopes on panel (b) exhibit spurious spikes only around 3.4, 4.2 and 4.35 ppm. This demonstrates stability of spectra averaging through the 2nd iterative step (Color figure online)

With this goal, we performed a number of iterations by extending further the procedure from Fig. 1. The 1st iteration from Fig. 1a consists of applying the FPT to the originally encoded FID (once zero filled) for a sequence of 31 values of model order K. Using the resulting 31 complex envelopes, their real parts are plotted in Fig. 1a and the ensuing real average envelope is re-displayed in Fig. 3a for further comparisons. What could be a part of the 2nd iteration has already been accomplished in Fig. 2b, but only for a single value of model order \((K=512).\) This time, however, we generalize the process from Fig. 2 by extending it to a sequence of values of K. Specifically, the reconstructed FID stemming from the IFFT-inverted 1st average complex envelope is subjected to the repeated Padé non-parametric shape estimation at the same interval of model orders as before, \(K\in [385,415].\) This constitutes the 2nd iteration which gives the new 31 complex envelopes whose real parts are shown in Fig. 3b. Such real-valued envelopes exhibit by far much less intense fluctuations in the lineshapes for the individual values of K than those in Fig. 1a. Averaging the complex envelopes from the 2nd iteration yields the 2nd average complex envelope whose real part is depicted in Fig. 3c. It is seen that the 1st and the 2nd real average envelopes from Fig. 3a and 3c are in excellent agreement.

Fig. 4
figure 4

Iterative spectra averaging with the first three iterations. Panel (a) is a copy borrowed from Fig. 1a (31 envelopes from the 1st iteration). Panels (b) and (c) represent the 2nd and the 3rd iterations, respectively. All the three panels (a), (b) and (c) deal with 31 envelopes generated by the non-parametric FPT for the same range of model order K,  i.e. \(K\in [385,415]\). The 2nd and the 3rd iterative envelopes are produced by using the FIDs obtained by the IFFT inversion of the 1st and the 2nd average complex envelopes, respectively. Marked spikes in the 1st iteration (panel a) are significantly reduced in the 2nd iteration (panel b), as already seen in Fig. 3b. The 3rd iteration (panel c) has no spikes at all. The corresponding real parts of the 1st, 2nd and 3rd average complex envelopes are depicted on panel (d). These three curves coincide with each other, showing robust stability in spectra averaging through the 3rd iterative step (Color figure online)

Figure 4 includes the 3rd iteration and the ensuing 3rd average spectrum. For a parallel presentation and comparisons of relative performance, panels (a) and (b) re-display the spectra due to the 1st and the 2nd iterations borrowed from Figs. 1 and 3. The 3rd iteration from panel Fig. 4c begins by first generating the next new FID. Prior to producing that FID, we perform arithmetic averaging of complex envelopes from the 2nd iteration (as we did in Fig. 3c). Afterward, this 2nd average complex envelope is IFFT-inverted to obtain the 2nd new complex FID. Finally, such a newly reconstructed FID is subjected to Padé non-parametric processing for \(K\in [385,415]\) yielding the 3rd iteration with 31 complex envelopes whose real parts are given on panel (c) of Fig. 4. Therein, even the last remaining fluctuations around 3.4, 4.2 and 4.35 ppm seen earlier in Fig. 3b have now disappeared altogether. There are no more noticeable spikes in any of these 31 envelopes. The new 31 complex envelopes from the 3rd iteration are averaged to yield the 3rd average complex total shape spectrum.

Finally, Fig. 4d compares the real parts of the 1st, 2nd and 3rd average complex envelopes. The ensuing three curves are practically indistinguishable from each other. This shows the stability of the reconstructed physical information through the 3rd order in the explained iterative process. Moreover, the same conclusion holds true also for further iterations, as we explicitly checked using every successive iteration of the orders \(4,5,\ldots ,12\). Thus, for a further analysis, it is safe to retain the 3rd average complex total shape spectrum.

In the illustrations from Figs. 1 and 3, we dealt exclusively with the non-parametric Padé processing. This concerns lineshape estimation alone. The exception is Fig. 2c showing the real part of the complex envelope generated by the parametric Padé processing. This latter parametric envelope employs the 1st new FID reconstructed from the 1st average complex envelope which is generated non-parametrically and whose real part is shown in Figs. 1b and 2a. Figure 2c was restricted to the analysis of lineshapes alone and, as such, carried no information about the underlying reconstructed parameters (via the peak positions, widths, heights, phases).

Fig. 5
figure 5

Component spectra from the 3rd average envelope. Panel (a) shows the real part of the 3rd average complex envelope (also present in Fig. 4d) from the non-parametric FPT. The corresponding 3rd average complex envelope is IFFT-inverted to deduce the reconstructed FID which is afterward subjected to the parametric FPT at a single model order \(K=512\, (N_{\mathrm{P}}=1024)\) associated with the full signal length \(N=1024\). Padé-based quantification with this FID yields the peak parameters. From these the complex components are generated in the “Usual”, U, and “Ersatz”, E, modes. The real parts of the complex U- and E-components are displayed in panels (b) and (c), respectively. These panels also show the corresponding peak heights. Absorption and dispersion components are mixed together in the U-spectra (panel b). Only pure absorptive, symmetric component Lorentzians appear in the E-spectra (panel c) (Color figure online)

However, Fig. 5 deals with the said peak parameters stemming from the Padé-based solution of the quantification problem. To quantify the information hidden in the input FID, reconstruction is done by parametric processing in the FPT. To this end, given that the 3rd iteration has fully stabilized, as stated earlier, we employ the 3rd new FID. This time signal results from the IFFT inversion of the 3rd average complex envelope which resulted from the non-parametric FPT. The real part of the 3rd average complex envelope due to the non-parametric FPT is shown in Fig. 5a (a repetition of one of the three curves in Fig. 4d). We then subject the 3rd reconstructed FID to the parametric FPT at \(K=512\) by using the full signal length (\(N_{\mathrm{P}}= N =1024\)). The results of quantification give all the sought peak parameters. These parameters are then used to generate the component lineshapes as the constituents of the envelopes.

Two kinds of such component spectra are shown in Fig. 5. These are the Usual, U, and Ersatz, E, component lineshapes. As stated, the U- and E-spectra are obtained by employing the complex amplitudes \(d_k\) and their absolute values \(|d_k|,\) respectively. The real parts of such complex component U- and E-spectra are shown in Fig. 5b, c, respectively. It is seen that the absorption and dispersion modes are mixed together in the component U-spectra from Fig. 5b. On the other hand, the component E-spectra from Fig. 5c are all in the purely absorptive mode and, hence, positive-definite.

Also displayed in panels (b) and (c) of Fig. 5 are the peak heights (symbolized by open circles) retrieved from the analytical expressions listed in this figure. Peak heights in the component U-spectra are generated from the complex-valued amplitude \(d_k\) of the \(k\,\)th harmonic in the FID. This is the reason for which the peak heights in the component U-spectra are phase modulated and, thus, appearing as \(|d_k|\cos {(\varphi _k)}\) where generally the phase \(\varphi _k\) takes on non-zero values (\(\varphi _k\ne 0\)).

This phase modulation of the magnitude \(|d_k|\) is present in the expression for the peak height written on the title of panel (b) in Fig. 5. As seen in Fig. 5b, such phase modulations of \(|d_k|\) cause that the open circles for the peak heights are close to the abscissa, where they become nearly zero-valued for \(\varphi _k\approx \pi /2\, (\cos {\varphi _k}\approx 0)\) in the case of dispersive lineshapes. On the contrary, some of the nearly absorptive U-component lineshapes in Fig. 5b, have their maximae closely coinciding with the associated peak heights because, in this case, \(\varphi _k\approx 0\, (\cos {\varphi _k}\approx 1)\) and, thus, no phase modulation occurs, i.e. \(d_k\approx |d_k|\ne 0\).

Of course, no phase modulation ever occurs in the component E-spectra on Fig. 5c since in these lineshapes the expressions for the peak heights deal with \(|d_k|\) from the onset (see the corresponding expression written in the title of Fig. 5c). Namely, in the component E-spectra, each complex amplitude, \(d_k=|d_k|\exp {(i\varphi _k)}\) is phase corrected through multiplication by the factor \(\exp {(-i\varphi _k)}.\) This makes the phase corrected amplitude \(d^{{({\mathrm{phase}}\, {\mathrm{corrected}})}}_k\) positive-definite via \(d^{{({\mathrm{phase}}\, {\mathrm{corrected}})}}_k=d_k\times \exp {(-i\varphi _k)}=|d_k|\exp {(i\varphi _k)}\times \exp {(-i\varphi _k)}=|d_k|>0.\) As a consequence, the open circles for the peak heights of the component E-spectra coincide exactly with the maximae of the purely absorptive lineshapes in Fig. 5c.

Peak heights from Fig. 5b, c depend on the peak widths as clear from the discussed analytical expressions listed in these two panels. Peak widths are connected to the imaginary parts of the reconstructed fundamental frequencies. The real parts of these complex frequencies are not explicitly contained in the expressions for the peak heights. However, each peak height is associated indirectly with the real parts of the complex fundamental frequencies, as well. This occurs because it is at the positions of the real parts of the complex fundamental frequencies that the maximae of e.g. absorptive Lorentzian lineshapes are located. To complement this indirect information, it is important to plot explicitly e.g. the magnitudes \(|d_k|\) as functions of the real fundamental frequencies. This type of graph is called the magnitude plot. Likewise, to illustrate the relationship between the real and complex fundamental frequencies, it would be useful to display the corresponding Argand plots. These plots or diagrams show the imaginary versus real fundamental frequencies on the ordinates and abscissae, respectively.

Fig. 6
figure 6

Signal-noise separation by way of the denoising Froissart filter, DFF (pole-zero cancellations). The spectral parameters reconstructed by the parametric FPT and used in Fig. 5 are here displayed as the Argand plot (panel a) and magnitude plot (panel b). On panel (a), the reconstructed imaginary frequencies are shown as the function of the corresponding real frequencies. Therein, spurious resonances are identified by pole-zero confluences (open circles, as spectral poles, coinciding with full dots, as spectral zeros). These are Froissart doublets describing noise. Such coincidences lead to zero-valued amplitudes (panel b) for noisy resonances. Consequently, the associated zero-valued peak heights also follow, as seen in the component “Ersatz” spectra (panel c). Pole-zero coincidences and the ensuing pole-zero cancellations represent the mechanism by which the underlying denoising Froissart filter, DFF, accomplishes signal-noise separation, SNS (Color figure online)

These discussed relationships are the topics on Fig. 6 where the Argand and the magnitude plots are shown on panels (a) and (b), respectively. In particular, two kinds of Argand plots are contained in Fig. 6a, one for the poles and the other for zeros. As discussed in the Theory section, these zeros and poles are the respective roots of the numerator (\(P_K\)) and denominator (\(Q_K\)) polynomials in the Padé spectrum \(P_K/Q_K.\) To correlate the parametric information from the Argand and magnitude plots to the spectral lineshapes, panel (c) of Fig. 6 shows the absorption lineshapes of the component E-spectra alongside with the pertinent peak heights.

The essence of Fig. 6 is grasped by two kinds of behavior of the peak parameters. On panel (a), coincidence or near-coincidence of poles and zeros in the complex frequency plane (as one of the two signatures of the appearance of Froissart doublets) indicates the presence of spurious, noisy information. The other signature of Froissart doublets is a set of zero or near-zero magnitudes \(|d_k|\) of complex amplitudes \(d_k\) of spurious resonances, as seen on panel (b) in Fig. 6. Overall, Fig. 6 illustrates the concept of signal-noise separation, SNS. As expounded in the Theory section, physical signal is separated from noise by identification of Froissart doublets. Once identified, noise which is disguised as Froissart doublets can be discarded from the linelist of the reconstructed parameters. This amounts to signal-noise separation, SNS, which in turn improves signal-to-noise ratio, SNR. In fact, this reduction or ultimately elimination of noise is achieved in a self-contained manner within the Padé spectrum \(P_K/Q_K\) because therein pole-zero coincidences automatically yield pole-zero cancellations. Hence, the denoising Froissart filter, DFF, is a reliable self-correcting procedure for retaining the true, physical information alone reconstructed from time signals encode by MRS.

5 Discussion and conclusions

This study covers several important grounds. It is within the realm of data analysis and interpretation by the methodology of reconstruction of information hidden in the measured input data. Mathematically, it deals with a difficult ill-conditioned (ill-posed) problem, called spectral analysis. This particular problem is alternatively called the quantification problem in nuclear magnetic resonance spectroscopy (or magnetic resonance spectroscopy, MRS, for short) from analytical chemistry. It is also known as the harmonic inversion problem.

In mathematics and quantum physics, spectral analysis solves the given eigen-value problem. The solution is a set of eigenvalues and eigenfunctions. When the underlying operator in an eigen-problem is a Hamiltonian (e.g. the sum of the kinetic and potential energy operators, as in physics) the eigen-values are the characteristic energies (eigen-energies or fundamental energies) representing a set of the possible states of the given system. In such a case, the absolute values of the squared complex scalar or inner products of the eigen-functions with the initial state of the system give the intensity of the lines in the energy spectrum. The eigen-energy spectrum is comprised of sticks or bars with zero-valued line widths for purely discrete stationary stable states. On the other hand, resonances of the system have non-zero line widths associated with finite lifetimes.

Quantification in MRS deals with resonances of molecules that are present in the examined substance or specimen. Solving this problem amounts to extraction of the unknown spectral parameters. These are the quantifiers of resonances or peaks in the frequency spectrum. Each peak is completely determined by reconstructing its position (location in the frequency spectrum), width, height and phase. Such peaks are associated with the components (complex harmonic functions) in the input data that are the time signals encoded by MRS from the examined specimen. This gives the name, harmonic inversion. Herein, an inverse problem is to be solved. It consists of retrieving the unknown components of the input time signal. Namely, a set of the digitized values of the encoded time signal is known and the task is to uncover the constituent components of these input data. This problem is ill-conditioned because of the lack of its continuous dependence on the input data points. As a result, even small perturbations in the input data can lead to large deviations in the reconstructed, output data. Noise in the encoded time signals is a prime example of such ubiquitous distortion of encoded time signals.

The main task in MRS is to improve resolution and reduce noise. Improved resolution would increase the specificity of magnetic resonance. Magnetic resonance imaging, MRI, has good sensitivity, but its specificity is insufficient. Optimized MRS would come to the rescue of the entire MR methodology. This could be achieved by improving specificity through reliable quantification of the physical information from the encoded data.

However, solving the quantification problem without knowing how to reliably separate physical signal from noise belongs to the category of guess work that cannot be of any practical use. Unfortunately, this is precisely the case with those signal processors in MRS that are based on fitting. Such fitting techniques usually employ total shape spectra (envelopes) from the fast Fourier transform, FFT. Another type of fitting the given spectrum is to perform adjustments of parameters in some linear combination of model spectra for each molecule assumed (prior to processing) to be present in the investigated substance.

Fitting approaches with any constraints are equivocal and non-unique. This occurs because the same least-square errors can be obtained by vastly different sets of the adjusted free parameters. Moreover, fitting procedures are always under- and/or over-modeling. In under-fitting, some of the true information is missed. In over-fitting, some of the wrong/misleading information is misclassified as true. Both cases impact adversely on whatever is thought to be worthwhile retaining after completion of fitting. This implies that the fitting parameters are either under- or over-estimated. When some of the genuine resonances have not been retrieved, or spurious resonances erroneously declared as true information, the constraints imposed on fitting stretch or squeeze particularly the peak heights and peak widths. The latter two parameters (assuming that the peak positions, i.e. chemical shifts, have properly been located) are the main quantities for determination of molecular concentrations. In other words, inadequate values of peak heights and peak widths would invalidate estimates of abundance of molecules. This amounts to failing to accurately reconstruct the molecular concentrations as the main information sought in MRS.

We pursue an alternative way by using an advanced signal processor without resorting to any fitting. It is the fast Padé transform, FPT. Its chief feature is non-linearity seeded in the mathematical form of the frequency spectrum given by quotient of two polynomials, \(P_K/Q_K.\) Rational polynomials are the best model functions in the mathematical branch called theory of approximations. The reason for their topping the list of all the existing approximations is primarily in the polar structure.

Poles as the only singularities of the meromorphic functions \(P_K/Q_K\) capture at once all the fundamental features of generic systems: the eigen-frequencies. The eigen-frequencies, or fundamental frequencies, are the trademarks of each system. They embody the features that make one system differ from the other. These eigen-frequencies or eigen-energies describe the states of the system with its characteristic oscillations of various degrees of freedom (e.g. vibrations, rotations of molecules, etc). The roots of the characteristic or secular equation \(Q_K=0\) give the sought fundamental frequencies. The corresponding intensities of the spectral lines are provided by an analytical expression for the Cauchy residues of the polynomial ratio \(P_K/Q_K\) taken at the eigen-frequencies from the said secular equation. It is in this straight and eminently clear manner that the parametric FPT completes exact quantification, the main goal of MRS. This holds true for synthesized and encoded time signals.

Mathematically, K is the polynomial degree which is also the model order. Physically, K represents the total number of reconstructed resonances. The number K is unknown prior to signal processing and this parameter is also unambiguously recovered by the FPT. The way this is accomplished is ingrained in signal-noise separation.

The true number of resonances is invariably larger than the estimated K. This is a direct consequence of the necessary over-determination of the system of equations for generation of the expansion of the numerator \(P_K\) and denominator \(Q_K\) polynomials. To eliminate the excess in the found K,  spurious resonances due to noise must be discarded from the output. This is accomplished by the denoising Froissart filter, DFF, which is the mechanism for signal-noise separation, SNS.

Spurious, noisy resonances emerge in pairs or doublets. They have first been detected by Froissart in numerical experiments with simulations. Hence the name Froissart doublets for spurious resonances associated with noise. Such random, noisy resonances are spotted by their pole-zero coincidences. Spurious zeros from \(P_K\) cancel spurious poles from \(Q_K\) because they appear in the spectrum as the polynomial ratio \(P_K/Q_K.\) This is the way in which pole-zero cancellation (as a self-correcting added value unique to the FPT) performs signal-noise separation, SNS. By this concept unstable, spurious resonances disappear altogether, whereas stable, genuine resonance stay. Subtracting the number of spurious resonances from the estimated total number K of resonances (genuine + spurious) yields the final, true number of molecules present in the specimen scanned by MRS.

The discussed features are internally validated within the two complementary and equivalent variants of the FPT, one initially convergent inside, and the other outside the unit circle in the plane of the complex harmonic variable. However, by the Cauchy analytical continuation, both variants converge everywhere in this complex plane with the exception of poles. The final linelist of spectral parameters is completed only after this cross-validation which retains the quantities jointly reconstructed by both versions of the FPT.

The parametric FPT improves resolution by splitting apart overlapping resonances. This is evidenced by the predicted component spectra generated from the reconstructed parameters. With this virtue at hand, alongside the concept of SNS (via DFF), the FPT becomes capable of simultaneously improving resolution and suppressing noise. These are the two major stumbling blocks that so far hampered progress of MRS. As such, it would seem that there could be a little doubt as to which signal processor should be most suitable for quantification in MRS.

All the expounded characteristics of the FPT are comprehensively illustrated in this work. The presented example is focused on exact reconstruction of spectral parameters from a time signal encoded by MRS with 1.5T scanner at short echo time, \(\hbox {TE} = 24\) ms. We here chose this short TE because in such a case most of short-lived resonances have not decayed yet. This offers an opportunity to extract maximal information since often short-lived resonances are most critical in interpretation of the encoded data. Another reason for selecting to analyze a time signal encoded at a short TE is the presence of many overlapping resonances in the frequency spectrum. With all the overlapped resonances resolved, the parametric FPT becomes the gold standard for testing performance of other estimators, including the derivative fast Padé transform [5,6,7, 42].

The present specific application of the FPT concerns MRS time signals encoded in vivo from a child who had suffered cerebral asphyxia. One of the reasons for choosing this problem is that MR spectra associated with asphyxia or ischemia are exceedingly dense and have heretofore been very difficult to resolve and interpret. The clinical context within MRS for pediatric, as well as adult neurodiagnostics is another reason for focusing on this particular problem.

Within neurodiagnostics, in vivo MRS has been primarily based upon a small number of metabolites and their concentration ratios. Among these is N-acetyl aspartate, NAA, resonating at approximately 2.0 ppm, which indicates the abundance and viability of neurons. With cerebral ischemia or asphyxia, the concentration of NAA is reduced [56]. This is due to the marked vulnerability of cerebral neurons to hypoxia. However, NAA can also be lowered with almost any damage to the brain, including brain tumors [57].

Another diagnostically important metabolite is choline, Cho, resonating around 3.2 ppm, which reflects phospholipid metabolism of cell membranes, and is a marker of membrane damage, cellular proliferation and cell density. In cerebral ischemia/hypoxia, a lactate, Lac, doublet, centered at 1.3 ppm, is expected, related to the predominance of anaerobic glycolysis. However, lactate can similarly be observed in brain tumors and sometimes in healthy brain tissue [56, 58,59,60]. Lipids, Lip, also resonating at around 1.3 ppm, often appear with reperfusion after hypoxia. Cerebral energy metabolism is reflected by creatine, Cr, which resonates at 3.0 ppm. Concentration of Cr in the brain is usually stable after the first year of life [60].

In the pediatric population, brain metabolite concentrations as well as metabolite concentrations ratios, depend upon the age. Myoinositol, m-Ins, is the dominant brain metabolite in neonates. In older infants, Cho is usually the largest peak. As the child’s brain matures, Cr and NAA concentrations increase. Concordantly, Cho to NAA and Cho to Cr concentration ratios normally fall with the child’s age [61].

One goal of the present study is to enhance neurodiagnostics through MRS, aiming towards greater accuracy in identifying cerebral hypoxia/ischemia versus other pathology, including brain tumors, which present differential diagnostic dilemmas. Broader implications for cancer diagnostics further motivate this investigation. Especially, it should be noted that hypoxic regions often occur within tumors. These regions are particularly resistant to radiation therapy as well as to chemotherapy. Moreover, hypoxia promotes genomic instability and is associated with the invasive/metastatic process [62]. Consequently, identifying hypoxic regions via MRS could contribute overall to better cancer treatment planning. Notably, the cancer biomarker phosphocholine, PC, also reflects hypoxia [63].