Synergism of spectra averaging and extrapolation for quantification of in vivo MRS time signals encoded from the ovary

Abstract

The fast Padé transform (FPT) is further optimized for encoded in vivo MRS time signals. This is achieved by a judicious combination of spectra averaging and time signal extrapolation. The motivation is in strengthening suppression of the over-sensitivity of signal processing to changes in model order K. Implementation is carried out in the FPT variant $\hbox {FPT}^{(+)}$ which, by numerically performed analytical continuation, converts divergent into convergent series. Convergence of reconstructions is monitored for a sequence of successive values of K. Comparison is made with the corresponding retrieval without spectra averaging and time signal extrapolation. Variances are dramatically reduced for the reconstructed parameters (complex frequencies and complex amplitudes) when spectra averaging and extrapolation are performed. Negligible variances imply convergence, which is accomplished herein using a single averaging procedure (no iterations). This has important implications in practice with encoded MRS data, providing remarkable efficiency, robustness and accuracy of Padé-based quantification. Algorithmically, spectra averaging and time signal extrapolation consist of four steps. First, the encoded time signal is used to compute total shape spectra (envelopes) for a sequence of model orders K, which is the number of resonances. Second, these envelopes are averaged (spectra averaging). All the spectra are computed at the same sweep frequencies whose number considerably exceeds the number of data points in the encoded time signal. Third, the complex average envelope is inverted to yield a new time signal which is longer than the encoded data (time signal extrapolation). Fourth, the Padé extrapolated time signal is quantified for a sequence of K to monitor convergence of the reconstructed parameters. All four steps are applied to in vivo MRS data encoded on a 3 T scanner from a borderline serous cystic ovarian tumor. We focus on this clinical problem to help develop effective methods for early detection of ovarian cancer, in order to improve survival for women afflicted with this malignancy. It is anticipated that the presently refined, multi-purpose Padé-based methodology with its practical advantages for in vivo MRS can contribute to this goal.

1 Introduction

Advanced signal processing methods are indispensable for handling data encoded via magnetic resonance spectroscopy (MRS). These methods are particularly important for early detection of ovarian cancer, for which the need is very great. In this paper we apply the fast Padé transform (FPT) to in vivo MRS time signals encoded from the ovary. We focus upon the combination of spectra averaging and time signal extrapolation as a further optimization of this advanced signal processing method for practical applications related to the mentioned public health problem of utmost importance. Before embarking on a detailed presentation of the mathematical features of the procedures for signal processing in MRS, we briefly review the medical aspects of the problem.

1.1 The clinical problem: the need for accurate and timely detection of ovarian cancer

Ovarian cancer is the sixth most frequently occurring malignancy among women worldwide and even more frequent in the USA, Scandinavia and Israel. Moreover, there appears to be a trend towards increasing incidence of this malignancy in many parts of the world [1,2,3,4]. The case fatality rate for ovarian cancer is very high; more than 14 000 women die each year from this disease in the USA alone [5, 6].

The high mortality rate is mainly due to late detection, at Stage III or IV with tumor spread outside the true pelvis [7]. When detected early, ovarian cancer has an excellent prognosis, especially when confined to a single ovary (Stage Ia) with five-year survival rates over 90% [8]. The difficulty is that early stage ovarian cancer is very frequently silent, i.e. without any symptoms. Further obstacles arise due to the small size of the ovary (6.1 to 1.8 cm$^{3}$, normally diminishing with age) and that it is not necessarily enlarged with early-stage cancer [9].

The most generally utilized diagnostic methods to screen for ovarian cancer are transvaginal ultrasound (TVUS) and serum cancer antigen (CA-125).¹ While there is some recent evidence that may possibly be to the contrary [10], the general consensus from large-scale randomized trials is that the combination of CA-125 and TVUS to screen for ovarian cancer in asymptomatic women does not reduce mortality nor contribute to earlier ovarian cancer detection and that this strategy is associated with a large percentage of false positive findings [11, 12]. The latter have a number of deleterious consequences, the most serious of which is that many women will undergo surgical intervention for benign ovarian lesions [13,14,15]. Thus, the overall view is that for women who are not at clearly high risk for ovarian cancer, the “harms” of routine screening for ovarian cancer outweigh the benefits [16]. Although the possibility of employing biomarkers other than CA-125 has been examined [17,18,19,20,21], none of these appear to sufficiently improve diagnostic accuracy to recommend their use for routine ovarian cancer screening [22].

Magnetic resonance imaging (MRI), through its high spatial resolution, has sometimes been found to be helpful in distinguishing benign from cancerous ovarian lesions that are indeterminate on TVUS [23,24,25]. Still, almost one-fourth of benign ovarian lesions were reportedly misdiagnosed as cancerous using TVUS with MRI as a second imaging technique [23]. In general, MRI is very sensitive, such that ovarian cancer will usually be detected, but many benign lesions will also be misdiagnosed as malignant. Diffusion-weighted imaging (DWI) can provide some additional help in assessing adnexal lesions, but many false positives also occur with DWI [26, 27].

Magnetic resonance spectroscopy, MRS, offers the possibility of going beyond anatomy, to evaluate the metabolic features of tissue or organs. Thereby, the molecular changes reflecting the cancer process, i.e. “hallmarks of cancer” [28] can potentially be revealed [29]. For almost two decades, it has been suggested that MRS could be the method of choice for early ovarian cancer detection [30, 31]. That potential has not yet been realized, however, due mainly to the limitations of the manner in which MRS time signals are analyzed. We now proceed to outline the specific limitations of the conventional Fourier-based signal processing method for MRS and to list the gains from an advanced approach to handling MRS time signals, as offered by the fast Padé transform, FPT.

1.2 Signal processing procedures in MRS

1.2.1 The conventional Fourier-based approach

It should be recalled that currently all clinical magnetic resonance (MR) scanners employ the fast Fourier transform (FFT) for converting the encoded time signal or free induction decay (FID) curve into its spectral representation in the frequency domain. The formula for generating the Fourier spectrum is:

$$\begin{aligned} \mathrm{FFT}: F_m =\sum _{n=0}^{N-1} {c_n \exp (-2\pi imn/N)}, \quad 0 \le m \le N - 1. \end{aligned}$$

(1)

This is seen to be a single polynomial, with the fixed mth Fourier grid frequency $2\pi m/T\left( {m=0,1,\ldots ,N-1} \right) .$ The expansion coefficients of the Fourier polynomial (1) are contained in the set of complex-valued time signal points $\left\{ {c_n } \right\}$. The total signal duration is T, where $T=N\tau$, with N being the total signal length and $\tau$ is the sampling time (dwell time, sampling rate). The inverse of $\tau$ is the bandwidth (BW). The variables $\hbox {exp}\left( {\pm 2\pi imn/N} \right)$ are the undamped sinusoids and cosinusoids $\left( {nm\tau /T=nm/N} \right)$. The continuous time variable t is discretized (digitized) according to $t=n\tau \left( {0\le n\le N-1} \right)$. Insofar as the signal lengths are in the composite form, $N=2^{k}\left( {k=1,2,3,\ldots } \right) ,$ only $N\hbox {log}_2 N$ multiplications are required. This is the basis of the computational efficiency of the FFT algorithm. However, whenever N is non-composite, i.e. any positive integer, the FFT from (1) becomes the discrete Fourier transform (DFT) requiring $N^{2}$ multiplications.

The time signal can be retrieved from the Fourier spectrum through the inverse Fourier transform (IFFT) for $N=2^{k}\left( {k=1,2,3,\ldots } \right)$:

$$\begin{aligned} \hbox {IFFT: } c_n =\frac{1}{N}\sum _{m=0}^{N-1} {F_m \exp (2\pi imn/N)},\quad 0 \le n \le N - 1. \end{aligned}$$

(2)

Insofar as N is non-composite $N\ne 2^{k}\left( {k=1,2,3,\ldots } \right)$, Eq. (2) will be the inverse discrete Fourier transform (IDFT).

Fourier-based processing has no possibility for interpolation, since the total shape spectrum is produced from pre-assigned frequencies whose minimal separation is fixed by the given acquisition time T. Attempts to improve resolution in the FFT inevitably deteriorate Signal-noise ratio (SNR). The reason is that this entails use of a longer T, which means that the physical part of the MRS time signal will have decayed and, thus, encoding would collect mainly noise. Such a worsening of SNR occurs particularly in clinical MR scanners (1.5 and 3 T) [32]. Poor resolution and low SNR are also due to the lack of extrapolation capabilities and to the linearity of the Fourier method. Due to the latter, noise is imported as intact from the time to the frequency domain. Regarding the former, without extrapolation, information is limited to that obtained from $c_0$ up until the final encoded signal point, $\hbox {c}_{N-1}$. It is common practice in the FFT to perform zero-filling of the original set $\left\{ {c_n } \right\} \left( {0\le n\le N-1} \right)$ to e.g. double the time signal length. Whereas this may somewhat enhance the appearance of the generated spectrum (notwithstanding the sinc-type artificial oscillations on the baseline), no new information is provided thereby and, consequently, there is no actual improvement in resolution. Crucially, Fourier-based processing is non-parametric and, thus, only a total shape spectrum can be forthcoming. Post-processing through fitting is often performed thereafter. The latter entails a presupposition about the number and nature of the resonances present in the total shape spectrum. Clearly, however, this guessing procedure can easily be incorrect, with consequent inaccuracies in estimating metabolite concentrations that are clinically most informative [33].

1.2.2 Advanced signal processing for MRS through the fast Padé transform

The fast Padé transform, FPT, is an advanced signal processor which has been shown to be highly suitable for handling MRS time signals [32,33,34,35,36]. Through the unique ratio of two polynomials, $P_K/Q_K$, of degree K in the diagonal form, the spectrum generated by the FPT is a non-linear response function. With the FPT, there is no requirement for a fixed Fourier mesh $2\pi m/T\left( {m=0,1,2,\ldots ,N-1} \right) .$ Consequently, the spectrum can be computed at any sweep frequency $\nu$, where the angular or circular ($\omega$) and linear ($\nu$) frequencies are related by $\omega =2\pi \nu$. Thus, with Padé-processing, the dilemma does not arise whereby attempts to improve resolution require increased T which, in turn, worsens SNR.

In the FPT, an additional degree of freedom is provided by the numerator ($P_K$) and denominator ($Q_K$) polynomials, helping to cancel noise from the Padé spectrum $P_K /Q_K$. Namely, $P_K$ and $Q_K$ through the expansion coefficients of these two polynomials, contain a similar amount of noise from {$c_n$}. This cancellation coheres with the typical circumstances in which for any two observables A and B generated either experimentally (or theoretically through numerical computations with finite precision), there is considerable noise cancellation in the ratio A / B. Overall, the non-linearity of the FPT improves resolution and SNR by suppressing noise [32, 33, 36].

In sharp contrast to Fourier processing, through the FPT, extrapolation beyond the acquisition time T is achieved. This extrapolation capability is based upon the unique polynomial quotient $P_K /Q_K$, extracted directly from the investigated FID.

Crucially, the FPT can reconstruct not only total shape spectra, but also their components through parametric processing (quantification). Thereby, the number of genuine resonances and their spectral parameters, i.e. the fundamental frequencies {$\omega _k$} and the associated amplitudes {$d_k$}, through the set $\left\{ {\omega _k,d_k } \right\} \left( {1\le k\le K} \right)$ inherently present in the specified time signal {$c_n$} ($0 \le n\le N-1$), are reconstructed with a high level of accuracy. It is only from these retrieved parametric data that metabolite concentrations can be correctly computed [33, 36].

The two variants of the FPT

There are two variants of the FPT denoted by the $\hbox {FPT}^{(+)}$ and $\hbox {FPT}^{\left( - \right) }$, that are defined in terms of variables z and $z^{-1}$, respectively. The former converges inside (${\vert }z{\vert }< 1$) and the latter outside (${\vert }z{\vert }> 1$) the unit circle in the complex plane of the harmonic variable z. Moreover, via the Cauchy analytical continuation, the $\hbox {FPT}^{\left( + \right) }$ and $\hbox {FPT}^{\left( - \right) }$ are convergent, as well, in their complementary domains (outside and inside the unit circle, respectively). The $\hbox {FPT}^{\left( - \right) }$ for ${\vert }z{\vert }> 1$ is an accelerator of the already convergent input series given by the Green function in the harmonic variable $z^{-1}$. In contrast, through the Cauchy analytical continuation, the $\hbox {FPT}^{\left( + \right) }$ must force convergence of the input series which diverges inside the unit circle, ${\vert }z{\vert }< 1$ [37].

Concomitant with achieving this more difficult task, the $\hbox {FPT}^{\left( + \right) }$ has a number of advantages for practical applications, as needed for the in vivo MRS setting. Namely, the $\hbox {FPT}^{\left( + \right) }$ is more powerful than the $\hbox {FPT}^{\left( - \right) }$ in handling noisy MRS time signals. Separation of noise from signal is more efficient in the $\hbox {FPT}^{\left( + \right) }$ through which genuine and spurious resonances are strictly partitioned into two opposite regions, inside and outside the unit circle, respectively. In contrast, with the $\hbox {FPT}^{\left( - \right) }$, both these resonance types are mixed together, since they are all located outside the unit circle, ${\vert }z{\vert }> 1$.

Internal cross-validation is provided by the mentioned two variants of the FPT. Whenever complete convergence is attained in the FPT $^{(\pm )}$, it would follow that $\omega _k^+ \approx \omega _k^-$ as well as $d_k^+ \approx d_k^-$, and these parameters can jointly be denoted as $\omega _k \hbox { and }d_k$, respectively. Here, the parameters {$\omega _k, d_k$}, reconstructed by both the $\hbox {FPT}^{\left( + \right) }$ and $\hbox {FPT}^{\left( - \right) }$, are the initially unknown complex frequencies and amplitudes from the input encoded time signal {$c_n$}, which is modeled by the geometric progression:

$$\begin{aligned} c_n =\sum _{k=1}^K {d_k } \hbox {e}^{in\tau \omega _k } \quad (\hbox {Input time signal or FID}). \end{aligned}$$

(3)

Mathematically, integer $K\ge 1$ is the model order, as well as the common degree of the polynomials $P_K$ and $Q_K$ in the spectrum $P_K /Q_K$, whereas in physics it represents the number of resonances. In the present paper, we will from here on, refer exclusively to the $\hbox {FPT}^{\left( + \right) }$, with the understanding that, as in all previous work, in the final analysis, both the $\hbox {FPT}^{\left( + \right) }$ and $\hbox {FPT}^{\left( - \right) }$ are always used for cross-checking.

The response function

The exact response function is given by the infinite-rank Green function $G(z^{-1})$, which is defined as the Maclaurin series:

$$\begin{aligned} G(z^{-1}) =\sum _{n=0}^\infty {c_n } z^{-n},\quad z=\hbox {e}^{i\tau \omega }\,\ (\hbox {Exact Green series}). \end{aligned}$$

(4)

Therein, the time signal points {$c_{n}$} are the expansion coefficients. In realistic situations, however, the total number N of available signal points {$c_{n}$} is finite ($N < \infty$). Thus, a truncated response function is needed. This is provided as the finite-rank Green function given by the Green polynomial $G_N (z^{-1})$:

$$\begin{aligned} G_N (z^{-1}) =\sum _{n=0}^{N-1} {c_n } z^{-n}\,\, (\hbox {Exact Green polynomial}). \end{aligned}$$

(5)

Another terminology can also be used, namely that of discrete time series. Then, the infinite- and finite-rank Green functions can alternatively be called the infinite and finite z-transform [33].

In the $\hbox {FPT}^{(+)}$, the input response function $G_N (z^{-1})$ from (5) is approximated by the causal Green–Padé function $G_K^+ (z)$ as, e.g. the diagonal rational polynomial in the harmonic variable z:

(6)

Here, the term “causal” means that the system must first be perturbed before giving its response. For example, if the excitation started at the initial instant $t_0=0$, the system’s response through the FID, $\{c_n\}(t_n=n \tau ,\tau >0)$, would not appear at any $t<t_0$, which implies that for $c_n=0$ and $n<0$. If such a time signal is used, the resulting frequency response function $G_K^{(+)} (z)$ from (6) would also be causal. Alternatively, $G_K^{(+)} (z)$ can be called the advanced Green–Padé function since it is associated with time evolution of the system along the positive portion of the time axis. Similarly, $G_K^{(-)} (z^{-1})$ from the $\hbox {FPT}^{(-)}$ is called the anti-causal Green–Padé function because of its association with time evolution of the system at the negative portion of the time axis. These associations in $G_K^{(\pm )} (z^{\pm 1})$ are rooted in the harmonic variables $z^{\pm 1}$ which act as operators that propagate the system at positive and negative times, respectively. The two descriptions by the $G_K^{(\pm )} (z^{\pm 1})$ are equivalent by reference to the so-called micro-reversibility of physical processes, as well-known in quantum-mechanical resonance scattering theory. The usual notion of time evolution is propagation in the future, i.e. at positive times. Propagation at negative times is also used as an equivalent notion for theoretical descriptions of time evolution in the past, as per the mentioned micro-reversibility principle.

Via $G_K^{(+)} (z)$, the $\hbox {FPT}^{(+)}$ employs the variable z and converges inside the unit circle (${\vert }z{\vert }< 1$), which is where the exact Green function $G(z^{-1})$ diverges. Then, as noted, via the Cauchy concept of analytical continuation, the $\hbox {FPT}^{(+)}$ induces convergence into the input divergent series from (4). The convergence radii $R_N$ of $G(z^{-1})$ and $R_K^+$ of $G_K^+ (z)$ differ greatly, where the former is exactly zero, $R_N=0$ at $N=\infty$ for ${\vert }z{\vert }< 1$, whereas the latter is strictly non-zero, $R_K^+> 0$, in the same region ${\vert }z{\vert }< 1$. Thereby, the $\hbox {FPT}^{(+) }$ extends the validity of the response function (spectrum) to ${\vert }z{\vert }< 1$, where the input Green series $G(z^{-1})$ does not exist, because of its divergence inside the unit circle.

Note that the $\hbox {FPT}^{(+)}$ from (6) has built-in both interpolation and extrapolation features. Interpolation is present because the sweep frequency $\nu$ in $z=\exp (2\pi i\tau \nu )$ can take any value. Extrapolation, which secures prediction in the mathematical modeling by the $\hbox {FPT}^{(+)}$, is implicit in $G_K^{(+)} (z)$ from (6), since the polynomial reciprocal $1/(\sum _{s=0}^K {q_s^+ } z^{s})$ is an infinite sum, i.e. a series: $1/\sum _{s=0}^K {q_s^+ } z^{s}=\sum _{s=0}^\infty {\eta _s^+ } z^{s}$ where $\eta _s^+ =-\sum _{{s}'=0}^{s-1} {\eta _{{s}'} } q_{s-{s}'}^+ (s\ge 1; q_0^+ \equiv 1, \eta _0^+ =1)$ [33].

Extraction of the expansion coefficients

The expansion coefficients of the polynomials $P_K^+$ and $Q_K^+$ are, respectively, $\{p_r^+\}$ and $\{q_s^+\}$, as per (6). It should be noted that there is no free term $p_0^+$ in the expansion for $G_K^{(+)} (z)$, i.e. $p_0^+ \equiv 0$. Using definition (6), the expansion coefficients $\{p_r^+, q_s^+\}$ of the numerator $P_K^+ (z)$ and denominator $Q_K^+ (z)$ polynomials are extracted uniquely from the time signal points {$c_{n}$} by solving a single system of linear equations. The $\hbox {FPT}^{(+)}$ usually requires a larger number of signal points compared to the $\hbox {FPT}^{(-)}$, i.e. a more over-determined system of linear equations for the polynomial expansion coefficients. This is related to its requirement to induce convergence into the divergent input expansion of $G(z^{-1})$ for ${\vert }z{\vert }< 1$.

Solutions of the characteristic equations

The solutions of the characteristic equations, $P_K^+ (z)=0$ and $Q_K^+ (z)=0,$ have the roots denoted by $z_{k,P}^+$ and $z_{k,Q}^+ (1 \le k \le K)$, respectively. Here, the second subscripts P and Q in the variables $z_{k,P}^+$ and $z_{k,Q}^+$ are introduced to distinguish the roots of the polynomials $P_K^+ (z)$ and $Q_K^+ (z)$, respectively. The Cauchy residues of the spectrum $P_K^+ (z)/Q_K^+ (z)$ taken at the fundamental harmonics $z\equiv z_{k,Q}^+$ are the fundamental amplitudes $d_k^+$. When the roots represent simple poles alone, which occurs for non-degenerate (unequal, i.e. non-coincident) roots of $Q_K^+ (z)$, these amplitudes are given by:

$$\begin{aligned} d_k^+ =\frac{P_K^+ (z_{k,Q}^+ )}{Q_K^{{+}'} (z_{k,Q}^+ )}, \quad Q_K^{{+}'} (z)\equiv \frac{\hbox {d}}{\hbox {d}z}Q_K^+ (z),\quad 1\le k\le K. \end{aligned}$$

(7)

A similar and more involved expression also exists for degenerate (coincident roots), as given in Ref. [33]. From the equivalent canonical representation of spectrum $P_K^+ (z)/Q_K^+ (z)$ [33], it follows that the amplitudes $d_k^+$ are proportional to the pole-zero distance (a metric):

$$\begin{aligned} d_k^+ {{\varvec{\propto }}}\,z_{k,Q}^+ -z_{k,P}^+. \end{aligned}$$

(8)

This Cauchy residue, consequently, reflects the behavior of a line integral of a meromorphic function around a specified pole. The Padé spectrum $P_K^+/Q_K^+$ has its poles as the only singularities and, thus, represents a meromorphic function. Hence, the zeros and poles of $P_K^+/Q_K^+$ are given by the roots of the characteristic equations $P_K^+ (z)=0$ and $Q_K^+ (z)=0$, respectively. It is through these outlined steps that the $\hbox {FPT}^{(+)}$ accomplishes reconstruction of the 2K complex fundamental parameters $\{\omega _{k,Q}^{+},d_k^{+}\} (1 \le k \le K)$ with the earlier notation $\omega _k^+$ relabeled as $\omega _{k,Q}^+$ where $\omega _{k,Q}^+ =\left[ {1/\left. {(i\tau )} \right] } \right. \ln z_{k,Q}^+$.

Modes of the component spectra

The component spectra can be displayed in two different modes. The absorption and dispersion components are mixed together in the “usual” (U) mode of component spectra. Since the phases $\varphi _k^+$ are non-zero, the amplitudes $\{d_k^+\}(1\le k\le K)$ are all complex-valued. When there are numerous overlapping resonances, the “absorption” components often appear as skewed absorptive Lorentzians. When setting the reconstructed phases $\varphi _k^+$ “by hand” to zero, $\varphi _k^+ \equiv 0\left( {1\le k\le K} \right) ,$ we generate the “ersatz” (E) mode of component spectra. Thereby, interference effects have been eradicated and purely absorptive Lorentzians are the result. The “ersatz” and “usual” modes of the component spectrum for the kth resonance are defined as:

$$\begin{aligned} \left( {\frac{P_K^+ (z)}{Q_K^+ (z)}} \right) _k^\mathrm{E}\equiv & {} \frac{\left| {d_k^+ } \right| z}{z-z_{k,Q}^+ } \quad (\hbox {Ersatz component } k), \end{aligned}$$

(9)

$$\begin{aligned} \left( {\frac{P_K^+ (z)}{Q_K^+ (z)}} \right) _k^\mathrm{U}\equiv & {} \frac{d_k^+ z}{z-z_{k,Q}^+ } \quad (\hbox {Usual component } k), \end{aligned}$$

(10)

respectively. When $\varphi _k^+=0$, we can return from (10) to (9), through substitution of $d_k^+ \equiv \vert d_k^+\vert \exp \,(i \varphi _k^+)$ by $\vert d_k^+\vert$ which is the magnitude of the amplitude $d_k^+$ [38].

The peak positions (chemical shift, $\text{ Re }(\nu _k^+)$) in the usual and ersatz mode coincide when the former, i.e. $\hbox {Re}(P_K^+ /Q_K^+ )_k^\mathrm{U}$, is in the absorption mode. On the other hand, when $\hbox {Re}(P_K^+ /Q_K^+ )_k^\mathrm{U}$ is a dispersive component, it will have two lobes. In such a case, the peak position $\hbox {Re}\left( {\nu _k^+ } \right)$ for $\hbox {Re}(P_K^+ /Q_K^+ )_k^\mathrm{E}$ will be located between the two lobes of $\hbox {Re}(P_K^+ /Q_K^+ )_k^\mathrm{U}$, as per visual inspection when juxtaposing the plots for $\hbox {Re}(P_K^+ /Q_K^+ )_k^\mathrm{U}$ and $\hbox {Re}(P_K^+ /Q_K^+ )_k^\mathrm{E}$.

The numerator $(P_K^+)$ and denominator $(Q_K^+)$ polynomials in (9) and (10), have the following explicit expressions, implied by (6):

$$\begin{aligned} P_K^+ (z)=\sum _{r=1}^K {p_r^+ } z^{r},\quad Q_K^+ (z)=\sum _{s=0}^K {q_s^+ } z^{s},\quad p_0^+ \equiv 0. \end{aligned}$$

(11)

The expansion coefficients $\left\{ {q_s^+ } \right\}$ for the polynomial $Q_K^+ (z)$ are extracted by solving the system of linear equations $\sum _{s=0}^K q_s^+ c_{{s'+s}}=0$ deduced from (6). The solutions $\left\{ {q_s^+ } \right\}$ are thereafter refined by singular value decomposition (SVD). Once the set $\left\{ {q_s^+ } \right\}$ becomes available, the expansion coefficients $\left\{ {p_r^+ } \right\}$ in $P_K^+$ are computed from the analytical expression $p_r^+ =\sum _{r'=0}^{K-r} c_{{r'}}q_{r'+r}^{+}$. The free term, $q_0^+$ can be set to e.g. 1 or −1 and this does not affect the spectra or the spectral parameters $\{ {\omega _{k,Q}^+,d_k^+ } \}$ ($1\le k \le K$) reconstructed by the $\hbox {FPT}^{\left( + \right) }$. Note that there is a coherence between the two sets $\left\{ {p_r^+ } \right\}$ and $\left\{ {q_s^+ } \right\}$ because the former depends on the latter.

The ersatz component spectra are often useful for visualizing the extent of the overlap of closely lying or hidden resonances, due to lack of interference between the absorption and dispersion modes. Whereas the number of component resonances $(P_K^+ /Q_K^+ )_k^\mathrm{U}$ and $(P_K^+ /Q_K^+ )_k^\mathrm{E}$ is the same, their full widths at half maximae (FWHM) differ. Thus, the peak areas of a given kth component are also different in the usual and ersatz modes. The parameters $\{ {\omega _{k,Q}^+ ,d_k^+ } \}$ with $\varphi _k^+ \ne 0$ from the usual components $(P_K^+ /Q_K^+ )_k^\mathrm{U}$ are to be utilized for computing metabolite concentrations, since interference effects occur for $\varphi _k^+ \ne 0$. The peak areas are affected by $\varphi _k^+ \ne 0$ and so are the metabolite concentrations. Consequently, the usual components $(P_K^+ /Q_K^+ )_k^\mathrm{U}$ together with $\{ {\omega _{k,Q}^+ ,d_k^+ } \}$ should be used to compute the metabolite concentrations, and not the ersatz components $(P_K^+ /Q_K^+ )_k^\mathrm{E}$ with $\{ {\omega _{k,Q}^+ ,\left| {d_k^+ } \right| } \}$.

Computation of the total shape spectra via the Heaviside partial fraction expansions

The total shape spectra from (6) can alternatively be computed by using the Heaviside partial fraction expansions:

$$\begin{aligned} \frac{P_K^+ (z)}{Q_K^+ (z)}=\sum _{k=1}^K {\frac{d_k^+ z}{z-z_{k,Q}^+ }} \quad (\hbox {Heaviside Partial Fractions}). \end{aligned}$$

(12)

The total shape spectra in the ersatz and usual modes are provided using (9) and (10) to construct the Heaviside partial fractions:

$$\begin{aligned} \left( {\frac{P_K^+ (z)}{Q_K^+ (z)}} \right) ^\mathrm{E}\equiv & {} \sum _{k=1}^K {\left( {\frac{P_K^+ (z)}{Q_K^+ (z)}} \right) _k^\mathrm{E} } =\sum _{k=1}^K {\frac{\left| {d_k^+ } \right| z}{z-z_{k,Q}^+ }}\quad (\hbox {Ersatz envelope}), \end{aligned}$$

(13)

$$\begin{aligned} \left( {\frac{P_K^+ (z)}{Q_K^+ (z)}} \right) ^\mathrm{U}\equiv & {} \sum _{k=1}^K {\left( {\frac{P_K^+ (z)}{Q_K^+ (z)}} \right) _k^\mathrm{U} } =\sum _{k=1}^K {\frac{d_k^+ z}{z-z_{k,Q}^+ }} \quad (\hbox {Usual envelope}), \end{aligned}$$

(14)

respectively. The difference between the lhs of (10) and (14) for the kth usual component $(P_K^+ /Q_K^+ )_k^\mathrm{U}$ and the usual envelope $\left( {P_K^+ /Q_K^+ } \right) ^{\mathrm{U}}$ is that the subscript k is omitted in the latter. Similarly, for the ersatz modes in (9) and (13).

Note that the polar structure in the Heaviside partial fraction representation of the $\hbox {FPT}^{(+)}$ in e.g. (14) is, by definition, best suited for functions with polar singularities such as those describing spectra in MRS. It is through the Heaviside representation that the $\hbox {FPT}^{(+)}$ can “swallow” (filter out) all the poles in the sense of properly decomposing the examined total shape spectrum into its true components. By contrast, all the response functions given by a single polynomial are utterly inadequate for describing spectra with peaks, as is clear from the FFT in (1) which is not a polar representation.

Efficiency of the FPT algorithms

The algorithms in the $\hbox {FPT}^{(+)}$ are both efficient and straight-forward. The only required numerical work is to solve a single system of linear equations for the expansion coefficients $\left\{ {q_s^+ } \right\} ,$ and then to root the characteristic polynomials $P_K^+ (z)$ and $Q_K^+ (z)$. The fundamental frequencies $\{ {\omega _{k,Q}^+ } \}$ are reconstructed from the roots $z_{k,Q}^+$ of $Q_K^+ (z)$. The other set of roots $\{z_{k,P}^+ \}$ from the characteristic equation $P_K^+ (z)=0$ serves to separate genuine from spurious resonances, depending on whether $\omega _{k,P}^+ \ne \omega _{k,Q}^+$ or $\omega _{k,P}^+ =\omega _{k,Q}^+$, respectively, where $\omega _{k,P}^+ =\left[ {1/\left. {(i\tau )} \right] } \right. \ln z_{k,P}^+.$ In contrast to e.g. the Hankel–Lanczos singular value decomposition (HLSVD) by which obtaining the amplitudes requires solving yet another system of linear equations, the $\hbox {FPT}^{(+)}$ generates the set $\left\{ {d_k^+ } \right\}$ from the analytical formulae as the Cauchy residues given by (7). The characteristic polynomial rooting is achieved (to machine accuracy) by solving the equivalent eigenvalue problem of the extremely sparse Hessenberg or companion matrix [33].

Signal-noise separation through the $\hbox {FPT}^{\left( + \right) }$ Via the spectral poles and zeros, one obtains the physical parameters of the system which generated the time signals as the system’s response to an external excitation. Since $d_k^+ \propto z_{k,Q}^+ -z_{k,P}^+$, as per (8), there is also a direct relation of the amplitudes with the spectral poles and zeros.

Stability or resilience against external perturbation indicates that the poles are physical. On the other hand, poles that show marked changes with exposure to even minimal perturbation are unstable and, consequently, unphysical. Moreover, unstable poles do not ever converge, as their behavior is noise-like, and just like random fluctuations, they do not ever stabilize due to lack of coherence.

The relation between spectral poles and zeros in the $\hbox {FPT}^{(+)}$ can also distinguish physical from unphysical poles. Unstable structures show pole-zero confluence, i.e. $z_{k,Q}^+ =z_{k,P}^+$ or $z_{k,Q}^+ \approx z_{k,P}^+$ and they are termed Froissart doublets. For stable structures, the poles and zeros are distinct, i.e. $z_{k,Q}^+ \ne z_{k,P}^+$. In keeping with (8), genuine resonances $(z_{k,Q}^+ \ne z_{k,P}^+ )$ have non-zero amplitudes, $d_k^+ \propto z_{k,Q}^+ -z_{k,P}^+ \ne 0$, whereas spurious resonances $(z_{k,Q}^+ =z_{k,P}^+$ or $z_{k,Q}^+ \approx z_{k,P}^+ )$ have zero or close to zero amplitudes $(d_k^+ =0$ or $d_k^+ \approx 0)$. Crucially, genuine and spurious resonances have positive and negative imaginary frequencies, $\hbox {Im}\omega _{k,Q}^+ >0$ and $\hbox {Im}\omega _{k,Q}^+ <0$, respectively. This implies that the exponentials in the reconstructed time signal:

$$\begin{aligned} c_{n}^{ + } \equiv \sum \nolimits _{{k = 1}}^{{K^{\prime }}} d_{k}^{ + } \text {e}^{in\tau \omega _{{k,Q}}^{ + }} = \sum \nolimits _{{k = 1}}^{{K^{\prime }}} d_{k}^{ + } \exp \left( {in\tau {\text {Re}}\omega _{{k,Q}}^{ + } - n\tau {\text {Im}}\omega _{{k,Q}}^{ + } } \right) ,\qquad \end{aligned}$$

(15)

are damped and exploding with increasing time $n\tau$ for genuine and spurious resonances, respectively. Unattenuated harmonics best justify their binning as the unphysical part of the retrieved time signal, $\{c_n^+ \}$. Occasionally, genuine resonances may also have very small amplitudes, $d_k^+ \approx 0$, but nevertheless, it is their feature $\hbox {Im}\omega _{k,Q}^+ >0$ (alongside stability) which allows us to bin them as the physical portion of the recovered FID from (15).

Subsequent to stabilization of the model order K in $P_K^+ /Q_K^+$, i.e. once all the genuine resonances have been reconstructed, further computation of the Padé spectra for a higher degree polynomial, $K+m\,\left( {m=1,2,3,\ldots } \right)$, would only produce more non-physical resonances. These would exhibit pole-zero coincidences ($z_{k,Q}^+ =z_{k,P}^+$) resulting in $d_k^+ =0$ for k = K + m ($m = 1, 2, 3,{\ldots }$). Consequently, pole-zero cancellation occurs, and this yields stabilization of the computed spectra:

$$\begin{aligned} \frac{P_{K+m}^+ (z)}{Q_{K+m}^+ (z)}=\frac{P_K^+ (z)}{Q_K^+ (z)} \left( {m=1,2,3,\ldots } \right) . \end{aligned}$$

(16)

It should be emphasized that the number of physical resonances, i.e. the number of fundamental harmonics K is also considered unknown. In the Padé methodology, this is a parameter whose value must also be reconstructed, rather than guessed. When the reconstructed frequencies and amplitudes have converged, K will then become determined. Simply stated, the model order ${K}'$ (or equivalently, the degree ${K}'$ of the Padé rational polynomial $P_{{K^{\prime }}}^{ + } /Q_{{K^{\prime }}}^{ + })$ for which the reconstructed frequencies and amplitudes have stabilized, will be the exact number K of harmonics contained in the input FID from (3).

To summarize, by gradually augmenting the running common degree ${K}'$ of the Padé polynomials in the diagonal $\hbox {FPT}^{(+)}$ until the set of reconstructed frequencies and amplitudes stabilize, any further augmentation of ${K}'$ would yield only spurious resonances. The latter are recognized by having negative imaginary frequencies, $\hbox {Im}\omega _{k,Q}^+ <0,$ as well as by their instability, together with coincidence of poles and zeros, that, in turn, yield zero or near-zero amplitudes. The Froissart concept is called “Signal-noise separation” (SNS) and has been extensively validated for MRS time signals [39, 40], with analytical confirmation of its mechanism given in Ref. [41]. This stabilization via pole-zero cancellation is a unique feature of the FPT, due to the special form of the rational polynomials for the Padé spectra. Overall, pole-zero coincidences can occur only in the quotients of two polynomials and subsequently, there are pole-zero cancellations. Identifying Froissart doublets through pole-zero confluences with the ensuing stabilization of the Padé spectra is the essential indication that the full information from the input time signal has been exhausted by the $\hbox {FPT}^{\left( + \right) }$. It is in this way that the genuine parameters from all the reconstructed data $\{ {\omega _{k,Q}^+ ,d_k^+ } \}\left( {1\le k\le {K}'} \right)$ in the $\hbox {FPT}^{\left( + \right) }$ can be identified as the true fundamental frequencies and amplitudes $\left\{ {\omega _k ,d_k } \right\} \left( {1\le k\le K} \right)$ from the encoded FID modeled by (3).

In the SNS concept, the phenomenon of coherence is behind the emergence of the system’s stability. The time signal $c_n$ is said to be built from the 2K stable complex pairs $\left\{ {\omega _k ,d_k } \right\}$ into a coherent sum (3) with the interference effect for the non-zero phases $(\varphi _k \ne 0,k\in \left[ {1,K} \right] ).$ Such a stable system can then be viewed as being closed. An attempt to add more configurations beyond the saturation number K would result in a temporary amalgamation of the new components in (3). These, as incoherent, would be unstable and as such, rejected by being assigned zero-valued amplitudes in (3) which correspond to pole-zero cancellation in the Padé spectrum. Quantum-mechanically, an amplitude $d_k$ is the probability of transition from one configuration to another. Thus, if the system assigns $d_k =0$, this means that there is zero chance for incorporating a new configuration (i.e. a new resonance or metabolite) into (3) beyond the K fully occupied states. In this way, by connecting the phenomena of coherence/incoherence with stability/instability and genuine/spurious, the $\hbox {FPT}^{\left( + \right) }$ is able to suppress and indeed eliminate the redundant and extraneous degrees of freedom of the system.

Reliable Padé-based extrapolation of in vivo encoded MRS time signals

The most important feature of mathematical modeling is prediction. The following argument explains how the fast Padé transform achieves prediction based upon extrapolation in the time domain where MRS encoding is made. In (15), the reconstructed time signal $c_n^+$, associated with spectrum $P_{{K^{\prime }}}^{ + } /Q_{{K^{\prime }}}^{ + }$, the running or sweep model order ${K}'$ is the total number of resonances (genuine K and spurious $K_S$, i.e. ${K}'=K+K_S$) extracted from the encoded data $\left\{ {c_n } \right\}$. For its part, $\left\{ {c_n } \right\}$ is modeled by (3) with K resonances in total which, alongside the relation ${K}'=K+K_S$ (i.e. ${K}' > K$), would suggest that over-modeling had occurred in the retrieved time signal $\left\{ {c_n^+ } \right\}$ from (15). However, the $\hbox {FPT}^{\left( + \right) }$ is safe-guarded against this obstacle, since all the $K_S$ extra resonances, being spurious (unstable), have zero-valued amplitudes that, in turn, reduce the sum in (15) to K terms alone, i.e. $c_n^+ \equiv \sum _{k=1}^K d_k^+ \exp ( {in\tau \omega _{k,Q}^+ } )$. Here, upon convergence, as stated, the set $\{ {\omega _{k,Q}^+ ,d_k^+ } \}$ can be denoted by $\left\{ {\omega _k ,d_k } \right\}$, so that $c_n^+ \equiv \sum _{k=1}^K d_k \exp \,\left( {in\tau \omega _k } \right)$. The latter exactly reconstructs the input data $c_n$ modeled by (3) and thus $\{c_n^+ \}=\left\{ {c_n } \right\}$ for the first N points $\left( {0\le n\le N-1} \right)$. The additional data points for $n\ge N$ in the full set $\{c_n^+ \}$ relative to $\left\{ {c_n } \right\}$ are the Padé-based extrapolations that would have been available had the encoding continued after $c_{N-1}$ beyond the total acquisition time T, i.e. at times $n\tau >T.$ Thus, the $\hbox {FPT}^{\left( + \right) }$ is a method with predictive power.

1.3 The fast Padé transform for processing in vivo MRS time signals encoded from the brain

The FPT has been used to process MRS time signals encoded in vivo on 1.5 T MR scanners from pediatric brain tumor [38], pediatric cerebral asphyxia [41, 42] and from healthy adult brain [43]. This method has also been applied to MRS time signals encoded in vivo from healthy adult human brain using 4 and 7 T MR scanners [32, 33, 36, 44,45,46]. The fast Padé transform has consistently demonstrated superior resolution of total shape spectra compared to the FFT. The parametric capability of the FPT has been the most noteworthy, especially for spectrally dense chemical shift regions of the brain, where very closely-overlapping resonances including cancer biomarkers, were clearly resolved and quantified [36, 38, 41, 42].

1.3.1 The problem of the giant water residual

Successful strategies should handle major problems arising with data encoded from clinical 1.5 T scanners. One of these problems is the (still) giant residual water resonance. Via a step function with the non-parametric FPT, an information-preserving windowing procedure for suppressing residual water was introduced and shown not to affect the spectral components within the spectral region of interest (SRI) [38]. There were some non-essential effects at the edges outside the SRI when comparing the water residual suppressed and unsuppressed Padé constructions. Full equivalence of the non-parametrically and parametrically generated total shape spectra was confirmed in the FPT and, therefore, we chose the latter in a more recent study [41]. Therein, we used only the components $(P_K^+ /Q_K^+ )_k$ with chemical shifts from the SRI selected to avoid the giant residual water resonance. We computed the parametrically generated envelopes via $P_K^+ /Q_K^+$ employing the Heaviside partial fraction sum (12). It was shown that with the parametric FPT and the given SRI, the water residual suppression problem could be entirely overcome without the need for windowing.

1.3.2 Spectra averaging

Another severe problem which arises with in vivo MRS is the destabilizing effect of marked changes in the sought model order K. To cope with this over-sensitivity to alterations in K, iterative averaging of spectra was introduced. This is an effective strategy for regularizing spectra by taking the arithmetic average of a set of envelopes computed for a string of selected values of K. The generation of an average envelope can be viewed as a counterpart to “signal averaging” which is routinely carried out by averaging some 200 encoded FIDs in the time domain, to improve SNR. However, there is more to spectra averaging because it can be iterated any number of times due to very efficient reconstructions in the fast Padé transform. No such advantage is afforded in signal averaging where repetition would only extend the patient’s examination time which is already one of the main concerns in encoding FIDs by MRS. Notably, for different values of model order K, numerous and relatively large noise-like spikes appear in the spectra. Using a sequence of values of the model order K, the 1st average envelope (arithmetic average) was produced. This complex average envelope was inverted via the IFFT or IDFT (depending on whether N is $2^{k}\left( {k=0,1,2,\ldots } \right)$ or not) to generate the 1st reconstructed FID. The latter was subjected to the FPT to obtain the next set of envelopes for the same sequence of values of K as considered in the previous iteration. This new sequence of envelopes was averaged and the procedure was repeated until the prescribed accuracy of the reconstructed spectral parameters was attained. With each iteration, there was progressively greater suppression of spurious spectral structures. In addition to the total shape spectra, all four Padé-reconstructed spectral parameters for each genuine resonance also showed progressively diminished fluctuations with consecutive iterations. The values of the spectral parameters were fully stabilized to the minimal level of variance consistent with stochasticity in encoded FIDs, indicating that convergence had indeed been achieved [41].

Spectra averaging (performed just once or iteratively), as a general concept, is an objective and unbiased way of mitigating the effect of redundancy and unphysical degrees of freedom from the reconstructions. This nuisance is a non-coherent part of the extracted information, present both in the input encoded FIDs (noise as well as other uncertainties) and in the reconstructed data (computational round-off errors, unstable recovered resonances, etc.). Since most such errors are random, their appearance is manifested in spectral envelopes through stochastic spikes roaming around from one chemical shift band to another. It is through spectra averaging of envelopes for a number of values of K that various errors can be corrected in the stabilized quotient $P_K^+/Q_K^+$. We are then talking here about the powerful concept of error self-correction through the unique coupling of averaging of Padé spectra and the rational function response of the examined system to external perturbations.

From these most recent in vivo studies on the brain [38, 41, 42], we concluded that this multi-faceted Padé-based strategy could have other clinical applications. Among these are areas of cancer diagnostics where in vivo MRS would be of greatest added value. We have highlighted early ovarian cancer diagnostics, where for nearly two decades, as noted, the need for an effective in vivo MRS-based screening method has been underscored [30, 31]. We proceed now to succinctly review the course of applications of in vivo MRS to the ovary, starting with conventional Fourier processing and then with the fast Padé transform.

1.4 MRS of the ovary: results to date

1.4.1 In vivo encoded MRS time signals using conventional Fourier-based processing

There is a fairly small number of published studies applying conventional FFT-based in vivo MRS to the ovary [47,48,49,50,51,52,53,54,55,56,57,58,59,60]. In Ref. [61], we performed a systematic review and meta-analysis of these studies. Altogether there were 134 malignant ovarian lesions, 114 benign ovarian lesions and 3 lesions of the ovary that were classified as “borderline”. Encoding was performed via clinical (1.5 or 3 T) MR scanners. For all these in vivo studies, the $^{1}\hbox {H}$ MRS time signals were processed by the FFT, with fitting-based post-processing in Refs. [49, 55, 56, 58, 59].

The main peaks resolved were: lipid (Lip) at 1.3 ppm (parts per million), a lactate (Lac) doublet also at around 1.3 ppm (J-modulated and appearing as inverted for TE = 136 ms), creatine (Cr) at 3.0 ppm, choline (Cho) at 3.2 ppm or total Cho between 3.14 and 3.34 ppm and Lip at 5.2 ppm. A peak whose resonance frequency was between 2.0 and 2.1 ppm was found in some of the studies [52,53,54,55, 59]. According to in vitro analysis, this peak was comprised of N-acetyl aspartate (NAA) as well as N-acetyl groups from glycoproteins and/or glycolipids [54]. In this context, the NAA component may be related to a molecular water pump [52]. Choline is a marker of membrane damage, cellular proliferation and cell density, reflecting phospholipid metabolism of cell membranes. Creatine is a marker of energy metabolism. Anaerobic glycolysis is indicated by Lac. The metabolic information was mainly described qualitatively (presence or absence of a given peak), and these are the data that were pooled for meta-analysis. In the literature [48, 49, 55, 56, 58, 59], diverse procedures were performed aimed at quantifying Cho, such that those data could not be pooled in Ref. [61]. Metabolite concentration ratios of Cho to Cr were reported in two of the studies, with the [Cho]/[Cr] ratio significantly higher in cancerous ovarian lesions than in the benign ovaries [55, 59].

Only two metabolites, Cho at 3.2 ppm and Lac at 1.3 ppm, were more often detected with statistical significance in cancerous compared to benign ovarian lesions. However, based on the detection of Cho alone, some 50 benign ovarian lesions would be erroneously classified as cancerous, i.e. false positive results, indicated by a positive predictive value (PPV) of 66%. Some 20 malignant ovarian lesions would be incorrectly considered benign based upon lack of detected Cho, i.e. false negative results, with a negative predictive value (NPV) of 57.4%. Lactate alone provided better PPV and NPV with a statistically more significant logistic regression model, but data were available for only 25% of the patients. Logistic regression models with adjustment for age and magnetic field strength, $B_{0}$, generated a stronger model for Cho with better NPV than did the unadjusted model, but with markedly fewer patients included. The adjusted model with both Lac and Cho, also with a total of fifty patients, provided the best PPV, NPV and overall accuracy. Even with this strongest logistic regression model, four of twenty-six patients with benign ovarian lesions were wrongly predicted to have ovarian cancer, and four of twenty-four patients with ovarian cancer were incorrectly predicted to have benign lesions.

The overall conclusion from this meta-analysis [61] was that, although some insights were provided with conventional Fourier-based processing, in vivo MRS still did not adequately distinguish cancerous from benign ovarian lesions. Motivations for application of the FPT were provided by the quite extensive data using in vitro MRS [30, 54, 62,63,64,65,66,67,68,69]. Further, as thoroughly reviewed in Ref. [61], there are many MR-observable compounds that do distinguish malignant versus benign ovarian lesions, and such metabolites need to be reliably quantified.

1.4.2 Studies applying the fast Padé transform to MRS time signals to the ovary

Proof of principle studies

The $\hbox {FPT}^{\left( - \right) }$ was first applied to synthesized noiseless time signals [70, 71] associated with MRS data for benign and cancerous ovarian cyst fluid analyzed in Ref. [64]. All the spectral parameters taken from Ref. [64] were correctly reconstructed by the $\hbox {FPT}^{\left( - \right) }$, for each of the input twelve true metabolites. This included the closely-lying resonances, isoleucine (Iso) at 1.023 ppm and valine (Val) at 1.042 ppm. The metabolite concentrations were accurately computed with only 64 signal points (N / 16) of the full time signal $N = 1024$ [70, 71]. At longer signal lengths, these results remained stable. In comparison, at the partial signal length $N_\mathrm{P} = 64$, the FFT provided exceedingly rough spectra that were entirely uninterpretable. In order to resolve all 12 resonances, some 8192 signal points were needed, but even then, several of the peak heights were incorrect in the FFT. To provide converged absorption total shape spectra for the noiseless data corresponding to benign and malignant ovary, the FFT needed a formidable 32768 signal points [64, 70, 71]. The superior resolving power of the $\hbox {FPT}^{\left( - \right) }$ for handling time signals from the ovary was clearly demonstrated thereby [70, 71].

The $\hbox {FPT}^{\left( - \right) }$ was applied thereafter to simulated noisy MRS time signals associated with ovarian cancer, for which the noise level was $\sigma = 0.01156$ RMS, where RMS is the root-mean-square of the noise-free time signal [72,73,74]. Convergence required 128 signal points ($N/8, N = 1024$) with this noise level in order to exactly reconstruct all the spectral parameters for the input twelve physical resonances. The 52 spurious resonances were identified by their pole-zero confluences as well as the associated zero-valued amplitudes [72].

At higher levels of noise ($\sigma = 0.1156$ RMS, $\sigma = 0.1296$ RMS and $\sigma = 0.2890$ RMS), the pole-zero coincidence in the $\hbox {FPT}^{\left( - \right) }$ was not complete in all instances and near-zero amplitudes were found for some of the spurious resonances [73, 74]. Noting that the concentrations of genuine metabolites could be low, their peak heights would be very small. Consequently, the question was how to distinguish non-physical from genuine resonances with full confidence? This question was critical when proceeding from synthesized FIDs to encoded time signals for which the number of resonances and their parameters are unknown prior to spectral analysis. By varying the partial signal length $N_\mathrm{P}$ and/or also by adding yet more noise, it was shown that a set of resonances, even those with very small amplitudes, was detected by their stability. These were classified as genuine resonances. Moreover, another set of resonances was identified by the $\hbox {FPT}^{\left( - \right) }$ as unstable even with a slight change in partial signal length or noise level $\sigma$. These were binned as spurious resonances that are, as such, discarded. In this way, all the true metabolic information was kept in the denoised spectrum [73, 74]. Detailed analysis [75] was next performed on the synthesized noise-corrupted benign ovarian cyst time signals reminiscent of those encoded in the study from Ref. [64]. At short total signal lengths, both FPT variants, the $\hbox {FPT}^{\left( + \right) }\hbox {and } \hbox { FPT}^{\left( - \right) },$ identified all the genuine resonances, and computed the correct metabolite concentrations. The $\hbox {FPT}^{\left( + \right) }$ provided the most effective Signal-noise separation, SNS, because in the complex z-plane the genuine and spurious resonances were separated inside and outside the unit circle, respectively. Equivalently, also in the complex frequency plane, $\hbox {FPT}^{\left( + \right) }$ distinctly separates genuine and spurious resonances, since they have $\hbox {Im}\omega _{k,Q}^+ >0$ and $\hbox {Im}\omega _{k,Q}^+ <0$, respectively. The pole-zero coincidence of spurious resonances remained complete in the $\hbox {FPT}^{\left( + \right) }$, with a denoised spectrum produced automatically for these simulated MRS data from benign ovary. These latter findings suggested that the $\hbox {FPT}^{\left( + \right) }$ could be especially helpful for processing in vivo MRS time signals encoded from the ovary.

Besides applying the FPT to theoretically synthesized FID data similar to time signals encoded from the ovary, several proof-of-principle studies have also been performed using data from other tissues. Padé-processing of MRS time signals provides quantitative information for many metabolites from cancerous, benign and/or normal brain, breast and prostate, as shown in studies on simulated time signals [34,35,36, 76,77,78,79,80,81,82,83,84]. A proof-of-principle MRS study was also performed on the signals encoded from the standard General Electric (GE) phantom head [85]. The convergence process was examined in detail via “parameter averaging”, confirming the accuracy and stability of the Padé-reconstructed spectral parameters, even for those resonances that were very closely overlapping. These parameters were the complex-valued fundamental frequencies $\{ {\omega _{k,Q}^+ } \}$ and the associated amplitudes $\left\{ {d_k^+ } \right\}$.

The first study applying the fast Padé transform to in vivo MRS time signals encoded from the ovary

The first study [61] has just been completed applying the fast Padé transform to MRS time signals encoded in vivo from a borderline serous cystic ovarian tumor on a 3 T MR scanner. At a relatively short partial signal length of $N_{\mathrm{P}}= 800$, the FPT-generated total shape spectrum was found to be better resolved compared to that produced via Fourier processing. This corroborates the mentioned benchmarking studies [70,71,72,73,74,75] on synthesized MRS time signals from the ovary, in which the high resolution capabilities of the FPT were demonstrated.

The spectra averaging procedure was further confirmed to be an effective way to stabilize shape estimation in face of a marked sensitivity to alteration in model order K. The total shape spectrum generated as the real part of the complex average of eleven usual envelopes was seen to be informative, and was very dense, due to the encoding at a short echo time (TE) of 30 ms, such that many short-lived metabolites had not yet decayed. Further, the complex average envelope was inverted to generate a new time signal. Using the latter FID, subsequent parametric analysis through the $\hbox {FPT}^{\left( + \right) }$ reconstructed dense component spectra in the usual, U, and ersatz, E, mode. Recall that in the former, the absorption and dispersion components are mixed. However, in the latter, only the absorptive Lorentzian components exist, since the reconstructed phases are set to zero in order to artificially eliminate interference effects. A large number of metabolites, including potential cancer biomarkers, were identified and quantified. These included Iso, Val, Lip, Lac, alanine (Ala), lysine (Lys), NAA, N-Acetyl neuraminic acid (AcNeu), glutamine (Glu), Cho, phosphocholine (PC), myoinositol (m-Ins). Many of these resonances are difficult if not impossible to detect with Fourier plus fitting of in vivo MRS data from the ovary.

1.5 Aim of the present study

In the first study [61] employing the FPT for in vivo MRS time signals encoded from the ovary [54], the reconstructed spectral parameters, the complex frequencies $\{ {\omega _{k,Q}^+ } \}$ and associated complex amplitudes $\left\{ {d_k^+ } \right\}$ were assessed for a single value of K. In the present paper, we aim to examine the convergence of these spectral parameters by scrutinizing Padé quantification at several values of K for these same in vivo MRS time signals encoded from the ovary [54]. Proving stability of reconstructions, i.e. finding minimal variance of the retrieved parameters for varying model order K, will be the most stringent test of accuracy, as was previously shown in our studies of phantom head data [85] and in vivo time signals encoded from the brain [41].

In Ref. [61], we employed a single averaging procedure with a wider increment among the model orders than in our previous studies [41, 42] in which iterative averaging was benchmarked. The former, i.e. a single averaging procedure would be the most practical for clinical implementation. The question then arises as to what is the true added value of the averaging procedure within the FPT, vis-à-vis convergence of the spectral parameters. This question will be scrutinized in the present paper. With the aim of further optimizing the practical implementation of MRS for in vivo encoded time signals, we will also examine the contribution of the extrapolation features of the FPT, i.e. its unique capability as a rational function, to glean salient information beyond the last encoded signal point $c_{N-1}$ $\left( \mathrm{i.e. \, at} \, \, {t>T} \right) .$

2 Methods

2.1 Acquisition of the MRS time signals

The encoded FID data are from a 56 year-old patient with an enlarged left ovary, as detected on TVUS, and who was included in the in vivo MRS study of Ref. [54]. These data were kindly provided to us by our colleagues from the Department of Obstetrics and Gynecology, Radboud University Nijmegen Medical Center in the Netherlands. The MRS time signals were encoded using a 3 T Magnetom Tim Trio, Siemens MR clinical scanner. Each FID contained 1024 data points. The bandwidth was BW = 1200 Hz, and the Larmor frequency was $\nu _{\mathrm{L}} = 127.732$ MHz corresponding to the magnetic field strength $B_{0} = 3$ tesla ($B_{0} = 3 \hbox { T}$). The sampling time was $\tau = 0.833$ ms ($\tau = 1/\hbox {BW} \approx 0.833 \hbox { ms}$).

A point-resolved spectroscopy sequence (PRESS) was used for this single-voxel proton MRS study. The voxel of interest (3 cm $\times$ 3 cm $\times$ 3 cm) was located in the inferior cystic part of the tumor. The repetition time (TR) was 2000 ms. In the study of Kolwicjk and colleagues [54] two echo times, TE = 30 and 136 ms were used. A total of 64 time signals were encoded and subsequently averaged to improve SNR, such that according to the standard terminology within MRS, the number of excitations (NEX) was 64. Herein, we examine only the FID encoded at 30 ms. The giant water peak was partially suppressed through encoding via WET (water suppression through enhanced $T_{1}$ effects) [54]. After the in vivo MRS encoding, the ovarian tumor was surgically removed. Histopathologic analysis revealed a borderline serous cystic lesion [54].

2.2 Reconstructions of the MRS time signals via the FPT

The encoded FID of length $N=1024$ was not corrected for the zero-order phase $\varphi _0$. Using the definition in (6), the expansion coefficients of the polynomials $P_{K}^{ + }$ and $Q_{K}^{ + }$ in the ${\text {FPT}}^{{\left( +\right) }}$ were calculated directly from the time signal {$c_{n}$}. Following this first step of the analysis, the non-parametrically computed total shape spectra were obtained at a chosen set of equidistant real-valued sweep frequencies $\nu$. If the phases $\varphi _{k}^{ + }$ of the reconstructed FID amplitudes $d_{k}^{ + } = \left| {d_{k}^{ + } } \right| {\text {exp}}(i\varphi _{k}^{ + } )$ were all equal to zero, $\varphi _{k}^{ + } = 0 \,(1 \le k \le K)$, the real and imaginary parts of the spectra, i.e. ${\text {Re}}\left( {P_{K}^{ + } /Q_{K}^{ + } } \right)$ and ${\text {Im}}\left( {P_{K}^{ + } /Q_{K}^{ + } } \right)$, would be purely absorptive and dispersive, respectively.

In encoded MRS time signals, the phases $\varphi _k$ of the FID amplitudes $d_{k}$ are most often non-zero ($\varphi _{k} \ne 0$), due to dephasing which takes place during encoding. Consequently, the reconstructed values $\varphi _{k}^{ + }$ are also such that $\varphi _{k}^{ + } \ne 0$. Thus, in both ${\text {Re}}\left( {P_{K}^{ + } /Q_{K}^{ + } } \right)$ and ${\text {Im}}\left( {P_{K}^{ + } /Q_{K}^{ + } } \right)$ absorption and dispersion lineshapes are mixed together. Parametric analysis through the ${\text {FPT}}^{{\left( +\right) }}$ was carried out with reconstruction of the component spectra in the Usual, U, and Ersatz, E, modes, as described in detail in Sect. 1.2.2. Note, that in the Sect. 3 all the envelopes and components will be plotted at 1024 sweep frequencies, irrespective of whether time signal extrapolation is used or not.

2.3 Averaging procedure and extrapolation through the FPT

As reviewed, arithmetic averaging of spectra has been demonstrated to overcome over-sensitivity of spectra (envelopes, components) as well as of the fundamental frequencies and amplitudes to changes in model order K [41, 42]. In Padé-processing of FIDs encoded from the brain, spectra averaging was previously done with three [42] and nine [41] iterations, to benchmark this procedure. In Ref. [61], on the in vivo MRS data encoded from the ovary, a single averaging of envelopes gave sufficiently accurate reconstructions. As done in Refs. [42, 61], prior to averaging, we shall presently also use only the non-parametric ${\text {FPT}}^{{\left( +\right) }}$ to generate a number of envelopes for a range of model orders K.

It was noted that in the Padé rational functions, as per (15), the spurious resonances cancel out with stabilization for systematically and gradually increased polynomial degree $K+m (m=1,2,3,\ldots )$. This is due to pole-zero cancellations that occur because spurious resonances display coincidence or near-coincidence of their poles and zeros [41]. These confluences (Froissart doublets) render spurious resonances highly unstable, particularly for changes in the model order K. Each envelope $P_{{K + m}}^{ + } /Q_{{K + m}}^{ + }(m = 1,2,3,\ldots )$ will show different spuriousness because of the random distributions of spurious poles and zeros in the complex frequency plane.

In Refs. [41, 42], the increment $\Delta K$ for equidistantly augmenting the value of K was equal to unity, $\Delta K = 1$. With this small $\Delta K,$ some 31 envelopes for $K\in$ [375, 415] were generated for iterative averaging [41, 42]. Alternatively, in our study on in vivo MRS time signals encoded from the ovary [61], a larger interval for K was chosen to be spanned, i.e. $K\in$ [575, 625], but with a longer increment, $\Delta K = 5$, in generation of 11 envelopes for spectra averaging. Any two adjacent spikes were likely to be more markedly different from each other in the case with $\Delta K=5$ than for $\Delta K=1$. In the present study, we also use the increment $\Delta K=5$, to produce 11 envelopes for spectra averaging. The complex 11 usual envelopes $\left( {P_K^+ /Q_K^+ } \right) ^{\mathrm{U}}$ will be computed for $K=575,580,\ldots ,625$ in increments of 5 (i.e. $\Delta K=5$) from the FID encoded at TE = 30 ms. Thereafter, we take the arithmetic average of these 11 envelopes, with the result denoted by $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$. The subscript Av denotes Average (Av). The complex average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ is then subjected to the IDFT to generate a new FID, to which the parametric $\hbox {FPT}^{\left( + \right) }$ is applied. This is the concept of spectra averaging used in Ref. [61], where the recovery of the spectral parameters was done at a single value of $K\left( {K=600,\hbox {the middle of the mentioned interval }575 \le K\le 625} \right)$. The current work goes beyond this by extending reconstructions of spectral parameters to several values of K with and without spectra averaging.

In the present study, we directly apply the extrapolation, as well as interpolation capabilities of the FPT. Recall that the encoded FID was of length $N=1024$. The length of the encoded FID determines the number of Fourier grid frequencies $\nu _m^F \equiv m/T\left( {0\le m\le N-1} \right)$ in the FFT. By contrast, in the FPT, the Padé spectrum $P_K^+ /Q_K^+$ can be computed at arbitrary sweep frequency $\nu$ between any two adjacent values of $\nu _m^F$ and for $\nu >\nu _{m_{\mathrm{max}} }$, amounting to interpolation and extrapolation, respectively. Herein, we compute the complex envelopes $\left( {P_K^+ /Q_K^+ } \right) ^{\mathrm{U}}\hbox {at }5N\left( {5120} \right)$ equidistant sweep frequencies $\nu$ for $K\in \left[ {575,625} \right]$ with increment $\Delta K=5.$ Subsequently, the arithmetic average of these spectra gives the complex average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ at the same 5N frequencies $\nu$. Inverting $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ by the IDFT yields the reconstructed FID of length 5N. Prior to quantification, for convenience, this reconstructed FID is truncated to $2N\left( {2048} \right)$. This is how Padé-based extrapolation, as well as interpolation, of the time domain data, are achieved, with extension of the time signal to twice the original T. The FID reconstructed from the complex average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ is quantified by the parametric $\hbox {FPT}^{\left( + \right) }$ for 6 values of model order $K\in \left[ {575,625} \right]$ with an increment of 10. The ensuing 6 group of spectral parameters are compared with their counterparts for $K=575,585,\ldots ,625$ retrieved from the encoded FID supplemented with 2K-1024 signal points with zero amplitudes. In other words, this latter group of 6 sets of spectral parameters do not stem from spectra averaging nor from Padé-based extrapolation. Comparisons of the two groups of these 6 sets of reconstructed spectral parameters (complex frequencies and amplitudes) will allow an assessment of the overall synergistic effect of spectra averaging and Padé-based extrapolation, which is a primary goal of the present investigation.

3 Results

3.1 Averaging of envelopes through the $\mathbf{FPT}^{\left( + \right) }$

Using the FID encoded at TE = 30 ms, the top panel (a) of Fig. 1 shows (within the spectral region of interest, SRI, between 0.75 and 3.75 ppm) the real parts of 11 usual envelopes $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^\mathrm{U}$ generated from the non-parametric $\hbox {FPT}^{\left( + \right) }$ for $K = 575,\, 580,{\ldots },625$, with an increment of 5. Here, $N_{\mathrm{P}} = 2K = 1150,\, 1160, {\ldots }, \,1250$, meaning that each of the partial lengths is longer than the length (1024) of the encoded FID. Recall that this is done by complementing the missing data by adding only 2K-1024 time signal points with zero amplitudes, and not by zero-filling to 2048 as done in the FFT. Note that the green color is used in panel (a). The largest structure is a spike at about 3.4 ppm, and there are numerous other tall spikes interspersed throughout the entire SRI. Panel (b) displays the real part $\hbox {Re}\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ of the complex arithmetic average for 11 envelopes. This average envelope is marked in blue. Recall that the term “chemical shift” refers to dimensionless real frequency, which is interchangeably denoted by $\hbox {Re}\nu _{k,Q}^+$ in the illustrations.

Fig. 1

The real parts of 11 usual complex envelopes, $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^\mathrm{U}$, marked in green, for $K = 575,580,{\ldots },625$, with increment $\Delta K = 5$, in the frequency band [0.75, 3.75] ppm, computed using the in vivo MRS time signal encoded from a borderline serous cystic ovarian lesion [54], with a 3 T MR scanner (a). Many large noise-like spikes are seen. The 11 complex envelopes are averaged and the real part of the result is denoted by $\hbox {Re}\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ (b), where a “clean” spectrum is generated, shown in blue. Metabolite assignments are presented in b, with full names given in the list of abbreviations. Plots from a and b are superimposed in c. (Color online)

The abbreviations for the metabolites are written above the corresponding peaks in Fig. 1b. Therein, the largest structure is in the region 2.0–2.1 ppm, with two peaks that have a deep splitting between them, delineating the taller narrower acNeu (2.06 ppm) resonance from the shorter and wider NAA (2.03 ppm). Many other structures are identified and marked on panel (b). These assignments are based upon Refs. [54, 64, 66, 67]. Proceeding from the right, i.e. from the lowest chemical shift, at about 0.94 ppm, a doublet and triplet of leucine (Leu) can be seen, followed by Iso, Val and glycine (Gly) at $\sim$1.02, 1.04 and 1.20 ppm, respectively. A Lip resonance, and then a threonine (Thr) peak, a second resonance of Lip and a Lac doublet are seen on a prominence centered at about 1.3 ppm. In the valley subsequently, a multiplet of Iso at 1.48 ppm appears and then a Lys multiplet with a prominent peak around 1.52 ppm. An Ala doublet is noted thereafter, with one of the peaks being rather tall. Within the next valley, two small Leu peaks appear, centered at about 1.73 ppm, and then a Lys peak. Within the next protuberance are a small acetic acid (Ace) doublet and then a Lys doublet, adjacent to the right side of the large NAA peak. Several peaks: glutamine (Gln), methionine (Met), Gly and pyruvate (Pyr) are seen to the left of AcNeu. Next, the large Glu and Gln peaks are seen, centered at about 2.45 ppm. Abutting on the left side is a small citrate (Cit) peak, followed by an m-Ins triplet centered at $\sim$ 2.6 ppm, which is then followed by betaine (Bet), Gln, NAA and then Cit. At 3.0 ppm, a small Cr peak is noticed and to its left is phosphocreatine (PCr), followed by a small tyrosine (Tyr) multiplet and a creatinine (Crn) peak. At 3.2 ppm, a Cho peak can be seen, adjacent to which is a smaller PC peak at $\sim$3.22 ppm, followed by a small glycerophosphocholine (GPC) resonance. Next, another Bet peak is observed at about 3.27 ppm, followed by histidine (His), mannose (Mann) and a multiplet of m-Ins. A very small glucose (Glc) peak appears at $\sim$3.54 ppm, after which is another m-Ins multiplet and Gly thereafter.

Panel (c) of Fig. 1 presents a superimposition of the results of panels (a) and (b), wherein there is a sharp contrast between the spurious structures that, in many cases, are much larger than the genuine peaks. Most prominently, the spike at about 3.4 ppm entirely overwhelms the miniscule physical peaks in that chemical shift region. Also, a large spike completely overrides the Glu peak at about 2.45 ppm. Similarly, but on a smaller scale, the Lac doublet centered at about 1.33 ppm is more difficult to identify in the presence of the spike in that region. The value of spectra averaging is patently observed to reduce the spurious content, such that mainly the genuine structures remain, many of which are very small, as seen on the total shape spectrum, i.e. envelope.

Overall, at several frequency bands, the most evident fluctuations in Fig. 1a occur at the level of peak heights when the values of the model order K are changed. This happens since increased values of K produce more spuriousness, whose randomness most noticeably alters the reconstructed amplitudes because of their global character. Such stochasticity is significantly damped by spectra averaging, as evidenced in Fig. 1b. The reason for this improvement is in the noise-like nature of the spurious part of reconstructions. The same argument is also behind improvements in SNR by time signal averaging, where signal and noise scale differently with the number M of the encoded FIDs, thus yielding a factor of $\sqrt{M}$ improved SNR.

3.2 Convergence of spectral parameters reconstructed by the $\mathbf{FPT}^{\left( + \right) }$ from spectra averaging and time signal extrapolation

We now examine convergence on the level of the parameters reconstructed by the $\hbox {FPT}^{\left( + \right) }$. For reference, panel (a) of Fig. 2 displays the average $\hbox {Re}\left\{ {\hbox {FPT}^{+}} \right\} _{\mathrm{Av}}^\mathrm{U}$ of the real parts of 11 usual envelopes $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^\mathrm{U}$ generated by the non-parametric $\hbox {FPT}^{\left( + \right) }$ for $K = 575, 580,{\ldots },625$, with an increment of 5, as was shown in panel (b) of Fig. 1. Next, a reconstructed FID is produced by the IDFT inversion of the complex average envelope with its real part from panel (a). This new FID (Padé-interpolated/extrapolated relative to the encoded time signal) is then processed by the parametric $\hbox {FPT}^{\left( + \right) }$. Panel (b) of Fig. 2 presents the Argand plot as the imaginary, $\hbox {Im}(\nu _{k,Q}^+ )$, versus real, $\hbox {Re}(\nu _{k,Q}^+ )$, frequencies for six sets of complex frequencies. These are for the interval of $K\in$[575, 625] with an increment of 10, color coded as black, green, cyan, red, magenta, and blue, associated with $K=575,585,595,605,615$ and 625, corresponding to the truncated, or partial signal lengths $N_\mathrm{P} = 1150, 1170, 1190, 1210, 1230$ and 1250, respectively. It can be seen that except for extremely few instances, at the chemical shift region at about 3.6 to 3.65 ppm, there is full agreement to the level of stochasticity among the six sets of reconstructed complex frequencies. Concordantly, in panel (c) of Fig. 2, for plot of magnitude $\left| {d_k^+ } \right|$ versus chemical shift, except for minimal variation in the chemical shift regions at about 3.6 to 3.65 ppm, there is full agreement among the six sets of reconstructed magnitudes. Further, the plot of phase $\varphi _k^+$ versus chemical shift is shown in panel (d). Therein, there is also full agreement among the six sets of reconstructed phases, except for a few slight discrepancies between 3.35 and 3.65 ppm.

Fig. 2

The real part of the usual complex average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ (a). The envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ is inverted via the IDFT producing an FID, which is quantified by the $\hbox {FPT}^{(+)}$ for the interval $K\in \left[ {575,625} \right]$ with increment $\Delta K=10$. The six sets of the results of quantification are color-coded, as black, green, cyan, red, magenta and blue, to give the Argand plot of imaginary, $\hbox {Im}\nu _{k,Q}^+$, versus real, $\hbox {Re}\nu _{k,Q}^+$, frequencies (b), magnitude $\left| {d_k^+ } \right|$ versus chemical shift (c), and phase $\varphi _k^+$ versus chemical shift (d). (Color online)

3.3 The joint impact of averaging and extrapolation on convergence of spectral parameters

Further, we proceed to reconstruct the spectral parameters by directly using the encoded FID to which the $\hbox {FPT}^{(+) }$ is applied for $K \in$ [575, 625] with an increment of 5. No averaging is performed. Moreover, no interpolation nor extrapolation by the Padé rational function is carried out, since the encoded 1024 FID data points are used with additional 2K-1024 zeros.

Fig. 3

The real parts of 11 usual complex envelopes, $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^\mathrm{U}$, for $K=575,580,\ldots ,625$, with increment $\Delta K$= 5 (a). Six of the FIDs (containing 2 parts: 1024 encoded time signal points and 2K-1024 zeros, $2K \ge 1150$) used in (a) are quantified by the $\hbox {FPT}^{(+)}$ for the interval $K\in \left[ {575,625} \right]$ with increment $\Delta K$= 10, also color-coded as black, green, cyan, red, magenta and blue, to give the Argand plot of frequencies $\hbox {Im}\nu _{k,Q}^+$ versus $\hbox {Re}\nu _{k,Q}^+$ (b), magnitude $\left| {d_k^+ } \right|$ versus chemical shift (c), and phase $\varphi _k^+$  versus chemical shift (d). (Color online)

Figure 3a shows the real parts of 11 usual envelopes, $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^{\mathrm{U}}$, generated from the non-parametric $\hbox {FPT}^{\left( + \right) }$ for $K = 575, 580,{\ldots },625$, with an increment of 5, as in Fig. 1a. Panel (b) of Fig. 3 presents the Argand plot as the imaginary, $\hbox {Im}(\nu _{k,Q}^+ )$, versus real, $\hbox {Re}(\nu _{k,Q}^+ )$, frequencies for six sets of complex frequencies reconstructed by the parametric $\hbox {FPT}^{\left( + \right) }$ from the same 6 FIDs used for every second envelope in panel (a). These 6 FIDs differ from each other only in the number (2K-1024) of added zeros, whereas the first 1024 time signal points were encoded. The frequencies from panel (b) are for the interval of $K \in$ [575, 625] with an increment of 10 and, as previously, color coded via black, green, cyan, red, magenta, and blue. Substantial discrepancy is observed throughout the chemical shift region. As such, each reconstructed imaginary frequency, $\hbox {Im}(\nu _{k,Q}^+ )$, at a given chemical shift, $\hbox {Re}(\nu _{k,Q}^+ )$, can very often be distinguished. This is most notable at about 3.6 ppm, where the black, green, cyan, red, magenta and blue circles are totally distinct. Although not quite as pronounced as for the Argand plot, in panel (c) of Fig. 3, for the plot of magnitude $\left| {d_k^+ } \right|$ versus chemical shift, there is also notable spread of the different colors of symbols for a given chemical shift at a number of sites along the displayed SRI. In the chemical shift region around 3.6 ppm, as was observed for the Argand plot, each of the individually reconstructed magnitudes is distinct. The plot of phase $\varphi _k^+$ versus chemical shift in panel (d) also shows marked discrepancies among the six sets of reconstructed phases throughout the entire SRI.

Fig. 4

The real parts of 11 usual complex envelopes, $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^\mathrm{U},$ for $K=575,580,\ldots ,625$, with increment $\Delta K$= 5 (a). Six of the FIDs used in (a) are quantified by the $\hbox {FPT}^{(+)}$ for the interval $K \in [ {575,625}]$ with increment $\Delta K=10$, for the Argand plot of frequencies $\hbox {Im}\nu _{k,Q}^+$ versus $\hbox {Re}\nu _{k,Q}^+$ (b). The real part of the complex usual average envelope $\left\{ {\hbox {FPT}^{( +)}} \right\} _{\mathrm{Av}}^\mathrm{U}$ (c) as in Fig. 1a. The envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ is inverted via the IDFT producing an FID which is quantified by the FPT$^{(+)}$ for the interval $K \in [ {575,625}]$ with increment $\Delta K=10$ to yield the Argand plot of frequencies $\hbox {Im}\nu _{k,Q}^+$ versus $\hbox {Re}\nu _{k,Q}^+$ (d). (Color online)

Fig. 5

The real parts of 11 usual complex envelopes, $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^\mathrm{U}$, for $K=575,580,\ldots ,625$, with increment $\Delta K = 5$ (a). Six of the FIDs used in (a) are quantified by the $\hbox {FPT}^{(+)}$ for the interval $K\in [ {575,625} ]$, with increment $\Delta K=10$ to give the magnitude $\left| {d_k^+ } \right|$ versus chemical shift (b). The real part of the complex usual average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ (c) as in Fig. 1a. The envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ is inverted via the IDFT producing an FID, which is quantified by the $\hbox {FPT}^{(+)}$ for the interval $K \in [575, 625]$ with increment $\Delta K=10$ to yield magnitude $\left| {d_k^+ } \right|$ versus chemical shift (d). (Color online)

Fig. 6

The real parts of 11 usual complex envelopes, $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^\mathrm{U}$, for $K=575,580,\ldots ,625$, with increment $\Delta K=5$ (a). Six of the FIDs used in (a) are quantified by the FPT$^{(+)}$ for the interval $K\in [ {575,625}]$, with increment $\Delta K$= 10, to yield the phase $\varphi _k^+$ versus chemical shift (b). The real part of the complex usual average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ (c) as in Fig. 1a. The complex average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$with $\hbox {Re}\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ from (c) is inverted via the IDFT producing an FID, which is quantified by the $\hbox {FPT}^{(+)}$ for the interval $K \in$ [575, 625] with increment $\Delta K=10$ to yield phase $\varphi _k^+$ versus chemical shift (d). (Color online)

Figures 4, 5 and 6 provide a more direct assessment of the coupling effect of averaging and extrapolation on the convergence of spectral parameters reconstructed by the $\hbox {FPT}^{\left( + \right) }$. As in Fig. 1a, the top panel (a) of Figs. 4, 5 and 6 show the real parts of 11 usual envelopes $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^{\mathrm{U}}$ generated from the non-parametric $\hbox {FPT}^{\left( + \right) }$ for $K = 575, 580,{\ldots },625$, with an increment of 5. Panel (c) of Figs. 4, 5 and 6 display the average envelope $\hbox {Re}\left\{ {\hbox {FPT}^{+}} \right\} _{\mathrm{Av}}^\mathrm{U}$  of the real parts of 11 usual envelopes $\hbox {Re}\left( {P_K^+ /Q_K^+ } \right) ^{\mathrm{U}}$ generated from the non-parametric $\hbox {FPT}^{\left( + \right) }$ for $K=575,580,\ldots ,625$, with an increment of 5, just as was shown in Fig. 1b. Panel (b) of Figs. 4, 5 and 6 present the respective 6 sets of spectral parameters reconstructed by the parametric $\hbox {FPT}^{\left( + \right) }$ from the same 6 FIDs used for every second envelope in panel (a). Panel (d) of Figs. 4, 5 and 6 show the six sets of spectral parameters reconstructed by the $\hbox {FPT}^{\left( + \right) }$ from the IDFT-generated FID given by inverting the complex average envelopes $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$, whose real parts $\hbox {Re}\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$ are from panels (c). In panels (b) and (d), as previously, the six sets are for the interval of $K\left[ {575,625} \right]$ with an increment of 10, and color coded as black, green, cyan, red, magenta and blue.

It should be re-emphasized that the difference between the two types of reconstructions within the $\hbox {FPT}^{\left( + \right) }$ is in the following. Panel (b) on Figs. 4, 5 and 6 use the encoded 1024 FID data points supplemented by 2K-1024 $(K \ge 575)$ zeros with no Padé extrapolation and no spectra averaging. By contrast, panel (d) on Figs. 4, 5 and 6 include the Padé extrapolation and spectra averaging. This is achieved in two steps. First, the 11 complex envelopes for $K\in [575,625]$ with increment of 5 are averaged. Prior to that, each of these 11 envelopes use the encoded 1024 FID data points augmented by 2K-1024$\ \hbox {zeros}\, ({K \ge 575} )$. Second, the obtained complex average envelope $\left\{ {\hbox {FPT}^{+}} \right\} _{\mathrm{Av}}^\mathrm{U}$ computed at ${N}'=5N=5 \hbox {x}1024\,(5120)$ equidistant sweep frequencies $\nu _m (\nu _m \equiv m/{T}', 0\le m \le N'-1, T'=N'\tau )$ is inverted by the IDFT to generate the reconstructed FID of length 5N, subsequently truncated to $2N\left( {2048} \right) ,$ to which the parametric $\hbox {FPT}^{\left( + \right) }$ is applied for producing the spectral parameters shown in Fig. 2, as well as in panel (d) of Figs. 4, 5 and 6. Note that adding 2K-1024 time signal points as pertinent to both panels (b) and (d) of Figs. 4, 5 and 6 is, of course, not a physical extrapolation based on the encoded 1024 FID data points. However, the Padé-based extrapolation, which is associated with panel (d) of Figs. 4, 5 and 6, is a physically motivated extension of the originally encoded FID length 1024 to 2048. This is the case because the average envelope computed at any ${N}'>N$ sweep frequencies $\nu$ is itself a Padé extrapolation in the frequency domain, compared to the FFT spectrum which can be computed only at 1024 Fourier grid frequencies from the encoded $N=1024$ FID data points. Further, when the said Padé complex average envelope at 5N sweep frequencies is inverted by the IDFT, the ensuing reconstructed FID is also of length 5N, which is subsequently truncated to $2N=2048$, for convenience in quantification. Here, the extra 1024 time signal points stem from the Padé rational polynomials, instead of just zero-padding in the FFT. This Padé procedure is a physical extrapolation. The reconstructed FID of length 2048 contains the full informational content from the encoded 1024 FID points. The additional 1024 time signal points are the Padé-predicted data that would have been available had the encoding continued beyond the last originally encoded point, $c_{1023}$. The extrapolated 1024 FID data points beyond $c_{1023}$ are genuine, i.e. physical because they are based upon the encoded 1024 time domain data. Once the reconstructed FID of length $N=2048$ becomes available, its further truncation is performed to the 6 partial signal lengths $N_{\mathrm{P}}\in \left[ {1150,1250} \right]$ with the 6 model orders $K \in \left[ {575,625} \right]$. In panels (b) and (d) of Fig. 4 the Argand plots as the imaginary, $\hbox {Im}(\nu _{k,Q}^+ ),$ versus real, $\hbox {Re}(\nu _{k,Q}^+ ),$ frequencies, are compared. In panel (b) the substantial spread among the six sets of complex frequencies generated from individual envelopes and without extrapolation is seen to be almost completely eliminated in panel (d). In the latter, there is nearly full agreement among the six sets of reconstructed complex frequencies, with very minimal deviations at the chemical shift region of 3.6 to 3.65 ppm. These minimal deviations in the latter chemical shift regions of panel (d) are practically inconspicuous when viewed alongside those same regions in panel (b).

The spread along the magnitude plot in panel (b) of Fig. 5 is also apparent, albeit not quite as marked as was the case for the Argand plot. The improvement provided by averaging plus extrapolation for the magnitudes is noteworthy in panel (d), where the six sets of magnitudes are practically converged except for very minor deviations at about 3.6 to 3.65 ppm.

The synergistic effect of averaging and extrapolation is most prominent for the phase plots seen in Fig. 6. Here, in panel (b) there are, in fact, relatively few chemical shift regions where concordance is observed. In sharp contrast, with the exception of a few slight discrepancies between 3.35 and 3.65 ppm, reconstruction by the $\hbox {FPT}^{\left( + \right) }$ together with averaging and extrapolation produced a converged phase plot.

3.4 Convergence of component spectra as reconstructed via the $\mathbf{FPT}^{\left( + \right) }$ with averaging and extrapolation

We now proceed in Figs. 7 and 8 to examine the convergence of the component spectra built from the Padé-reconstructed parameters that were generated through spectra averaging and extrapolation. To aid in visualization, by correlating the reconstructed chemical shifts and the other parameters with the component spectral lineshapes, panels (b), (c) and (d) of Figs. 7 and 8 reiterate panels (b), (c) and (d) of Fig. 2. Namely, these are the Argand plots as imaginary, $\hbox {Im}(\nu _{k,Q}^+ )$, versus real, $\hbox {Re}(\nu _{k,Q}^+ )$, frequencies, on panels (b), magnitude plots as $\left| {d_k^+ } \right|$ versus chemical shift on panels (c) and phase plots as $\varphi _k^+$ versus chemical shift on panels (d). As previously, these, as well as panel (a) on Figs. 7 and 8, are for the six sets of reconstructed parameters in the interval of $K\in \left[ {575,625} \right]$ with an increment of 10, color coded as black, green, cyan, red, magenta and blue.

Fig. 7

The complex average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$is inverted via the IDFT producing an FID, which is quantified by the $\hbox {FPT}^{(+)}$ for the interval $K\in [ {575,625} ]$ with increment $\Delta K=10$, to generate six sets of usual complex component spectra whose real part is shown in (a) color-coded as black, green, cyan, red, magenta and blue. The absorption and dispersion lineshapes are mixed. The Argand plot of frequencies $\hbox {Im}\nu _{k,Q}^+$ versus $\hbox {Re}\nu _{k,Q}^+$ (b), magnitude $\left| {d_k^+ } \right|$ versus chemical shift (c) and phase $\varphi _k^+$ versus chemical shift (d). (Color online)

Fig. 8

The complex average envelope $\left\{ {\hbox {FPT}^{\left( + \right) }} \right\} _{\mathrm{Av}}^\mathrm{U}$is inverted via the IDFT producing an FID, which is quantified by the FPT$^{(+)}$ for the interval $K\in \left[ {575,625} \right]$ with increment $\Delta K=10$, to generate six sets of ersatz complex component spectra whose real part is shown in (a) color-coded as black, green, cyan, red, magenta and blue. Many closely-overlapping positively-oriented resonances are seen, all in the absorption mode. The Argand plot of frequencies $\hbox {Im}\nu _{k,Q}^+$ versus $\hbox {Re}\nu _{k,Q}^+$ (b), magnitude $\left| {d_k^+ } \right|$ versus chemical shift (c) and phase $\varphi _k^+$ versus chemical shift (d). (Color online)

The six sets of usual component spectra, as per (10), are displayed in panel (a) of Fig. 7. It is noteworthy that there is nearly complete concordance among these six sets. Namely, other than very small deviations of red and green in the dense chemical shift region around 3.4 ppm, the usual component spectra drawn in the other colors appear to be completely merged with the spectrum plotted in the blue color. In other words, these six sets of usual component spectra seem to be as if they were almost a single set. The usual component spectra on panel (a) of Fig. 7 displays a very dense admixture of absorption and dispersion spectra; recall that this mixing is due to the amplitudes $\{d_k^+\}$ being complex-valued with non-zero phases $\varphi _k^+ \left( {1\le k\le K} \right)$. In the usual component spectra the NAA peak appears to be larger than that of acNeu. Such a pattern is opposed to the total shape spectra and the average spectrum from Fig. 1. This can be explained by the argument which runs as follows.

In the usual component spectra, the lineshape of acNeu is in the dispersive mode with a narrow positive and a sizable negative lobe. It is the negative lobe which is responsible for the peak height reduction relative to the case if acNeu were a pure absorptive Lorentzian. Further, in the usual total shape spectrum, constructive interference of the two lobes in the dispersive Lorentzian for acNeu increases the peak height of this metabolite as per Fig. 1. On the other hand, in the same usual component spectra from panel (a) of Fig. 7, the lineshape of NAA is a slightly skewed absorptive Lorentzian. This circumstance is behind the occurrence that the peak heights of NAA are very similar in the usual mode for both the total and component shape spectra in Figs. 1 and 7, respectively. The outlined synthesis/analysis of the peak heights of acNeu and NAA in the usual component spectra explains not only their ratio, but also provides the separate interference mechanisms within each of these two metabolites, in the resulting usual total shape spectrum.

In the chemical shift region around 1.3 ppm, two resonances in the dispersion mode are the most prominent. These correspond to lipids, responsible for the protuberance on the total shape spectrum in the region around 1.3 ppm. Large, wide resonances are observed, as well, in the chemical shift region from about 3.45 to 3.75 ppm. They are also in the dispersion mode, with the negative lobes substantially larger than the positive lobes. These resonances correspond to macromolecules and they contribute to the marked oscillations in the baseline in this chemical shift region. It is seen in panel (c) of Fig. 7 that these macromolecules have the largest magnitudes compared to all the other metabolites in the SRI: $\nu \in \left[ {0.75,3.75} \right]$ ppm. Yet, their peak heights in the envelopes from Fig. 1 are quite small relative to e.g. those of acNeu or NAA. This occurs because these macromolecules, as seen in panel (b) of Fig. 7, also have the widest linewidths in the same SRI. Wide linewidths yield considerable attenuation of peak heights that are proportional to the ratios of the magnitudes and the linewidths [61]. This explanation gives the rationale for the contribution of the discussed macromolecules to the envelopes in Fig. 1 by reference to panels (b) and (c) of Fig. 7.

Similarly, in panel (a) of Fig. 8, the six sets of ersatz component spectra, as per (9), show practically full concordance among themselves. As noted, interference effects are eliminated in the ersatz components through the replacement of $\{d_k^+ \}$ by $\left\{ {\left| {d_k^+ } \right| } \right\}$ which amounts to setting all the phases $\{\varphi _k^+ \}$ of $\{d_k^+ \}$ to zero for 1 $\le k \le K$. This yields purely absorptive Lorentzians throughout. The extremely small discrepancies among the lineshapes from the set of 6 ersatz spectra are localized to a few peaks at around 3.3 to 3.4 ppm. These very minimal discrepancies are mainly magenta and green-coded. With the removal of phases via $\{\varphi _k^+ \}=0\left( {1\le k\le K} \right)$, the right negative lobe of acNeu, has merged with the left positive lobe, such that acNeu peak now appears taller than NAA. Similarly, the Lip resonances centered at 1.3 ppm now appear as larger, symmetrical Lorentzians. Overall, the large number of resonances, well over 90 within this SRI, can clearly be distinguished in these ersatz spectra, although many peak heights are very small.

4 Discussion and conclusions

Several striking results emerge in the present paper using in vivo MRS time signals encoded from the ovary. First is the precision with which the reconstructed six sets of spectral parameters from the FID generated through inversion of the complex average envelope plus extrapolation. The exceedingly small variance among these six sets is a finding which is consistent with our previous studies applying the iterative averaging procedure to in vivo MRS time signals encoded from the brain of a pediatric patient with cerebral asphyxia [41], as well as to our earlier investigation on the standard GE phantom head [85]. In both of these papers, we exhaustively examined the convergence of the Padé-reconstructed spectral parameters, finding a very high level of precision, as seen in minimal variance obtained from the analysis of the reconstructed spectral parameters for the genuine resonances from several consecutive values of the model order (or rank) K, after convergence had been reached.

In the current study, the contribution of the averaging procedure together with extrapolation is clearly demonstrated by our detailed comparison of the results for Padé reconstruction using six FIDs with the same model orders K generated without averaging and without extrapolation. For the latter, convergence of the reconstructed spectral parameters was definitely not achieved. On the other hand, it is herein indicated that a single averaging combined with extrapolation appears to be sufficient to attain convergence of spectral parameters coherent with multiple averaging through iterations. This has important implications for practical implementation, expediting the entire Padé-based methodology.

The present results also show that both component spectra and total shape spectra generated from the Padé-reconstructed spectral parameters are trustworthy. Notwithstanding the fact that the spectral parameters are the most vital for quantitative evaluation in MRS, the usual and ersatz component spectra, as well as the total shape spectra are also helpful by way of visual inspection for clinical interpretation.

The current findings further corroborate our previous contention that this Padé-based strategy would render the use of short echo times, TE, to be reliably exploitable for quantification and, in fact, advantageous, in light of the rich spectral information which can be garnered from short-lived resonances. The vital prerequisite, especially when using short TEs, is the unequivocal disentangling of overlapping resonances, a task for which the parametric FPT is fully capable, as seen particularly herein and in our previous study [61] on FIDs encoded in vivo from the ovary at a TE of 30 ms. In contradistinction, it should be noted that the common practice with Fourier plus fitting estimation of MRS time signals from the ovary (as well as other organs), has been to encode at longer TEs in order to avoid “spectral crowding”, i.e. to have simpler spectra with fewer peaks to fit. Notably, only 3 of the 13 studies included in our meta-analysis [61] used a short TE of 30 ms. Obviously this practice has been due to the problems inherent in Fourier-based processing, alongside subsequent fitting of tightly overlapping peaks. Diagnostic dilemmas arise as a consequence. A notable example is the uncertainty due to the overlap of Lip and Lac at short TEs. The general practice has been to encode at a TE of 136 ms, so that an inverted Lac doublet appears due to J-modulation, whereas Lip has already decayed. Even if Lip did not decay, the Lip-Lac overlap would still be split apart at TE = 136 ms, since at that echo time, Lac and Lip are peaked downward and upward, respectively. However, with the high-resolution capability of the FPT, both Lac and Lip doublets, if present, can be clearly distinguished at any TE. This encourages encoding of MRS time signals at short TEs that provide richer clinical information. Of course, given that the main goal of MRS is quantification, the FPT does not stop at reconstruction of the envelopes alone, as the FFT does. Quite the contrary, the FPT completes its signal processing by reconstructing the spectral parameters from which the components are generated. It is the set of component spectra that are critical to disentangling any overlap of resonances, including Lac and Lip. Thus, indeed, it is only after the components become available via quantification that distinct spectral structures in the total shape spectra can be assigned with certainty to actual metabolites, especially for the average envelopes, as seen in Figs. 1b, 2a and 4c, 5c, 6c of the present study.

A number of other clinically important insights can be gleaned from a review of the total shape spectra together with the components and the spectral parameters. Continuing our consideration of the region around 1.3 ppm, the usual and ersatz component spectra clarify the overlap among Lip, Lac as well as Thr and other resonances. The lifting of the baseline in that chemical shift region is due to wide Lip resonances. Through such insights, it may well be possible to resolve the uncertainty as to whether or not the presence of Lip distinguishes benign from cancerous ovarian lesions. It should be noted that in our meta-analysis [61], Lip at 1.3 ppm was more often identified in malignant lesions. However, this difference was not statistically significant. As mentioned, we did find in our meta-analysis [61], that the presence of Lac at 1.3 ppm was a significant predictor of cancerous lesions, but the data from the existing MRS literature were very scant regarding this metabolite. It can be anticipated that Padé-based parametric analysis using a short TE could identify and quantify both Lip and Lac in the chemical shift region around 1.3 ppm. This would help in clarifying the actual diagnostic importance of Lip and Lac for distinguishing benign from malignant ovarian lesions.

The chemical shift region between 3.20 and 3.24 ppm is also of particular diagnostic importance for identifying ovarian cancer. The components of total Cho lie therein in extremely close proximity. Studies comparing human epithelial ovarian carcinoma cell lines with normal or immortalized ovarian epithelial cells, have shown that phosphocholine, PC, is of three to eight-times higher levels in the cancer cells [69]. Other investigators have identified PC as marker of malignant transformation [86]. In the present and in our previous studies [61], for the first time phosphocholine, PC, as well as GPC were both clearly delineated via in vivo MRS time signals encoded from the ovary. Examining the 6 sets of spectral parameters in that chemical shift region, as reconstructed by the FPT from the extrapolated time signal generated by inverting the complex average envelope, it can be seen that there is complete convergence of all the spectral parameters (complex frequencies and amplitudes) for all the genuine resonances in that chemical shift region. This finding holds promise for non-invasively assessing and quantifying PC with full reliability. With such an achievement, improvement in identifying ovarian cancer can be anticipated. Quantification of PC could also help detect progression of ovarian cancer, based upon studies of ovarian tumor cell lines [69].

Another chemical shift region of interest for ovarian cancer diagnostics is between 2.0 and 2.1 ppm. Complete convergence of all spectral parameters of all the physical resonances was also achieved in that region, with Padé-based reconstruction using the extrapolated time signals generated by inverting the complex average envelope. This advance could help clarify the existing dilemmas regarding the presence of NAA from Fourier-based analysis of in vivo MRS of the ovary [52,53,54,55, 59]. Via this Padé-based strategy, the two resonances between 2.0 and 2.1 ppm corresponding to N-acetyl aspartate, NAA, and N-acetyl neuraminic acid, acNeu, were not only clearly distinguished (via a deep splitting between NAA and acNeu even in the total shape spectra), but also quantified with full accuracy. With these cutting-edge advances, clarification can be anticipated as to the true significance of NAA versus acNeu in identifying malignant as opposed to benign ovarian lesions.

Extensive multivariate exploration will be needed to ascertain the metabolite patterns from MRS that best distinguish benign from borderline or clearly cancerous ovary. Having the rich, quantitative metabolic information that can be garnered with full confidence through this multi-faceted Padé-based strategy, we expect improvement compared to the results thus far achieved with conventional Fourier-based processing of in vivo MRS time signals encoded from the ovary. The overall paucity of in vivo MRS studies on the ovary has been attributed to the difficulty of acquiring good quality time signals from this small, moving organ, and then of reliably processing such time signals to glean meaningful diagnostic information [30, 61]. The present study, with the practical advances needed for clinical implementation, strongly motivates further applications of the Padé-designed methodology for in vivo MRS, with a primary focus upon early ovarian cancer detection. It is also anticipated that this strategy will aid in better identifying benign ovarian lesions, thus helping to avoid unnecessary invasive procedures. It should be recalled that due to the harms related to the latter, ovarian cancer screening is not currently recommended for the general population. Overall, when considering screening platforms within MR and more broadly, a key issue is to avoid harm of any kind [87]. Along with ultrasound, MR-based diagnostic modalities hold promise in this regard, considering that they entail no exposure to ionizing radiation. Such exposure is particularly relevant for women at increased risk of ovarian cancer. Notably, exposure to diagnostic radiation may be associated with further elevation in risk for radiation-induced ovarian cancer [88]. In addition to X-ray-based diagnostic imaging of the pelvic region, radiation exposure may occur with treatment of cervical cancer. Whether due to ionizing radiation exposure, hereditary and/or other risk factors, Padé-optimized MRS could be particularly helpful for surveillance of women with an elevated ovarian cancer risk. It is important to reemphasize the underlying motivation of this work from a clinical perspective, namely that early detection of ovarian cancer would contribute substantially to improved survival for women afflicted with this malignancy. Padé-based signal processing of proven validity with its practical advantages for in vivo MRS is poised to be a prime candidate which can contribute to this long-sought goal.