Robust identification of the cancer biomarker phosphocholine through partitioned envelopes in noisy magnetic resonance spectroscopic data by the non-parametric fast Padé transform

Abstract

The recently introduced concept of envelope partitioning for shape estimation within the fast Padé transform (FPT) is presently further explored and solidified. Earlier, noise-free time signals were used and the results were reported for a single model order K . Currently, partitioned envelopes are computed for several values of model orders K by employing noise-corrupted time signals of increasing standard deviations, $\sigma =$ 0.0289, 0.289, 2.89 in units of root-mean-square of the noise-free time signal. Moreover, spectra averaging is exploited to stabilize shape estimation in face of sensitivity to changes in model order. The main goal of this study is to establish the robustness of the non-parametric FPT for reconstructions of partitioned average envelopes computed with noisy time signals. Both the previous and present illustrations concern synthesized time signals typically encountered in single-voxel magnetic resonance spectroscopy (MRS), akin to in vitro encoding from malignant breast tissue. This particular problem area is chosen for a twofold reason: clinical urgency in cancer medicine, and a huge challenge to reliably identify a key cancer biomarker (phosphocholine), completely hidden underneath a dominant peak (phosphoethanolamine), with a separation of mere 0.001 parts per million of chemical shift. Upon successful benchmarking of partitioned envelopes for noisy simulated MRS time signals, the road would be paved for applications of this special shape estimation to the associated data from in vivo encodings. Partitioned envelopes are important since they offer possibilities to peer into the tightly overlapped resonances by splitting their components apart already at the level of sole total shape spectra. Although constrained by non-parametric estimations, they can still qualitatively decompose the regions of higher spectral density. This would enable subsequent focusing on the most critical spectral regions of interest when solving the local quantification problem by parametric estimations. Such a stepwise strategy is expected to be especially beneficial for multi-voxel magnetic resonance spectroscopic imaging (MRSI), where thousands of noisy spectra need to be processed. The FPT-based partitioned envelopes, followed by accurate local spectral analysis in narrow frequency intervals are poised to help MRSI become an efficient diagnostic modality for everyday clinical practice.

1 Introduction

Time signals are ubiquitous in experimental measurements across vastly different interdisciplinary research fields in basic and applied sciences alike. They also appear everywhere in technological and industrial advances. The great majority of dynamic phenomena in all these branches, and in all of nature for that matter, undergoes time evolution during its development. It is precisely time signals that capture the essence of such evolution to house the entire content of the underlying information. However, the measured time signals do not generally deliver this information to the user “on a silver platter” in a ready-made mode to exploit it. Rather, they very often hide it, fold it and otherwise mask this information. There is a myriad of reasons for such intricate convolutions. The complexity of the examined phenomena, external perturbations and unavoidable noise are among the prominent nuisances that compromise the resilience of the system to disturbances. Add to this the appearance of many time signals, and the challenge of analyzing and interpreting them becomes overwhelming. Many time signals appear as tightly packed waveforms that decay with increasing time. Their mathematical function is invariably a linear combination of exponentially damped sinusoids buried in noise. Little or nothing can be deciphered by inspection of such waveforms in the original time domain of measurement. The reason is that most systems usually have many degrees of freedom and modes of internal motions, encapsulated in oscillations, and manifested through the mentioned complex exponentials that all blend together. This is why usually we cannot simply peer into a time signal waveform and read off its content.

The mentioned exponential damping implies instability as the system responds to a perturbation through time signals whose decay is associated with unsteady, transient phenomena. To discern, at least qualitatively, the structure of a given time signal, computations are vital using a mathematical transformation to generate an equivalent representation called the frequency spectrum. This dual representation is possible thanks to the fact that time and frequency are conjugate variables. It is a spectrum which can transcend some of the time signal’s structure, and uncover a part of the hidden information. A spectrum exhibits a number of relatively discernable peaks or resonances, each of which may correspond to one or more decaying oscillatory waveforms of the time signal. As such, the frequency domain through the computed spectra becomes more amenable to direct interpretation of what has actually been measured in the time domain. A total shape spectrum or envelope provides qualitative information. Quantitative information is given by component shape spectra of an envelope. These qualitative and quantitative interpretations are the subject of a separate branch called signal processing. Therein, one encounters the so-called “shape estimators” for envelopes and “parameter estimators” for components.

The fast fourier transform (FFT) is the prime example of a shape estimator. There are many parameter estimators, one of which is the fast Padé transform (FPT). Advantageously, the FPT can also perform shape estimations through non-parametric processing. Standard envelopes are useful only if there are no overlapping peaks. For overlapping peaks, parameter estimation is needed to carry out the spectral analysis (also called quantification) which splits the given envelope into its component shape spectra. This is how traditionally signal processing has been done.

Of late, however, the possibility was explored to decompose an envelope into its components by shape estimation alone [1]. To this end, the non-parametric version of the FPT has been used as a shape estimator to qualitatively extract the components of compound peaks that appear as single isolated resonances in standard envelopes. This was achieved by exploiting the fact the time domain evolution of a generic system is described in the FPT by way of the auto-regressive moving average (ARMA) process. The $z-$transform (or the transfer function) of the ARMA time domain process is a quotient of two polynomials and this is a spectrum (or response function) in the FPT. Prior to inspecting such a complex-valued spectrum $P_{K}/Q_{K}$ of model order K, it is possible to disentangle it and visualize the separate contributions from its two constituents, the moving average (MA) and auto-regression (AR) given by the numerator $P_{K}$ and denominator $Q_{K}$ polynomial, respectively. The main reason for the overwhelming use of ARMA estimations in vastly different fields is that the MA part efficiently models noise, while the AR part provides the spectral poles. It turns out that the so-called “partitioned envelopes”, as the unique combinations of absorptive and dispersive envelopes from the MA and AR processing, are able to visualize the hidden components of compound peaks [1]. The first illustrations of the concept of spectra partitioning was in a medical application by reconstruction of the Padé-based partitioned envelope spectra for a single model order K. As the entries, synthesized noiseless time signals were generated to be reminiscent of the corresponding data, alternatively called a free induction decay (FID), encoded from cancerous breast tissue via magnetic resonance spectroscopy (MRS).

We chose this example for two reasons. First, it is a case of utmost relevance for breast cancer diagnostics, which could be greatly improved by timely identification of a recognized cancer biomarker, phosphocholine (PC). Second, in an envelope spectrum, the PC peak is completely hidden underneath of a nearby dominant resonance, phosphoethanolamine (PE), which is not a cancer biomarker. The partitioned absorptive envelopes were shown to be able of visualizing both PC and PE that in the associated total shape spectrum ${\mathrm {Re}}(P_{K}/Q_{K})$ appear as a single perfectly symmetrical Lorentzian [1]. By construction, the partitioned envelopes are unique. This is also confirmed in numerical computations by proving that the appropriate combinations of the MA and AR pathways yield exactly the result from direct computations of $P_{K}/Q_{K}$ (with no recourse to separate evaluations of moving average and auto-regression). Guided by the initial success, it was deemed necessary to further test the findings reported in Ref. [1] so as to firmly establish the robustness of visualization of hidden resonances by using noise-contaminated time signals, as well. Moreover, it would be important to carry out reconstructions by the non-parametric FPT for a sequence of model orders K, followed by arithmetic averaging. This is the subject of the present work.

2 Advanced signal processing by the fast Padé transform

Effective strategies for the non-invasive detection of cancer biomarkers through molecular imaging are being urgently sought. Magnetic resonance spectroscopy, MRS, and spectroscopic imaging (MRSI) hold particular promise to realize this aim. In line with this overall goal, we investigated the advantageous properties of rational polynomials for handling functions containing resonances. Due to its uniqueness for the power series expansion of the given function, the Padé approximant is the most important of these rational polynomials. Very recently, we explored the heretofore unrealized possibilities of non-parametric analysis with partitioning via the fast Padé transform, FPT, as an expedient first step in processing MRS time signals. The convenience of this initial analysis is that it can be accomplished once the expansion coefficients of the polynomials are generated from the time signal, without polynomial rooting. Our focus in Ref. [1] was upon the breast cancer biomarker phosphocholine, PC, a hidden component, which had not previously been visualized with non-parametric analysis of MRS time signals via the FPT. The present paper continues with this sub-topic in further depth. We begin with salient background for this multi-faceted problem.

It should be recalled that time signals, i.e. free induction decay, FID, curves in MRS and MRSI must be converted into their spectral representation in the frequency domain. In conventional practice this is done by the fast Fourier transform, FFT, which is built-in to practically all magnetic resonance (MR) scanners. The advantages and drawbacks of the FFT for MRS have been exhaustively reviewed in e.g. Refs. [1,2,3,4]. Herein, we note that the Fourier spectrum is conveyed as a single polynomial. The lack of a polar structure is one of the important reasons for the inadequacy of the FFT in processing MRS data, since these are representations of functions with peaks [2, 5]. Thus, we now turn our attention to the FPT starting with a brief review of its basic features.

2.1 Basic features of the fast Padé transform

Through the FPT, a spectrum is generated as a non-linear response function via the unique ratio of two polynomials [2]. In the diagonal form, where K is the model order or polynomial degree, this spectrum is $P_{K}/Q_{K}$. The exact response function is given by the infinite-rank Green function $G(z^{-1})$, which is defined by the Maclaurin series:

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

(1)

where i is imaginary unity ($i=\sqrt{-1})$ and the time signal points $\left\{ c_{n} \right\}$ are the expansion coefficients. Whenever only a finite number $N\,(N < \infty )$ of signal points $\{c_{n}\}$ is available, as is always the case in practice, a truncated response function is needed. This is the finite-rank Green function or the Green polynomial $G_{N} (z^{-1})$:

$$\begin{aligned} G_{N} (z^{-1}) =\sum \limits _{n=0}^{N-1} {c_{n} } z^{-n}\quad (\hbox {exact Green polynomial}). \end{aligned}$$

(2)

A major advantage of rational polynomials, as a quotient of two polynomials, is that the polar representation is automatically built-in by design [5]. For that reason, rational polynomials are the most appropriate for describing functions with peaks, such as MR spectra. Among the rational polynomials, the Padé approximant is the most important because of its uniqueness for the power series expansion of the given function.

There are two variants of the FPT with respect to the complex harmonic variable z. These are defined inside $(\left| z \right| <1)$ and outside $(\left| z \right| >1)$ the unit circle for the causal and anti-causal representations, and are termed ${\mathrm {FPT}}^{(+)}$ and ${\mathrm {FPT}}^{(-)}$, respectively. Both FPT variants are frequency-dependent polynomial quotients extracted from the common exact Green polynomial (2) [2, 4]. Consequently, the single input response function $G_{N} (z^{-1})$ from (2) in the ${\mathrm {FPT}}^{\left( \pm \right) }$ is approximated by the two Green-Padé functions $G_{K}^{\pm } (z^{\pm 1})$, defined by the diagonal rational polynomials in their respective harmonic variables $z^{\pm 1}$:

$$\begin{aligned} \text{ FPT }^{(+)}: G_{N} (z)\approx G_{K}^{(+)} (z)= & {} \frac{P_{K}^{+} (z)}{Q_{K}^{+} (z)}\equiv \frac{\sum \nolimits _{r=1}^K {p_{r}^{+} z^{r}} }{\sum \nolimits _{s=0}^K {q_{s}^{+} z^{s}} } ; \quad \,\, p_{0}^{+} =0, \end{aligned}$$

(3)

$$\begin{aligned} \text{ FPT }^{(-)}: G_{N} (z^{-1})\approx G_{K}^{-} (z^{-1})= & {} \frac{P_{K}^{-} (z^{-1})}{Q_{K}^{-} (z^{-1})}\equiv \frac{\sum \nolimits _{r=0}^K {p_{r}^{-} z^{-r}} }{\sum \nolimits _{s=0}^K {q_{s}^{-} z^{-s}}}. \end{aligned}$$

(4)

Thus, for the same input $G_{N} (z^{-1})$, the two equivalent Green-Padé representations (3) and (4) for spectra $P_{K}^{+} (z)/Q_{K}^{+} (z)$ and $P_{K}^{-} (z^{-1})/Q_{K}^{-} (z^{-1})$ are the mentioned causal and anti-causal response functions in the ${\mathrm {FPT}}^{(+)}$ and ${\mathrm {FPT}}^{(-)}$, respectively. In fact, many more variants of the ${\mathrm {FPT}}^{(\pm )}$ can be generated for the more general case of non-different degrees K and L of the numerator and denominator polynomials $P_{K}^{\pm }$ and $Q_{L}^{\pm }$ according to the definitions:

$$\begin{aligned} G_{N} (z^{-1})\approx G_{K,L}^{\pm } (z^{\pm 1})=\frac{P_{K}^{\pm } (z^{\pm 1})}{Q_{L}^{\pm } (z^{\pm 1})}+O(z^{\pm K\pm L\pm 1}). \end{aligned}$$

(5)

Here, the remainders $O(z^{\pm K\pm L\pm 1})$ are the Maclaurin series in power of $z^{\pm 1}$ beginning with the first terms $z^{\pm K\pm L\pm 1}$. There are two main important features that emanate from (5). First, $O(z^{\pm K\pm L\pm 1})$ are the errors of the approximations $G_{N} (z^{-1})\approx G_{K,L}^{\pm } (z^{\pm 1})$, as evident from the difference of the input $G_{N} (z^{-1})$ and the output (or model) $G_{K,L}^{\pm } (z^{\pm 1})$ functions via:

$$\begin{aligned} \left. \begin{array}{l} \Delta G_{K,L}^{\pm } (z^{\pm 1})\equiv \{G_{N} (z^{-1})\}-\{G_{K,L}^{\pm } (z^{\pm 1})\}=\{O(z^{\pm K\pm L\pm 1})\}\\ \{\hbox {Input}\}-\{\hbox {Outputs}\}=\{\hbox {Errors}\} \end{array}\right\} . \end{aligned}$$

(6)

Second, the errors $\Delta G_{K,L}^{\pm } (z^{\pm 1})$ are the smallest via $O(z^{\pm K\pm L\pm 1})$ for the diagonal ($L=K$) case, $G_{K,K}^{\pm } (z^{\pm 1})\equiv G_{K}^{\pm } (z^{\pm 1})$, and this is why we opt to work with the diagonal forms of the ${\mathrm {FPT}}^{(\pm )}$ from (3) and (4). The expansion coefficients of the numerators $P_{K}^{\pm } (z^{\pm 1})$ and denominators $Q_{K}^{\pm } (z^{\pm 1})$ are $\left\{ p_{r}^{\pm } \right\} \, \mathrm {and} \left\{ q_{s}^{\pm } \right\} ,$ respectively. By solving a single system of linear equations from definitions (3) or (4), these latter coefficients are extracted uniquely from the time signal points $\{c_{n}\}$ for fixed values of $q_{0}^{\pm }$ . The ${\mathrm {FPT}}^{(-)}$ operates with the reciprocal variable 1 / z and, thus, it is an accelerator of convergence of the input slowly converging expansion in powers of $z^{-1}$ [2, 3]. The ${\mathrm {FPT}}^{(+)}$ works with variable z and, hence, performs analytical continuation of the same input development (3), which is in powers of $z^{-1}$. The ${\mathrm {FPT}}^{(+)}$ is algorithmically more difficult, since it must induce convergence into divergent series [2, 6]. Therefore, for convergence, the ${\mathrm {FPT}}^{(+)}$ generally requires more signal points than the ${\mathrm {\, FPT}}^{(-)}$. By definition, the ${\mathrm {FPT}}^{(+)}$ and ${\mathrm {FPT}}^{(-)}$ converge for $\left| z \right| <1$ and $\left| z \right| >1$, respectively. However, by the Cauchy concept of analytical continuation, the ${\mathrm {FPT}}^{(+)}$ and ${\mathrm {FPT}}^{(-)}$ also converge in their complementary domains $\left| z \right| >1$ and $\left| z \right| <1$, respectively. Thus, the ${\mathrm {FPT}}^{(\pm )}$ converge everywhere in the $z^{\pm 1}$ complex planes, with the exceptions of the poles given by the zeros of the denominator polynomials, $Q_{K}^{\pm } (z^{\pm 1})=0$. The ${\mathrm {FPT}}^{(+)}$ and ${\mathrm {FPT}}^{(-)}$ are always employed simultaneously for a fully self-contained cross validation, using different algorithms [2, 4, 7]. It is for this reason of cross-validation within the same type of methodology that we use both variants of the FPT in practice.

Extensive investigations demonstrate that the FPT achieves high resolution [2,3,4,5, 7,8,9,10,11], for which there is a number of reasons. Firstly, unlike the FFT for which there is a sharp cut-off of the time signal at the end of the total acquisition time, T, the FPT can extrapolate beyond T via the unique polynomial quotient $P_{K} /Q_{K}$. Secondly, in the FPT, the fixed equidistant Fourier mesh frequencies $2\pi m/T\, (m\, =\, 0,\, 1,\, 2,\, \ldots \,)$ are not required. Once the Padé polynomials $P_{K}$ and $Q_{K}$ become available, the non-parametric total shape spectrum $P_{K}/Q_{K}$ can be reconstructed by the FPT at any sweep frequency. Thirdly, due to the non-linearity of the FPT, as a rational response function, noise is suppressed and, thus, signal-to-noise ratio (SNR) is enhanced [2,3,4].

The FPT can be applied as both a non-parametric and parametric processor. Regarding the latter, extensive controlled study with and without added noise indicates that exact reconstruction of the input spectral parameters is achieved by the FPT [4, 9,10,11,12,13,14,15,16,17,18,19,20,21,22]. Extremely closely overlapping resonances are fully resolved in that way by the FPT. Studies of in vivo encoded MRS time signals processed by the parametric FPT concordantly indicate that full convergence of spectral parameters is achieved, even for resonances that are entirely overlapping [5, 23,24,25,26,27,28]. Quantification through the FPT amounts to solving exactly the problem of spectral analysis to determine four real parameters per resonance (position, height, width and phase) [4, 24, 29]. Herein, we focus upon the non-parametric capabilities of the FPT, as an important initial step for processing MRS time signals, particularly in the clinical setting, as mentioned, and will be elaborated further on in this paper.

2.2 Non-parametric computation of envelopes through the FPT with and without partitioning

Notwithstanding the well-documented capabilities of the FPT to carry out quantification in the most reliable way, it was nevertheless deemed important to consider whether shape estimation alone by this method could be explored to gain the first hint about the components of the envelopes and spectral density, prior to precise parametrization of the investigated system. The answer to this query is in the affirmative, as was initially demonstrated in Ref. [1] for synthesized noise-free time signals. To this end, a novel notion of the so-named “partitioned envelopes” was introduced to show that the non-parametric FPT can qualitatively decipher the constituents of a given compound resonance. This is remarkable, especially for e.g. a composite resonance which appears as a single absorptive Lorentzian peak in the conventional envelope computed without partitioning. Throughout, partitioning refers to splitting the absorption and dispersion lineshapes into two distinct sections uniquely extracted analytically (from a complex-valued spectrum) prior to numerical computations. The encouraging results from Ref. [1] with noise-free time signals motivate an extension to also encompass noise-corrupted input data. This is what we are set to do in the present investigation.

Once the expansion coefficients $\left\{ p_{r}^{\pm } \right\} \, \mathrm {and}\, \left\{ q_{s}^{\pm } \right\}$ of the polynomials $P_{K}^{\pm }$ and $Q_{K}^{\pm }$ , respectively, are extracted using the intact input time signal $\{c_{n}\}$, non-parametric analyses by the ${\mathrm {FPT}}^{(\pm )}$ are carried out automatically. The non-parametric Padé envelopes $P_{K}^{\pm } /Q_{K}^{\pm }$ are thereby generated. The standard way to do so has been by directly feeding the set of the extracted complex numbers $P_{K}^{\pm } /Q_{K}^{\pm }$ into the computer to generate the spectral absorptions ${\mathrm {Re}}(P_{K}^{\pm } / Q_{K}^{\pm })$ and dispersions ${\mathrm {Im}}(P_{K}^{\pm } / Q_{K}^{\pm })$. Such lineshapes will be purely absorptive and dispersive, respectively, if in a synthesized FID, all the amplitude phases are set to be equal to zero.

Alternatively, these absorption and dispersion lineshapes can also be computed from the explicit expressions for the real and imaginary parts, respectively, of the complex spectra $P_{K}^{\pm } /Q_{K}^{\pm }$. When ${\mathrm {Re}}(P_{K}^{\pm } /Q_{K}^{\pm })$ and ${\mathrm {Im}}(P_{K}^{\pm } / Q_{K}^{\pm })$ are calculated analytically, i.e. taken “by hand” from $P_{K}^{\pm } / Q_{K}^{\pm }$, the partitioned envelopes are obtained. Polynomials $P_{K}^{\pm }$ and $Q_{K}^{\pm }$ can produce only zeros and poles, respectively, in the spectra $P_{K}^{\pm } /Q_{K}^{\pm }$ that are the meromorphic complex functions. A meromorphic function is a function in which the only singularities are its poles. Importantly, each of the spectra $P_{K}^{\pm } / Q_{K}^{\pm }$ is itself already composed of two parts as the $z-$transforms of the moving averages, MA via $P_{K}^{\pm }$ and that of the auto-regressions, AR, through $Q_{K}^{\pm }$. The combined spectra $P_{K}^{\pm } / Q_{K}^{\pm }$ of both AR and MA processes are the auto-regressive moving averages, ARMA. As such, the $z-$transforms of the ARMA processes are equivalent to the ${\mathrm {FPT}}^{(\pm )}$ [2, 4]. Note that the expansion coefficients $\{q_{s}^{\pm }\} \,\,(1\le s\le K)$ of polynomial $Q_{K}^{\pm }$ correspond to the backward and forward prediction coefficients, respectively, in the AR processes.

Use of the $P_{K}^{\pm }$ alone produces spectra corresponding to an “All-zero model”, whereas an “All-pole model” results when spectra are built only with $Q_{K}^{\pm }$. Thus, the reciprocals ${1/Q}_{K}^{\pm }$ would contain only poles and these are determined by the roots of the characteristic or secular equations $Q_{K}^{\pm }=0.$ By considering the quotients $P_{K}^{\pm } / Q_{K}^{\pm }$ as the equivalent products $[P_{K}^{\pm }]\cdot {[1/Q}_{K}^{\pm }]$, the separate contributions to ${\mathrm {Re}}(P_{K}^{\pm } / Q_{K}^{\pm })$ and ${\mathrm {Im}}(P_{K}^{\pm } / Q_{K}^{\pm })$ could be deduced by the various products of the real and imaginary parts of $P_{K}^{\pm }$ (MA) and ${[1/Q}_{K}^{\pm }]$ (AR). And this is how the partitioned envelopes are generated.

2.2.1 Formulae for the partitioned envelopes

From here on in this paper, we will be directly referring to the ${\mathrm {FPT}}^{(+)}$. We now proceed to create the separate products that follow automatically from the formula for the ratio of any two complex numbers $z_{1}$ and $z_{2}$ as $z_{1} /z_{2} =z_{1} z_{2}^{*} /\left| {z_{2} } \right| ^{2}$. Here, the star superscript specifies the standard operation of complex conjugation. With this rule applied to $P_{K}^{+} / Q_{K}^{+}$, we obtain the mentioned separate products in the form of the quantities ${A}_{K}^{+},\, {\, B}_{K}^{+},\, {C}_{K}^{+}$ and ${D}_{K}^{+}$ that represent the partitioned envelopes:

$$\begin{aligned} \frac{P_{K}^{+}}{Q_{K}^{+}}={\mathrm {Re}}\left( P_{K}^{+}/Q_{K}^{+}\right) +\, {i\cdot \mathrm {Im}}\left( P_{K}^{+} / Q_{K}^{+}\right) , \end{aligned}$$

(7)

where,

$$\begin{aligned} {\mathrm {Re}}\left( P_{K}^{+} / Q_{K}^{+}\right) = {A}_{K}^{+}+B_{K}^{+},\ {\mathrm {Im}}\left( P_{K}^{+} / Q_{K}^{+}\right) = C_{K}^{+}+D_{K}^{+}, \end{aligned}$$

(8)

with

$$\begin{aligned} A_{K}^{+}= & {} \left[ \hbox {Re}\left( P_{K}^{+}\right) \right] \left[ \hbox {Re}\left( Q_{K}^{+}\right) \right] /\left| Q_{K}^{+} \right| ^{2}=\left[ \hbox {Re}\left( P_{K}^{+}\right) \right] {\mathrm {Re}}\left( 1/Q_{K}^{+}\right) , \end{aligned}$$

(9)

$$\begin{aligned} B_{K}^{+}= & {} \left[ \hbox {Im}\left( P_{K}^{+}\right) \right] \left[ \hbox {Im}\left( Q_{K}^{+}\right) \right] /\left| Q_{K}^{+} \right| ^{2}= -\left[ \hbox {Im}\left( P_{K}^{+}\right) \right] {\mathrm {Im}}\left( 1/Q_{K}^{+}\right) , \end{aligned}$$

(10)

and

$$\begin{aligned} C_{K}^{+}= & {} -\left[ \hbox {Re}\left( P_{K}^{+}\right) \right] \left[ \hbox {Im}\left( Q_{K}^{+}\right) \right] /\left| Q_{K}^{+} \right| ^{2}= \left[ \hbox {Re}\left( P_{K}^{+}\right) \right] \hbox {Im}\left( 1/Q_{K}^{+}\right) , \end{aligned}$$

(11)

$$\begin{aligned} D_{K}^{+}= & {} \left[ \hbox {Im}\left( P_{K}^{+}\right) \right] \left[ \hbox {Re}\left( Q_{K}^{+}\right) \right] /\left| Q_{K}^{+} \right| ^{2}= \left[ \hbox {Im}\left( P_{K}^{+}\right) \right] \hbox {Re}\left( 1/Q_{K}^{+}\right) . \end{aligned}$$

(12)

Here, we use the property ${\mathrm {Im}}\left( 1/Q_{K}^{+}\right) =-\left[ \mathrm {Im}\left( Q_{K}^{+}\right) \right] /\, \left| Q_{K}^{+} \right| ^{2}$. As in an “All pole model”, the quantity $1/\left| Q_{K}^{+} \right| ^{2}$ in $A_{K}^{+}$, $B_{K}^{+}$, $C_{K}^{+}$ and $D_{K}^{+}$ is the power spectrum. In the Result Section, we shall present both partitioned spectra and their arithmetic average values for a selected sequence of model orders K. The average partitioned spectra will be denoted as ${{\{X}_{K}^{+}\}}_{\mathrm {Av}}$, where $X=A,\, B,\, C,\, D$ and the complete average envelopes shall be labeled by ${{\{A}_{K}^{+}+B_{K}^{+}\}}_{\mathrm {Av}}$ and ${{\{C}_{K}^{+}+D_{K}^{+}\}}_{\mathrm {Av}}$. The numerical results from the latter two quantities are expected to coincide with the corresponding conventional (non-partitioned) average absorption and dispersion lineshapes, ${{\mathrm {Re}}\left( P_{K}^{+} /Q_{K}^{+}\right) }_{\mathrm {Av}}$ and ${\mathrm {Im}}\left( P_{K}^{+} /Q_{K}^{+}\right) _{\mathrm {Av}}$, respectively.

Overall, we see that there is an additional degree of freedom with the non-parametric complex spectrum $P_{K}^{+} /Q_{K}^{+}$ consisting of an alternative way of computing ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$ and ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$. As stated, customarily, these two latter absorption and dispersion spectra are obtained directly from the computer using the complex-valued entry $P_{K}^{+} / Q_{K}^{+}$. Alternatively, however, the analytical expression for e.g. ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$ can be derived first in the form of two partitions $A_{K}^{+}$ and $B_{K}^{+}$. Similarly, the analytical expression for ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$ also contains its own two partitions, $C_{K}^{+}$ and $D_{K}^{+}$. The sum of the partitioned spectra $A_{K}^{+}$ and $B_{K}^{+}$ is the complete absorption partitioned spectrum, $A_{K}^{+}+B_{K}^{+},$ so that, theoretically, we must have $A_{K}^{+}+B_{K}^{+}= {\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$. This is because $A_{K}^{+}$ and $B_{K}^{+}$ are uniquely extracted from ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+}).$ Likewise, when the partitioned spectra $C_{K}^{+}$ and $D_{K}^{+}$ are added together, the complete dispersion partitioned spectrum is $C_{K}^{+} + D_{K}^{+}.$ This latter sum, in theory, must satisfy the relation $C_{K}^{+} + D_{K}^{+}= {\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$, since the unique partition of ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+}),$ in fact, led to $C_{K}^{+}$ and $D_{K}^{+}$. Thus, our nomenclature is to call ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$ the complete absorptive partitioned envelope when generated by way of the sum $A_{K}^{+}+B_{K}^{+}.$ By the same token, ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$ is termed the complete dispersive partitioned envelope if it is obtained from the sum $C_{K}^{+} + D_{K}^{+}.$ On the other hand, when ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$ and ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$ are computed directly without partitioning, they will be referred to as the conventionally (or customarily) generated non-partitioned envelopes. In the case of envelope partitioning, it is the analytical expressions for $A_{K}^{+},\, {B}_{K}^{+},\, {C}_{K}^{+}$ and $D_{K}^{+}$ that we feed separately into the computer. The ensuing numerical results are graphed to visualize the partitioned absorption envelopes $A_{K}^{+}$ and $B_{K}^{+}$, as well as the partitioned dispersion envelopes $C_{K}^{+}$ and $D_{K}^{+}$. The partitions $A_{K}^{+}$ and $B_{K}^{+}$ in ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$ as well as $C_{K}^{+}$ and $D_{K}^{+}$ in ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$ redistribute the interference effect. This is the mechanism behind separating overlapping peaks in composite resonances. It is the interference of $A_{K}^{+}$ and $B_{K}^{+}$ in $A_{K}^{+}+B_{K}^{+}$ that, in fact, prevents splitting of adjacent overlapping resonances in a composite peak. In a rearranged interference, followed by plotting $A_{K}^{+}$ and $B_{K}^{+}$ separately, the individual resonances have a chance to “pop up” and, thus, split apart the compound peaks in $A_{K}^{+}+B_{K}^{+}$. We therefore computed the partial envelope spectra from $A_{K}^{+},\, {\, B}_{K}^{+},\, C_{K}^{+}$ and $D_{K}^{+}$. As a check, the results for the complete absorptive partition $A_{K}^{+}+B_{K}^{+}$ and the complete dispersive partition $C_{K}^{+} + D_{K}^{+}$ must be shown to coincide with the conventional non-partitioned absorption ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$ and conventional non-partitioned dispersion ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$ envelopes, respectively.

2.2.2 Partitioning for peering into the inner structure of total shape spectra

When $P_{K}^{+}$ and $Q_{K}^{+}$ become available and their quotient $P_{K}^{+} / Q_{K}^{+}$ directly programmed by this indicated division, the computer would give the entire contributions $\mathrm {Re}(P_{K}^{+} / {Q_{K}^{+})}$ and $\mathrm {Im}(P_{K}^{+} / {Q_{K}^{+})}$ as in (7), without the explicit information contained in (8)–(12). Customarily, computed in that way, all the intact interference effects in $\mathrm {Re}(P_{K}^{+} / {Q_{K}^{+})}$ and likewise in $\mathrm {Im}(P_{K}^{+} /Q_{K}^{+})$ are at play, but no insight can be gleaned into the contribution of the constituent parts of ${\mathrm {Re}}(P_{K}^{+} /Q_{K}^{+})$ and ${\mathrm {Im}}(P_{K}^{+} / Q_{K}^{+})$.

In contrast, insofar as $\mathrm {Re}(P_{K}^{+} /Q_{K}^{+})$ and $\mathrm {Im}(P_{K}^{+} / {Q_{K}^{+})}$ are identified by (8)–(12), before performing numerical computations, it becomes possible to look into the inner structure of these absorptive and dispersive envelopes, respectively. This inner structure is apparent in the partitioned envelopes $A_{K}^{+}$, $B_{K}^{+}$, $C_{K}^{+}$ and $D_{K}^{+}$ that are built from the products of the two spectra (at a time) in the MA and AR sequences. Consequently, the absorptive portion $\mathrm {Re}(P_{K}^{+} / {Q_{K}^{+})}\,$of the entire complex spectrum $P_{K}^{+} / Q_{K}^{+}$ is the sum $A_{K}^{+}+ B_{K}^{+}$ of the absorption-absorption ${(A}_{K}^{+})$ and dispersion-dispersion ${(B}_{K}^{+})$ products via $A_{K}^{+}={\mathrm {[Re}(P}_{K}^{+})]{\mathrm {Re}(1/Q}_{K}^{+})$ and ${B_{K}^{+}=-\mathrm {[Im}(P}_{K}^{+})]{\mathrm {Im}(1/Q}_{K}^{+}),$ as per (9) and (10), respectively. By the same token, the dispersion $\mathrm {Im}(P_{K}^{+} / {Q_{K}^{+})}$ of $P_{K}^{+} / Q_{K}^{+}$ is the sum $C_{K}^{+} + D_{K}^{+}$ of the cross-products (or mixed products) through absorptive-dispersive $C_{K}^{+}={\mathrm {[Re}(P}_{K}^{+})]{\mathrm {Im}(1/Q}_{K}^{+})$ and dispersive-absorptive $D_{K}^{+}={\mathrm {[Im}(P}_{K}^{+})]{\mathrm {Re}(1/Q}_{K}^{+})$ lineshapes, as per (11) and (12), respectively.

The full interference between the two partitioned envelopes is, in fact, redistributed via this compartmentalization of both $\mathrm {Re}\left( P_{K}^{+} / Q_{K}^{+}\right)$ and $\mathrm {Im}\left( P_{K}^{+}/Q_{K}^{+}\right)$. Consequently, the interference effect in $A_{K}^{+}$ and $B_{K}^{+}$ is diminished when each of these are examined separately. Hidden resonances can thereby be identified in compound peaks. Similarly, when $C_{K}^{+}$ and $D_{K}^{+}$ are taken as individual, partitioned envelopes, some concealed resonances can be seen that otherwise can be absent or obscured in $C_{K}^{+} + D_{K}^{+}$.

In order to check the correctness of the expressions (7)–(12), the numerical values of the spectrum $\mathrm {Re}\left( P_{K}^{+} / Q_{K}^{+}\right)$ computed conventionally, i.e. directly, must fully coincide with the sum $A_{K}^{+} + B_{K}^{+}$. Likewise, the direct computation of $\mathrm {Im}\left( P_{K}^{+} / Q_{K}^{+}\right)$ and that via $C_{K}^{+} + D_{K}^{+}$ must be entirely coincident.

Either of these two computations would merely result in superpositions of peaks with minimal interference effects for resonances that are single and isolated due to weak or virtually non-existent interactions. However, overlapping resonances can interact strongly, with the amount of their interference being related to the extent of overlap. With augmented interference, individual lineshapes of closely spaced resonances would be masked. Partitioned envelopes are then found to be helpful in disentangling the hidden spectral content of the component peaks. This expectation from the theory has been confirmed by the ${\mathrm {FPT}}^{(+)}$ in our initial study on noise-free synthesized time signals for a single value of model order K in $P_{K}^{+} / Q_{K}^{+}$ [1]. This finding is encouraging. Nevertheless, it needs to be further scrutinized in light of the known sensitivity of reconstructions to the presence of noise and to changes in model order K. Moreover, the reconstructed noisy envelopes for a sequence of values of model order K need to be averaged to improve SNR and damp the outliers (spikes in spectra).

It should be emphasized that, in principle, the concept of envelope partitioning is not limited to the FPT and, thus, can be tried with some other processors. However, what is particularly appealing in the FPT is that its envelope partitioning comes naturally as the unique decomposition into the two already ingrained constituents, the moving average and auto-regression. This great advantage of the Padé-based envelope partitioning permits realistic interpretations of absorptions and dispersions directly in terms of the familiar MA and AR processing, thus attaching physical meanings and rationale to the spectra partitioning concept.

We re-emphasize that both sets of Padé-based partitioned envelopes $\{A_{K}^{+},\, {B}_{K}^{+}\}$ and $\{C_{K}^{+},\, D_{K}^{+}\}$ represent the unique and exact decomposition of $\mathrm {Re}\left( P_{K}^{+} / Q_{K}^{+}\right)$ and $\mathrm {Im}\left( P_{K}^{+} / Q_{K}^{+}\right)$, respectively. In other words, there are no freely adjustable parameters in either of the two sets of partitioned envelopes. Of course, envelope partitioning is restricted to qualitative estimation alone. Such a restriction is expected from the onset since non-parametric processing can provide only shape estimation. This is a qualitative evaluation, which is descriptive regarding the components: i.e. “present” or “absent”. It is important to note this feature in order to avoid any attempt to view envelope partitioning as a substitute (or alternative) to quantification. It is not. Instead, it should be considered as a practical aid (or supplement) to quantification. Scanning, i.e. monitoring envelopes in a pre-quantitative stage could be used to spot the regions of high spectral density and determine whether compound peaks are decomposed. The partitioned envelopes would thus play the role of an adaptive pre-processor, laying the ground for a more focused local quantification. Taken in this way, both envelope partitioning and local quantification have clearly defined, complementary tasks.

We now proceed to briefly describe the spectra averaging procedure, an important application of the FPT, which can be used to further refine the information gleaned from the total shape spectra.

2.3 Spectra averaging through the FPT

A serious problem which has been noted to arise with spectra in MRS is the destabilizing effect of marked changes in the sought model order K [24,25,26,27,28]. Namely, for different values of model order K, many large noise-like spikes appear in the spectra. Earlier, we referred to these spikes as outliers, because they single themselves out from the more or less coherent pattern of most of the spectral peaks. We introduced an averaging procedure of spectra to handle this over-sensitivity to alterations in K [24]. This procedure regularizes spectra by taking the arithmetic average of a set of envelopes computed for selected values of K. Using a sequence of values of the model order K, an average envelope (arithmetic average) is produced.

This spectra averaging procedure is carried out within the FPT. Through the FFT, averaging of spectra is not possible since Fourier vectors in the frequency domain are not of the same length for different truncation of the total acquisition time T [24].

Spectra averaging can be performed iteratively. Thereby, the complex average envelope is inverted to generate the 1st reconstructed complex FID. The latter is then subjected to the FPT to generate the next set of envelopes for the same sequence of values of K as considered in the previous iteration. This new sequence of envelopes is averaged and the procedure can be repeated until the prescribed accuracy of the reconstructed spectral parameters is attained. With each iteration, there is progressively greater suppression of spurious spectral structures [24].

Whether performed once or iteratively, spectra averaging is an effective and unbiased strategy for alleviating the effect of redundancy and unphysical degrees of freedom from the reconstructions. This problem stems from the non-coherent part of the extracted information, present in the input encoded FIDs and also stemming from the reconstructed data (computational round-off errors, unstable recovered resonances, etc.). Through spectra averaging for a number of values of K, these random errors can be corrected in the stabilized quotient $P_{K}^{+}/Q_{K}^{+}$. This reflects error self-correction through the unique coupling of averaging of Padé spectra and the rational function response of the examined system to external perturbations. Heretofore, spectra averaging through the FPT has been performed for in vivo encoded MRS signals [24,25,26,27,28].

We now orient the above-outlined considerations concerning signal processing to a specific clinical problem in MRS, for which the envelope-partitioning capabilities of the FPT are particularly salient. Namely, through the ${\mathrm {FPT}}^{(+)}$, we address the question of how non-parametric analysis can visually identify the presence of a recognized breast cancer biomarker.

2.4 Clinical motivation: Towards an effective strategy for detection of the breast cancer biomarker phosphocholine through MRS

The clinical motivation for the present study is to help develop an effective strategy for identifying phosphocholine, PC, which is a molecular indicator of breast cancer. The importance of this problem is underscored by the fact that breast cancer is the most frequently occurring cancer and cause of cancer-related deaths among women worldwide [30, 31]. Early detection unequivocally contributes to better survival [31,32,33]. Further improvements in diagnostic accuracy over the anatomic imaging methods via magnetic resonance imaging (MRI), which is very sensitive but insufficiently specific, are anticipated by molecular imaging through MRS [34]. With the FFT processing for in vivo MRS and MRSI, somewhat higher specificity has been achieved. This has mainly been through assessment of a single peak, total choline (tCho), located at chemical shift $\sim$ 3.2 parts per million (ppm). Reports have been published on assessment of tCho via in vivo MRS for several hundred breast lesions, but pooled estimates of sensitivity and specificity have not ever surpassed 90% [35,36,37,38]. As summarized in more detail in Ref. [1], based upon tCho assessments through in vivo MRS, there are still no sufficiently trustworthy standards to diagnose a breast lesion as cancerous versus benign.

In vitro nuclear magnetic resonance (NMR) spectroscopic studies of extracted specimens show that within the tCho peak are not only phosphocholine, PC, at $\sim$ 3.22 ppm, but also glycerophosphocholine (GPC) at $\sim$ 3.23 ppm and free choline (Cho) resonating at $\sim$ 3.21 ppm [39]. A high PC to GPC concentration ratio and increased PC levels are associated with malignant transformation of the breast, as well as of other tissues [40,41,42], possibly related to loss of tumor suppressor p53 function [43]. The much more abundant phosphoethanolamine, PE, which is not a cancer biomarker, also resonates at $\sim$ 3.22 ppm [44]. For synthesized noiseless and noise-corrupted FIDs akin to in vitro MRS time signals encoded from normal breast, fibroadenoma and breast cancer, as per Ref. [44], the standard non-parametric formula for $\mathrm {Re}(P_{K}/Q_{K})$ was applied to visualize the total shape spectrum in the FPT [11, 16, 19,20,21]. Thereby, the large peak centered at $\sim$ 3.22 ppm appeared as a pure Lorentzian, without any suggestion whatsoever that underneath PE was a PC resonance. With the parametric FPT, applied to the same synthesized noise-free and noise-corrupted MRS time signals from Ref. [44], both PC and PE were identified and precisely quantified, as was also the case with the seven other resonances for all three types of breast tissue [11, 16, 19,20,21, 45]. In Ref. [11], we explored in detail how these results could contribute to a more individualized approach to breast cancer diagnosis and treatment, coherent with the aims of “personalized cancer medicine” (PCM).

2.4.1 PC and PE resolved in partitioned envelopes through the FPT for breast cancer: The noiseless case

Most recently in Ref. [1], we computed absorptive $\mathrm {Re}\left( P_{K}^{+} / Q_{K}^{+}\right)$ and dispersive $\mathrm {Im}(P_{K}^{+} / Q_{K}^{+})$ envelope spectra from their corresponding partitioned analytical expressions as described in subsection 2.2. The question was raised as to whether a clear indication of the presence of underlying resonances in some composite spectral structures could be ascertained from these alternative, non-parametric representations despite their default restriction to estimations of lineshapes alone. This was carried out using a single value of model order K in the spectral region between 3.2 and 3.3 ppm for synthesized, noise-free MRS time signals associated with breast cancer according to the in vitro data of Ref. [44]. In the absorptive partitioned envelopes, both PC and PE were clearly distinguished as two separate peaks. In fact, they were so well delineated that in the absorption mode, a dip between them descended all the way down to the baseline. In the dispersion mode via the partitioned envelopes, both PC and PE were also readily identifiable. In sharp contrast, without partitioning, the non-parametric FPT generated a single composite smooth Lorentzian peak (PC $+$ PE) in the absorptive envelopes, without any indication that two resonances PC and PE were present therein. This was the first study applying the non-parametric FPT in the partitioned manner, i.e. without resorting to quantification at all. The conclusion was that this line of investigation was fruitful. A qualitative comparison with the kth component shape spectrum $\left( 1\le k\le K \right) ,$ parametrically computed via ${\mathrm {Re}(P_{K}^{+} / Q_{K}^{+}\mathrm {)}}_{k}$ and ${\mathrm {Im}(P_{K}^{+} / Q_{K}^{+}\mathrm {)}}_{k},$ confirmed the reliability of this approach for the noiseless case at a single model order K. Here, as before, the adjective “qualitative” is used to refer to the situations “PC and PE visualized as two separate lineshapes” without regard to their correct relative abundance, which is the subject of quantification.

2.5 Aims of the present study

This is the second study on the sub-topic of partitioned envelopes, where we continue our investigation of the capabilities of the non-parametric ${\mathrm {FPT}}^{(+)}$ prior to quantification of noisy synthesized time signals. We seek to ascertain whether phosphocholine, PC, and phosphoethanolamine, PE, are still identifiable on the envelopes generated non-parametrically with partitioning via the ${\mathrm {FPT}}^{(+)}$, if progressively greater levels of noise were added to the noiseless MRS time signals associated with breast cancer. We also examine whether PC and PE are recognizable across a wide range of model orders K on the partitioned envelopes. Spectra averaging is implemented to improve SNR and to determine whether the findings remain robust. In particular, it is necessary to check if the clear splitting of PC and PE found in each of the partitioned envelopes would persist once spectra have been averaged over different model orders K. Such a check is needed because spectra averaging could mask or wash out the splitting of PC and PE if the lineshapes of the partitioned envelopes would exhibit huge changes for varying K. These aims constitute the most thorough validation of the concept of non-parametric envelope partitioning within the ${\mathrm {FPT}}^{(+)}$ in the controlled setting with simulated FIDs. The clinical incentive of this study is to justify implementing this strategy for in vivo MRS for a wide range of concerns within the realm of diagnosis and treatment of breast cancer and beyond.

3 Methods

3.1 The input data: Time signals based on in vitro encoding from breast cancer

Following the quantum-mechanical form of autocorrelation functions (arising directly from the general time evolution operator):

$$\begin{aligned} c_{n} =\sum \limits _{k=1}^K {d_{k} z_{k}^{n} } , z_{k} =\text{ e }^{i\tau \omega _{k} }, 0\le n\le N-1,\quad \hbox { Im}(\omega _{k} )>0,\quad \omega _{k} =2\pi \nu _{k} , \end{aligned}$$

(13)

we synthesized the MRS time signals as the sum of the indicated complex-valued attenuated exponentials. Here, the positive integer N denotes the total signal length. These data are based upon the encoded MRS time signals from cancerous breast, as per Ref. [44]. For concreteness, and following Ref. [44], each signal point $c_{n}$ from (13) is the sum of $K = 9$ damped complex exponentials $\exp \left( in\tau \omega _{k} \right) \, (1\le k\le 9)$. The angular frequencies $\omega _{k}$ are also complex, with an exponential decrease in $c_{n}$ over time $n\tau \, ( n=0,1,2,\ldots , N-1).$ The amplitudes $d_{k}$ are generally complex, as well. The corresponding FIDs of total length $N = 65536$ in Ref. [44] were encoded at a Larmor frequency $(\nu _{\mathrm {L}}$) of 600 MHz, which corresponds to static magnetic field of strength $B_{0}\approx 14.1\mathrm {T}$. The bandwidth (BW) was 6 MHz, where the inverse of this bandwidth is the sampling time $\tau$.

We use only one quarter of the above full signal length, namely $N/4=16384$. The median metabolite concentrations were based upon 14 samples (with two exceptions) taken from twelve patients (2 samples were taken from two of the patients). However, for $\beta$-glucose ($\beta$-Glc) and myoinositol (m-Ins) metabolite concentrations were available for only 6 and 9 samples, respectively [44]. The values of $\left| d_{k} \right|$ were derived from the median concentrations for the breast cancer input data, using the relation $\left| d_{k} \right| =2\mathrm {C}_{\mathrm {met}}/\mathrm {C}_{\mathrm {ref}}$, where $\mathrm {C}_{\mathrm {ref}}$ is the reference concentration ($\mathrm {C}_{\mathrm {ref}} = 0.05\,\hbox {mM}/\hbox {g}$). We denoted the mean metabolite concentration by $\mathrm {C}_{\mathrm {met}}$. The $T_{2}^{*}$ relaxation times were not specified in Ref. [44], so we set the imaginary parts of the complex frequencies $\nu _{k}$ to be 0.0008 ppm. The peaks are all Lorentzian, in accordance with the time signal from (13). The input amplitudes are defined as real, $d_{k}=\left| d_{k} \right|$, since the phases $\varphi _{k}\, (1\le k\le 9)$ from generally complex-valued $d_{k}$, were all set to zero. The string of the input metabolites and their fundamental parameters are:

$$\begin{aligned}&\left. \begin{array}{l} {\text{ Met }_{k} =\{\hbox {Lac}, \text { Ala,} \text { Cho,} \text { PC,} \text { PE,} \text { GPC,}\,\beta \text{-Glc, } \text { Tau,} \text { m-Ins}\}} \\ {\text{ Re( }\nu _{k} )=\, \{1.332,1.471,3.212,3.220,3.221,3.232,3.251,3.273,3.281\}\, \text{ ppm }}\\ {\text{ Im( }\nu _{k} )=0.0008\, \text{ ppm }\, (1\le k\le \text{9) }} \\ {d_{k} =\left| {d_{k} } \right| =\{0.325,0.032,0.004,0.012,0.090,0.009,0.029,0.112,0.036\}\, \text{ au }} \\ \end{array}\right\} \quad \nonumber \\&\quad \left( 1\le k\le 9 \right) , \end{aligned}$$

(14)

where the acronym “au” denotes arbitrary units. The internal reference (ref) was TSP (3-(trimethylsilyl-) 3,3,2,2-tetradeutero-propionic acid), a molecule which is, in fact, not present in the tissue. Thus, $\left| d_{k} \right| =\mathrm {C}_{\mathrm {met}}/(\mathrm {25}\mu \mathrm {M/g})$ of wet weight (ww).

We focus on the frequency band between 3.2 and 3.3 ppm, wherein PC and PE resonate at 3.220 and 3.221 ppm, respectively, as per (14) i.e. they are separated by only 0.001 ppm. Seven of the mentioned resonances lie within this spectral region of interest (SRI) for sweep frequencies $\nu \in [ 3.2,\, 3.3 ]$ ppm. The remaining two of the nine resonances, lactate (Lac) at 1.332 ppm and alanine (Ala) at 1.471 ppm, reported in Ref. [44] are outside our current SRI. Within this SRI, the largest concentration is that of taurine (Tau).

3.2 Addition of noise

Prior to reconstructions, the noiseless input time signals were mixed with added complex valued random zero-mean Gauss-distributed white noise of selected levels or standard deviations ($\sigma =$ 0.0289, $\sigma =$ 0.289 and $\sigma =$ 2.89 RMS). Here, the acronym RMS refers to the ordinary root-mean-square, RMS $=(1/N)\, \sum \nolimits _{n=0}^{N-1} \left| c_{n} \right|$, where $\{c_{n}\}$ is the noiseless FID from (13). We express the noise level in a noisy FID by percentage of the RMS of the noiseless time signal $\{c_{n}\}$. For example, stating that a noisy FID contains 3% of noise means that “the noisy RMS” is 3% of “the noiseless RMS”, or equivalently, RMS (FID with 3% noise level) $=$ 0.3 RMS (noiseless FID). We measure the noise level in terms of RMS because this noise quantifier minimizes the bias relative to the actual (sought) value and the variance of noise. Furthermore, RMS converts signal oscillations to variations of the power of the signal across a given bandwidth. Thus, RMS is a measure of the dynamics of the signal [21, 22]. The ${\mathrm {FPT}}^{(+)}$ is herein applied to the noise-free and noise-corrupted MRS time signals.

3.3 Non-parametric computation of the envelopes with and without spectra partitioning

The envelopes will be computed as in the standard way, by directly feeding the set of extracted complex numbers $P_{K}^{+} / Q_{K}^{+}$ into the computer to generate the pure absorptive, ${\mathrm {Re}}\left( P_{K}^{+} / Q_{K}^{+}\right)$, and dispersive, ${\mathrm {Im}}\left( P_{K}^{+} / Q_{K}^{+}\right) ,$ spectral lineshapes. In this way, the computer generates both spectra ${\mathrm {Re}}(P_{K}^{+} / Q_{K}^{+})$ and ${\mathrm {Im}}\left( P_{K}^{+} / Q_{K}^{+}\right)$ in their entirety as non-partitioned, i.e. composite envelopes. These conventionally computed spectra are the reference envelopes for establishing the validity of their partitioned counterparts.

As described in Sect. 2.2.1, we will obtain the separate partitioned products in the form of the quantities $A_{K}^{+}$, $B_{K}^{+}$, $C_{K}^{+}$ and $D_{K}^{+}$, corresponding, respectively to (9), (10), (11) and (12). The spectra $\mathrm {Re}\left( P_{K}^{+} / Q_{K}^{+}\right)$ computed directly and via $A_{K}^{+}\, + B_{K}^{+}$ will be checked for their expected coincidence. Similarly, the direct computation of $\mathrm {Im}\left( P_{K}^{+} / Q_{K}^{+}\right)$ and via the sum $C_{K}^{+} + D_{K}^{+}$ will also be compared to verify agreement of these two quantities.

3.4 Examination of several model orders and spectra averaging

Herein, we examine 6 model orders K with an increment of 500, ranging from 2500 to 5000. This is written from here on according to the mathematical convention $K=2500\left( 500 \right) 5000$ with the increment $\Delta K=500$ written in the parentheses between $K_{\mathrm {min}}=2500$ and $K_{\mathrm {max}}=5000$ [28]. The very large increment was chosen to scrutinize sensitivity of reconstructions to vastly different values of the model orders because their huge changes are likely to produce markedly different lineshapes of partitioned envelopes [26]. The arithmetic average of all the envelopes (partitioned and non-partitioned) for these 6 model orders will be computed, and denoted with the subscript Av, e.g. ${{\{X}_{K}^{+}\}}_{\mathrm {Av}}\, \left( X=A,B,C,D \right) ,\, {\{\mathrm {Re}(P_{K}^{+}/Q_{K}^{+}\}}_{\mathrm {Av}}$ and ${\{\mathrm {Im}(P_{K}^{+}/Q_{K}^{+}\}}_{\mathrm {Av}}$. It should be noted that the reconstructions are performed using partial signal lengths $N_{\mathrm {P}} =$ 5000 ($K = 2500$) to $N_{\mathrm {P}} =$ 10 000 ($K\, =\, 5000)$ of full signal length $N=16384$, and these are approximately thirty to sixty percent of its encoded counterpart of 65536 FID data points from Ref. [44]. In order to check the dependence of the results upon the choice of the interval for K and $\Delta K$, we have also carried out computations of partitioned envelopes for $\Delta K=250$ in the following three sets: $K=250\left( 250 \right) 5000,\,\, K=1000\left( 250 \right) 5000\, \mathrm {and} \,\, K=2500\left( 250 \right) 5000$ with some 20, 17 and 11 values of K,  respectively. The findings from each of these sets were fully consistent with those from the set $K=2500\left( 500 \right) 5000$ presented in the Results section.

4 Results

The level of added noise is $\sigma =0.0289\,\hbox {RMS}$ for Figs. 1–8, whereas for Fig. 9, it is $\sigma =0.289\,\hbox {RMS}$ and $\sigma =2.89\,\hbox {RMS}$. Figure 1 presents the partitioned (noisy) and non-partitioned (noiseless) absorption and dispersion envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$, within the critical frequency window [3.215, 3.225] ppm containing the overlapping PC and PE resonances. All six model orders $K=2500\left( 500 \right) 5000$ are displayed with color coding as black $(K=2500)$, green $(K=3000)$, cyan $(K=3500)$, red $(K=4000)$, magenta $(K=4500)$ and blue $(K=5000)$. The absorption spectra are presented in the left column, whereas the right column shows the dispersion spectra. Along the abscissae of each panel are the input chemical shifts, that are symbolized by filled circles only when the PC and/or PE lineshape peaks are located practically at their correct chemical shifts 3.220 and 3.221 ppm, respectively.

Fig. 1

Partitioned and non-partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.215, 3.225] ppm containing PC and PE. Use of the FID sampled at $N=16384$ with noise level $\sigma =0.0289\,\hbox {RMS}$ and model orders $K=2500\left( 500 \right) 5000,\,$color coded as black $(K=2500)$, green $(K=3000)$, cyan $(K=3500)$, red $(K=4000)$, magenta $(K=4500)$ and blue $(K=5000)$. Panel (a): the noisy partitioned absorptions for $A_{K}^{+}$ (PC correct, filled circle). Panel (b): the noisy partitioned absorptions for $B_{K}^{+}$ (PE correct, filled circle). Panel (c) for the 6 values of K: the noisy complete partitioned absorptions $A_{K}^{+}+B_{K}^{+}$ and noiseless non-partitioned (conventional) absorptions $\mathrm {Re}( P_{K}^{+}/Q_{K}^{+} )$. These 12 curves are indistinguishable. Panel (d): the noisy partitioned dispersions for $C_{K}^{+}$ (PC correct, filled circle). Panel (e): the noisy partitioned dispersions for $D_{K}^{+}$ (PE correct, filled circle). Panel (f) for the 6 values of K: the noisy complete partitioned dispersions $C_{K}^{+}+ D_{K}^{+}$ and noiseless conventional dispersions $\mathrm {Im}( P_{K}^{+}/Q_{K}^{+} )$. These 12 curves are indistinguishable. (Color online)

In panel (a) of Fig. 1 are the noisy partitioned absorption envelope spectra for $A_{K}^{+}.$ Resonance PC at 3.220 ppm is indicated with a filled circle, and clearly seen at nearly the correct chemical shift position for all the model orders, albeit with slight variations among these six different values of K. Greater variation is observed at 3.221 ppm, with the PE peak marked by an open circle signifying a slight displacement. There is also more variation in the PE lineshapes among the 6 model orders. Still, the PE peak is clearly present for each of the model orders. In panel (b) are the noisy partitioned absorption envelope spectra for $B_{K}^{+}$. Resonance PE at 3.221 ppm is a filled circle at the practically correct position and the variation among the 6 model orders is somewhat less pronounced than for the corresponding peaks in panel (a). This time, the circle at 3.220 associated with PC is open since that resonance is slightly shifted to the right in panel (b). Nevertheless, also in this panel, PC is well visualized for all the 6 model orders. In panel (c), the complete noisy absorption envelopes $A_{K}^{+}+ B_{K}^{+}$ are displayed for the 6 model orders. In addition, the noiseless non-partitioned (conventional) absorption envelopes $\mathrm {Re}(P_{K}^{+}/Q_{K}^{+})$ are presented therein for the 6 model orders. Thus, altogether there are 12 curves that fully coincide, appearing as an entirely blue single curve, which is drawn last for $K=5000$. However, in the noisy complete partitioned absorption envelopes $A_{K}^{+}+ B_{K}^{+}$, just like in the noiseless conventional absorption envelopes $\mathrm {Re}(P_{K}^{+}/Q_{K}^{+})$, there is no suggestion whatsoever of the two resonances PC and PE. Rather, a single smooth Lorentzian appears. Thus, in panel (c), both PC and PE are symbolized by their open circles.

Panel (d) of Fig. 1 shows the noisy partitioned dispersion envelope spectra for $C_{K}^{+}$ with the 6 model orders. Therein, the PC dip at 3.220 ppm is a filled circle, with slight discordance among the 6 model orders. Also, the PE is slightly pushed relative to 3.221 ppm and, thus, the associated circle is open. Although there is more noticeable variance associated with the displaced PE compared to PC, the former resonance can still be seen for each of the 6 model orders. In panel (e) are the noisy partitioned dispersion envelopes for ${D}_{K}^{+}$, where PE at 3.221 ppm is a filled circle and the PC circle at 3.220 is open, due to a displacement. Nevertheless, the PC resonance is still identifiable for all the model orders. In panel (f), the noisy complete dispersion envelopes $C_{K}^{+}+ D_{K}^{+}$ are shown with the 6 model orders. The noiseless conventional, non-partitioned dispersion envelopes $\mathrm {Im}(P_{K}^{+}/Q_{K}^{+})$ are also depicted therein for the 6 model orders. In panel (f), there are actually 12 curves that appear as one single blue line (the last drawn for $K=5000)$, without indication of the two component resonances PC and PE. Consequently, in panel (f), both PC and PE are symbolized by open circles. Importantly, full agreement between the noisy (partitioned) and noiseless (non-partitioned) spectra, as seen in panels (c) and (f) of Fig. 1, demonstrates the noise suppression ability of the ${\mathrm {FPT}}^{(+)}.$

Fig. 2

Partitioned and non-partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.205, 3.290] ppm. Use of the FID sampled at $N=16384$ with noise level $\sigma =0.0289\,\hbox {RMS}$ and model orders $K=2500\left( 500 \right) 5000$, color coded as in Fig. 1. Panel (a): the noisy partitioned absorptions for $A_{K}^{+}$ (PC correct, filled circle); the other input chemical shifts are represented by open circles. Panel (b): the noisy partitioned absorptions for $B_{K}^{+}$ (PE correct, filled circle). Panel (c) for the 6 values of K: the noisy complete partitioned envelope $A_{K}^{+}+ B_{K}^{+}$ and noiseless non-partitioned (conventional) absorption envelopes $\mathrm {Re}( P_{K}^{+}/Q_{K}^{+} )$. These 12 curves are indistinguishable. Panel (d): the noisy partitioned dispersions for $C_{K}^{+}$ (PC correct, filled circle). Panel (e): the noisy partitioned dispersions for $D_{K}^{+}$ (PE correct, filled circle). Panel (f) for the 6 values of K: the noisy complete partitioned dispersions $C_{K}^{+}+ D_{K}^{+}$ and noiseless conventional dispersions $\mathrm {Im}( P_{K}^{+}/Q_{K}^{+} )$. These 12 curves are indistinguishable. (Color online)

Figure 2 displays the partitioned (noisy) and non-partitioned (noiseless) absorption and dispersion envelopes computed non-parametrically by the ${\mathrm {FPT}}^{(+)}$ for a wider SRI which is enlarged to include chemical shifts between 3.205 and 3.290 ppm. Again, the six model orders $K=2500\left( 500 \right) 5000$ are used, with color coding as in Fig. 1. Along the abscissae of each panel are the input chemical shifts in this extended SRI. These are symbolized by circles that are mainly open. The two exceptions with the filled circles relate to PC and PE. Circles for 5 isolated non-overlapping resonances (Cho, GPC, $\beta$-Glc, Tau, m-Ins) are systematically kept as open, irrespective of whether or not their positions of the reconstructed metabolites are close to the input chemical shifts. This practice is introduced because the primary focus here is on the overlapping PC and PE resonances. Thus, in panels (a), (b), (d) and (e), resonances PC and PE are marked with filled circles. Therein, two filled circles are not simultaneously present on any panel. For example, on panel (a), the location of the center of the reconstructed lineshape for PC almost coincides with the input chemical shift (3.220 ppm) and, thus, appears as a full circle. However, on the same panel (a), the PE resonance is seen as being slightly displaced from its exact position (3.221 ppm) and, therefore, the location of this lineshape is marked with an open circle. This pattern is reversed on panel (b), where the PE resonance is retrieved at its correct location (3.221 ppm) and, hence, marked by a filled circle. On the other hand, panel (b) shows an open circle for the PC resonance due to its slight displacement from the correct chemical shift (3.220 ppm). Panels (a) and (b) are for the noisy absorption spectra $A_{K}^{+}$ and $B_{K}^{+}$, respectively. A similar situation with the appearance of filled and open circles for the PC and PE resonances is also encountered for the noisy dispersion spectra $C_{K}^{+}$ and $D_{K}^{+}$ in panels (d) and (e), respectively. Overall, in panels (a), (b), (d) and (e), there is quite close concordance among the partitioned envelopes for the 6 model orders, with noticeable, but fairly minimal discrepancies seen at PE, PC, and elsewhere along the SRI. Nevertheless, all the resonances, including PC and PE are clearly resolved in each of the panels displaying noisy partitioned envelopes for every model order K. Taurine and $\beta -\mathrm {Glc}$ show much smaller peak heights in the partitioned envelope spectra for $B_{K}^{+}$ compared to that for $A_{K}^{+}$ and in $D_{K}^{+}$ compared to that for $C_{K}^{+}$. For the 6 model orders $K=2500\left( 500 \right) 5000$ in panel (c), the noisy complete partitioned absorption envelopes $A_{K}^{+}+ B_{K}^{+}$ and the noiseless non-partitioned (conventional) absorption envelopes $\mathrm {Re}(P_{K}^{+}/Q_{K}^{+})$ are displayed. The resultant 12 curves are fully concordant, and appear as a single blue curve (the last plotted for $K=5000)$. The peak at $\sim$3.22 ppm is a smooth, regular Lorentzian. Analogously, in panel (f) for the 6 model orders $K=2500\left( 500 \right) 5000$, the noisy complete dispersion envelopes $C_{K}^{+}+ D_{K}^{+}$, and the noiseless conventional dispersion envelopes $\mathrm {Im}(P_{K}^{+}/Q_{K}^{+})$ are shown, appearing as a single blue curve, which is also the last plotted for $K=5000$. In panels (c) and (f), there is no suggestion at all of the two resonances PC and PE. Thus, both PC and PE are open circles in both these panels. For this larger SRI, [3.205, 3.290] ppm, the ${\mathrm {FPT}}^{(+)}$ efficiently suppresses noise as evidenced by full agreement between the noisy (partitioned) and noiseless (non-partitioned) envelopes in panels (c) and (f) of Fig. 2.

Figure 3 shows only noisy partitioned envelopes. Therein, displayed are the absorption spectra, sequential (left column) and averaged (right column), reconstructed by the ${\mathrm {FPT}}^{(+)}$ on the critical frequency window [3.215, 3.225] ppm. On the left side, for the six model orders $K=2500\left( 500 \right) 5000,$ color coded as previously, are the partitioned absorption envelope spectra for $A_{K}^{+}$ in panel (a), for $B_{K}^{+}$ in panel (b) and for the complete envelope $A_{K}^{+}+ B_{K}^{+}$ in panel (c). In other words, we are basically repeating the left column of Fig. 1, except that in panel (c) only the complete envelopes $A_{K}^{+}+ B_{K}^{+}$ for the 6 model orders are actually present. This repetition is done to remind us about the source used in spectra averaging. On the right side of Fig. 3 are the averaged partitioned spectra for ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f). On panel (d) for ${{\{A}_{K}^{+}\}}_{\mathrm {Av}},$ the PC lineshape is correctly placed at 3.220 ppm. Therein, the PE peak, while slightly displaced to the left, is a well-defined resonance separated from PC almost all the way down to the baseline. Analogously, on panel (e) for ${{\{B}_{K}^{+}\}}_{\mathrm {Av}},\,$both PC and PE are seen as entirely separate peaks, with PE correctly placed at 3.221 ppm and larger than PC, which is pushed slightly rightward. Panel (f) with the complete envelope ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ coincides with the complete envelopes $A_{K}^{+}+ B_{K}^{+}$ for the six model orders $K=2500\left( 500 \right) 5000$ in panel (c).

Fig. 3

Partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.215, 3.225] ppm containing PC and PE. Use of the FID sampled at $N=16384$ all with noise level $\sigma =0.0289\,\hbox {RMS}$ and model orders $K=2500\left( 500 \right) 5000$, color coded as in Fig. 1. Panel (a): the partitioned absorptions for $A_{K}^{+}$ (PC correct, filled circle). Panel (b): the partitioned absorptions for $B_{K}^{+}$ (PE correct, filled circle). Panel (c) for the 6 values of K: the complete partitioned envelope $A_{K}^{+}+ B_{K}^{+}$. On the right side are the average envelopes ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for the average complete envelope ${{\{A}_{K}^{+}+\mathrm {\, }B_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f). These latter averages are from the envelopes on the left column. (Color online)

Fig. 4

Partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.205, 3.290] ppm. Use of the FID sampled at $N=16384$ with noise level $\sigma =0.0289\,\hbox {RMS}$ and model orders $K=2500\left( 500 \right) 5000$, color coded as in Fig. 1. Panel (a): the partitioned absorptions for $A_{K}^{+}$ (PC correct, filled circle). Panel (b): the partitioned absorptions for $B_{K}^{+}$ (PE correct, filled circle). Panel (c): the complete partitioned envelope $A_{K}^{+}+ B_{K}^{+}$ for the 6 model orders of K. On the right side are the average envelopes ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for the average complete absorptions ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f). These latter averages are from the envelopes on the left column. (Color online)

Figure 4 also depicts only noisy partitioned envelopes. Herein, we continue with spectra averaging of the absorption envelopes reconstructed by the ${\mathrm {FPT}}^{(+)}$ in the extended SRI between 3.205 and 3.290 ppm. On the left side for the 6 model orders $K=2500\left( 500 \right) 5000,$ color coded as before, are the partitioned absorption envelopes for $A_{K}^{+}$ in panel (a), for $B_{K}^{+}$ in panel (b) and for the complete envelopes $A_{K}^{+}+ B_{K}^{+}$ in panel (c). This is a repeat of the left column of Fig. 2, except that here in panel (c), the noiseless case for the conventional absorption envelopes is not included. The average partitioned spectra for ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f) are displayed on the right side of Fig. 4. In panel (d), both PC (correctly placed at 3.220 ppm) and PE (slightly displaced to the left) are well-defined, sharply separated resonances. All the other resonances are also clearly delineated, with Tau by far the most prominent. For ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e), both PC and PE are completely distinct peaks, with PE correctly placed at 3.221 ppm and more prominent than PC, which is displaced slightly rightward. Taurine, m-Ins and $\beta$-Glc are attenuated in panel (e), compared to panel (d). The average complete envelope ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f) is indistinguishable from the sequential complete envelopes $A_{K}^{+}+ B_{K}^{+}$ for the 6 model orders $K=2500\left( 500 \right) 5000$ from panel (c).

Fig. 5

Partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.215, 3.225] ppm containing PC and PE. Use of the FID sampled at $N=16384$ with noise level $\sigma =0.0289\,\hbox {RMS}$ for model orders $K=2500\left( 500 \right) 5000$, color coded as in Fig. 1. Panel (a): the partitioned dispersions for $C_{K}^{+}$ (PC correct, filled circle). Panel (b): partitioned dispersions for $D_{K}^{+}$ (PE correct, filled circle). Panel (c) for the 6 values of K: the complete dispersions $C_{K}^{+}+ D_{K}^{+}$. On the right column are the average dispersions for ${{\{C}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{D}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for the average complete dispersions ${{\{C}_{K}^{+}+D_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f). These latter averages are from the envelopes on the left column. (Color online)

The noisy partitioned dispersion envelopes reconstructed by the ${\mathrm {FPT}}^{(+)}$ on the frequency window [3.215, 3.225] ppm are presented in Fig. 5. On the left side, for 6 model orders $K=2500\left( 500 \right) 5000,$ color coded as previously, are the noisy partitioned dispersion envelopes for $C_{K}^{+}$ in panel (a), for $D_{K}^{+}$ in panel (b) and for the complete envelopes $C_{K}^{+}+ D_{K}^{+}$ in panel (c). This repeats the right column of Fig. 1, with the exception that here in panel (c), the noiseless case for the conventional (non-partitioned) dispersion envelopes is not included. On the right side of Fig. 5 are the average partitioned spectra for ${{\{C}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{D}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for ${{\{C}_{K}^{+}+ D_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f). Panel (d) shows PC as being correctly located at 3.220 ppm. Therein, the PE lineshape displaced a bit to the left, is also a well-defined dip which is separated from PC. Similarly, in panel (e) for ${{\{D}_{K}^{+}\}}_{\mathrm {Av}}$ , both PC and PE are entirely distinct lineshapes, now with PE correctly placed at 3.221 ppm and larger than PC, which is pushed slightly rightward. Panel (f) with the complete average envelope ${{\{C}_{K}^{+}+ D_{K}^{+}\}}_{\mathrm {Av}}$ is precisely the same as the complete sequential envelopes $C_{K}^{+}+ D_{K}^{+}$ for the 6 model orders $K=2500\left( 500 \right) 5000$ in panel (c).

Fig. 6

Partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.205, 3.290] ppm. Use of the FID sampled at $N=16384$ with noise level $\sigma =0.0289\,\hbox {RMS}$ for model orders $K=2500\left( 500 \right) 5000$, color coded as in Fig. 1. Panel (a): the partitioned dispersions for $C_{K}^{+}$ (PC correct, filled circle). Panel (b): partitioned dispersions for $D_{K}^{+}$ (PE correct, filled circle). Panel (c) for the 6 values of K: the complete dispersions $C_{K}^{+}+ D_{K}^{+}$. On the right column are the average dispersions for ${{\{C}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{D}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for the average complete dispersions ${{\{C}_{K}^{+}+D_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f). These latter averages are from the envelopes on the left column. (Color online)

Figure 6 displays the noisy partitioned dispersion spectra reconstructed by the ${\mathrm {FPT}}^{(+)}$ without and with averaging, for the wider SRI between 3.205 and 3.290 ppm. On the left side are the 6 model orders $K=2500\left( 500 \right) 5000,$ color coded as earlier, with the partitioned dispersion envelopes for $C_{K}^{+}$ in panel (a), for $D_{K}^{+}$ in panel (b) and for the complete envelopes $C_{K}^{+}+ D_{K}^{+}$ in panel (c). This is a repeat of the right column of Fig. 2, except that here on panel (c), the noiseless case for the conventional (non-partitioned) dispersion envelopes is not given. On the right side of Fig. 6 are the averaged partitioned spectra for ${{\{C}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (d), for ${{\{D}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e) and for ${{\{C}_{K}^{+}+ D_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f). Panel (d) shows the PC resonance as being accurately placed at 3.220 ppm, with the PE resonance pushed a bit to the left. The two resonances, PC and PE, are entirely distinguishable. The other resonances in this fuller SRI are also distinct, with Tau being the dominant lineshape. For ${{\{D}_{K}^{+}\}}_{\mathrm {Av}}$ in panel (e), the PC and PE lineshapes are separate resonances. Therein, PE is correctly placed at 3.221 ppm and larger than PC, which is pushed rightward somewhat. Compared to panel (d), resonances Tau, m-Ins and $\beta$-Glc are attenuated in panel (e). The complete envelopes are identical with and without averaging, namely the average spectrum ${{\{C}_{K}^{+}+ D_{K}^{+}\}}_{\mathrm {Av}}$ in panel (f) is indistinguishable from the complete sequential envelopes $C_{K}^{+}+ D_{K}^{+}$ for the 6 model orders $K=2500\left( 500 \right) 5000$ in panel (c).

Fig. 7

Partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.205, 3.290] ppm. Use of the FID sampled at $N=16384$ for model orders $K=2500\left( 500 \right) 5000$. Panel (a): the partitioned absorptions for $A_{K}^{+}\, \hbox { with}\hbox { noise}\hbox { level }\sigma =0.0289\,\mathrm { RMS}$. Panel (b): the noiseless partitioned absorptions for $A_{K}^{+}$. Color coding for panels (a) and (b) as in Fig. 1. Panel (c): the average absorptions for ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ with $\hbox { noise}\hbox { level }\sigma =0.0289\,\mathrm { RMS}$ (blue) and noiseless (red). Averages in panel (c) computed from absorptions from panels (a) and (b). Panels (a−c): PC correctly located (filled circles). (Color online)

Figure 7 compares noisy and noiseless partitioned absorption envelopes. Therein, in panels (a), (b) and (c), we present the partitioned absorption envelopes $A_{K}^{+}$ and their averages ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ reconstructed by the ${\mathrm {FPT}}^{(+)}$ for the SRI of chemical shifts [3.205, 3.290] ppm with six model orders $K=2500\left( 500 \right) 5000$. In panel (a) are the partitioned absorption envelopes for $A_{K}^{+}$ with noise of level $\sigma =0.0289\,\hbox {RMS}$. The PC resonance is correctly placed and its spread is quite minimal, whereas the PE peak is pushed a bit leftward with more noticeable spread. Panel (b) shows the noiseless partitioned absorption envelopes for $A_{K}^{+}$ for the 6 model orders. Again, PC is correctly placed and PE is pushed leftward. Good agreement between the two sets of envelopes from panels (a) and (b) illustrates the powerful noise suppression capability of the ${\mathrm {FPT}}^{(+)}$. Spectra averaging for the 6 model orders is displayed in panel (c) with the results for the partitioned average absorption envelopes ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ indicated in red (noiseless) and blue (noisy, $\sigma =0.0289\,\hbox {RMS}$). These averages for the noise-free and noise-contaminated data, especially at 3.220 ppm corresponding to the correctly placed PC, are almost perfectly concordant, with also very good agreement seen for the other resonances in the SRI.

Fig. 8

Partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$ in the SRI with [3.205, 3.290] ppm. Use of the FID sampled at $N=16384$ for model orders $K=2500\left( 500 \right) 5000$. Panel (a): the partitioned absorptions for $B_{K}^{+} \hbox { with }\hbox { noise} \hbox { level }\sigma =0.0289\,\mathrm { RMS}$. Panel (b): the noiseless partitioned absorptions for $B_{K}^{+}$. Color coding for panels (a) and (b) as in Fig. 1. Panel (c): the average absorption for ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ with $\hbox { noise } \hbox { level }\sigma =0.0289\,\mathrm { RMS}$ (blue) and noiseless (red). Panels (a−c): PE correctly located (filled circles). Panel (d): the average complete absorption ${{\{A}_{K}^{+}+B_{K}^{+}\}}_{\mathrm {Av}}$ for $\hbox { noise }\hbox { level}\ \sigma =0.0289\,\mathrm { RMS}$ (blue), noiseless (red) and the noiseless non-partitioned (conventional) absorptions ${\mathrm {Re}( P_{K}^{+}/Q_{K}^{+} )}_{\mathrm {Av}}$ (green). All three curves are indistinguishable. Averages computed from the corresponding envelopes for the 6 values of model order K. (Color online)

Figure 8 compares the noisy and noiseless partitioned absorption envelopes $B_{K}^{+}$ as well as their averages ${{\{B}_{K}^{+}\}}_{\mathrm {Av}},$ as reconstructed by the ${\mathrm {FPT}}^{(+)}$ for the SRI of [3.205, 3.290] ppm with 6 model orders $K=2500\left( 500 \right) 5000$. For a cross–validation, also compared are the averaged partitioned ${{\{A}_{K}^{+}+B_{K}^{+}\}}_{\mathrm {Av}}$ and non-partitioned $\{{\mathrm {Re}(P_{K}^{+}/Q_{K}^{+})}\}_\mathrm{Av}$ absorption envelopes, where ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ is taken from Fig. 7c. In panel (a) are the partitioned absorption envelopes for $B_{K}^{+}$ with added noise of level $\sigma =0.0289\,\hbox {RMS}$. The PE resonance is correctly placed, whereas PC is slightly displaced rightward. Throughout the SRI, the spread among the lineshapes for different values of K is quite minimal. In the noiseless case in panel (b), the PE peak is also correctly placed, whereas PC is shifted rightward. Therein, the spread for PE is slightly more noticeable than that for the noisy case in panel (a). Spectra averaging for the 6 model orders is shown in panel (c), with the results for the partitioned average absorption envelopes ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ indicated in red (noiseless) and blue (noisy, $\sigma =0.0289\,\hbox {RMS}$). These averages for the noise-free and noise-contaminated data, especially at 3.221 ppm corresponding to the correctly placed PE are almost entirely aligned. For the other resonances in the SRI, concordance is also quite close in panel (c). Panel (d) of Fig. 8 represents the complete envelopes ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ for the noiseless (red) and noisy (blue, $\sigma =0.0289\,\hbox {RMS}$) data. These two curves are indistinguishable testifying to the noise-suppression feature of the ${\mathrm {FPT}}^{(+)}$. As such, whatever discrepancies in ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ and ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ might exist between the noiseless and noisy partitions, these disappear when the sum of ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ and ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ is taken via ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}+\, {{\{B}_{K}^{+}\}}_{\mathrm {Av}}$, for which the shortened notation ${{\{A}_{K}^{+}+B_{K}^{+}\}}_{\mathrm {Av}}$ is used. There is yet a third curve (green) plotted in panel (d) for the non-partitioned conventionally computed noiseless envelope ${\mathrm {Re}(P_{K}^{+}/Q_{K}^{+}\mathrm {\, )}}_{\mathrm {Av}}$. This latter result too blends together with its two partitioned counterparts. All the 3 curves entirely coincide and, thus, appear as a blue single curve, which is the last drawn. In the noisy absorption partitioned average envelope, ${{\{A}_{K}^{+}+B_{K}^{+}\}}_{\mathrm {Av}},$ and the noiseless conventional non-partitioned envelope, ${\mathrm {Re}(P_{K}^{+}/Q_{K}^{+})}_{\mathrm {Av}}$, there is no suggestion whatsoever of the two resonances PC and PE. Rather, a single smooth Lorentzian appears. Thus, both PC and PE are symbolized by open circles in panel (d). Spectra averaging was carried out to ascertain whether the findings about the PC splitting from PE remain robust for noisy input time signals. The outcome is that the partitioning strategy ensures unequivocal identification of PC and PE as two separate resonances for each of the 6 model orders. Spectra averaging for the 6 model orders provides clear delineation of the two resonances, PC and PE, in the average partitioned envelopes for both noiseless and noisy input time signals.

Figure 9 displays the absorption envelopes reconstructed by the ${\mathrm {FPT}}^{(+)}$ for the extended SRI between 3.205 and 3.290 ppm with the 6 model orders $K=2500\left( 500 \right) 5000$, color coded as previously, with added noise of two levels differing by a factor of 10. Namely, the added noise is with standard deviations $\sigma =0.289\,\hbox {RMS}$ and $\sigma =2.89\,\hbox {RMS}$, where the latter level is 100 times larger than $\sigma =0.0289\,\hbox {RMS}$ in Figs. 1–8. In panel (a) is the partitioned absorption envelope spectra for $A_{K}^{+}$ with $K=2500\left( 500 \right) 5000$ and the noise level of $\sigma =2.89\,\hbox {RMS}$. The PC resonance is correctly placed and its spread is quite minimal, although the wing to its right side shows substantial spread among the 6 model orders K. Though identifiable for all the model orders, PE is displaced somewhat leftward, and the spread for the 6 model orders is apparent. The four single peaks to the left of PE are all clearly defined, with the distinction among the 6 model orders being noticeable, but fairly minimal. Panel (b) of Fig. 9 shows the partitioned absorption envelope spectra for $A_{K}^{+}$ with the 6 model orders with added noise of $\sigma =\,0.289\,\hbox {RMS}$. In panel (b), the lineshape spread for the 6 model orders is rather small. Next, spectra averaging for the 6 model orders is performed in partitioned envelopes, with the results for the two noise levels displayed in panel (c) of Fig. 9. Therein, ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ for $\sigma =\,0.289\,\hbox {RMS}$ is indicated in red, and for $\sigma =2.89\,\hbox {RMS}$ in blue. Their concordance is remarkably close throughout the SRI. In panel (d), we zoom into the critical frequency window [3.215, 3.225] ppm to compare the average absorption envelope from the other partition, namely ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ for $\sigma = 0.289\,\hbox {RMS}$ and $\sigma =2.89\,\hbox {RMS}$ marked in red and blue, respectively. Therein, while both PC and PE are resolved for the two noise levels, it is the PE chemical shift which is correctly reconstructed and with less discrepancy between the envelopes for $\sigma =0.289\,\hbox {RMS}$ and $\sigma =2.89\,\hbox {RMS}$. Pushed to the right, PC shows slightly closer concordance with its actual chemical shift position at the higher noise level, $\sigma =2.89\,\hbox {RMS}$. Panel (e) of Fig. 9 presents the average complete envelopes ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ for both noise levels, again with the color coding of red for $\sigma = 0.289\,\hbox {RMS}$ and blue for $\sigma = 2.89\,\hbox {RMS}$. In panel (e), the noiseless average complete envelope ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ is drawn in green. It is seen that all the 3 curves entirely coincide and, hence, they appear as a blue single curve, which is the last drawn for $\sigma =2.89\,\hbox {RMS}$. Hence, the remarkable noise-suppression capacity of the ${\mathrm {FPT}}^{(+)}$ also in this more demanding case with 10 and 100 times higher noise levels than in Figs. 1–8. In the absorption mode of the complete ${{\{A}_{K}^{+}+ B_{K}^{+}\}}_{\mathrm {Av}}$ from panel (e), there is no suggestion whatsoever of the two resonances PC and PE. Rather, a single smooth Lorentzian appears. Thus, both PC and PE are symbolized by their open circles in panel (e).

Fig. 9

Partitioned envelopes computed non-parametrically in the ${\mathrm {FPT}}^{(+)}$, using the FID sampled at $N=16384$ for model orders $K=2500\left( 500 \right) 5000$. Panel (a): the partitioned absorptions for $A_{K}^{+}\, \hbox { with }\hbox { noise }\hbox { level}\ \sigma =2.89\,\mathrm { RMS}$. Panel (b): the partitioned absorptions for $A_{K}^{+}\, \hbox { with }\hbox { noise }\hbox { level }\sigma =0.289\,\mathrm { RMS}$. Panels (a) and (b) color coded as in Fig. 1. Panel (c): the average absorptions ${{\{A}_{K}^{+}\}}_{\mathrm {Av}}$ with noise levels $\sigma =2.89\,\mathrm { RMS}$ (blue) and $\sigma =0.289\,\mathrm { RMS}$ (red). Panels (a−c): PC correctly located (filled circles). Panel (d): the partitioned envelopes for ${{\{B}_{K}^{+}\}}_{\mathrm {Av}}$ with noise levels $\sigma =2.89\,\mathrm { RMS}$ (blue) and $\sigma =0.289\,\mathrm { RMS}$ (red); PE correctly located (filled circle). Panel (e): The complete absorptions ${{\{A}_{K}^{+}+\mathrm { }B_{K}^{+}\}}_{\mathrm {Av}}$ with noise levels $\sigma =2.89\,\mathrm { RMS}$ (blue), $\sigma =0.289\,\mathrm { RMS}$ (red) and noiseless ${{\{A}_{K}^{+}+B_{K}^{+}\}}_{\mathrm {Av}}$ (green). The averages computed from the corresponding absorptions obtained for the 6 model orders K. All three curves are indistinguishable. (Color online)

5 Discussion

5.1 General considerations

Major challenges were herein placed upon the non-parametric FPT, applied in the partitioned manner, to identify the potential underlying components in the total shape spectra or envelopes. Firstly, in addition to solely processing the noiseless case, as in our previous study [1], progressively greater levels of noise were added. Secondly, instead of using a single model order, as done earlier [1], some 6 model orders ranging from 2500 to 5000 with a very large interval ($\Delta K=500)$ were examined. Thirdly, the partitioned envelopes for these 6 model orders were averaged to assess the stability of estimations within envelope partitions. In other words, we inquired as to whether there was a reasonable degree of concordance between the average spectra and the sequential envelopes for the 6 model orders? For each of these challenges, it can be said that the partitioned non-parametric ${\mathrm {FPT}}^{(+)}$ demonstrated remarkable robustness. Simply stated, in this controlled setting, the partitioned non-parametric ${\mathrm {FPT}}^{(+)}$ passed all the necessary tests that were particularly focused on the qualitative identification of PC and PE in the sense of “present versus absent splitting of PC and PE”.

As in the previous study with the noiseless case [1], the interference mechanism present in $A_{K}^{+}$ and in $B_{K}^{+}$ for the partitioned absorption envelopes as well as in $C_{K}^{+}$ and in $D_{K}^{+}$ for the partitioned dispersion envelopes was such that in all four of these spectra, PC and PE were clearly resolved in the presence of added noise of standard deviations differing by 10 ($\sigma =0.289\, \mathrm {RMS})$ and 100 $(\sigma =2.89\, \mathrm {RMS})$ from the beginning value ($\sigma =0.0289$ RMS: Figs. 1–8). On the other hand, with the sums $A_{K}^{+} + B_{K}^{+}$ for the real part, $\mathrm {Re}\left( P_{K}^{+}/Q_{K}^{+} \right) ,$ or $C_{K}^{+} + D_{K}^{+}$ for the imaginary part, $\mathrm {Im}\left( P_{K}^{+}/Q_{K}^{+} \right)$, the interference effect caused PC and PE to appear as a single unresolved peak. Once again, the absorption and dispersion envelopes for $A_{K}^{+}+ B_{K}^{+}$ and $C_{K}^{+} + D_{K}^{+}$ were found to be identical to the conventional, non-partitioned total shape envelopes $\mathrm {Re}\left( P_{K}^{+}/Q_{K}^{+}\right)$ and $\mathrm {Im}\left( P_{K}^{+}/Q_{K}^{+} \right)$, respectively. This latter finding is further confirmation of the remarkable stability of the sums from the partitions ($A_{K}^{+} + B_{K}^{+}, {C}_{K}^{+} + D_{K}^{+})$ with respect to the associated conventional non-partitioned reconstruction of the total shape spectra. Overall, neither of the two sums $A_{K}^{+} + B_{K}^{+}$ and $C_{K}^{+} + D_{K}^{+}$ nor their customary non-parametric counterparts $\mathrm {Re}\left( P_{K}^{+}/Q_{K}^{+} \right)$ and $\mathrm {Im}\left( P_{K}^{+}/Q_{K}^{+} \right)$, respectively, resolved PC and PE. This, in turn, underscores the pivotal importance of the partitioning strategy.

Sensitivity to model order K with the very wide interval of 500 between successive values of K was apparent, particularly at the highest noise level of $\sigma =2.89\,\hbox {RMS}$. Nevertheless, the partitioning strategy is seen to be sufficiently hardy so as to still ensure that PC and PE were identified for each of the 6 model orders, despite the noticeable variance among them. Moreover, spectra averaging for the 6 model orders was found to provide patently clear delineation of the two resonances, PC and PE, in the partitioned envelopes from Fig. 9 for noise levels $\sigma \, =\, 0.289\,\hbox {RMS}$ and $\sigma \, =\, 2.89\,\hbox {RMS}$, as was also the case for $\sigma \, =\, 0.0289$ from Figs. 1–8.

It should be recalled that within the FPT, the concept of spectra averaging was first introduced for in vivo encoded MRS time signals [24,25,26,27,28], where its effectiveness in both non-parametric and parametric processing was clearly demonstrated. In those studies with in vivo encoded FIDs, validation of the spectra averaging procedure was achieved by examining elimination of spurious spikes in total shape spectra and by investigating convergence of the spectral parameters from which the component spectra are built. Further confirmation of the effectiveness of spectra averaging for in vivo MRS was provided by scrutinizing the poles and zeros, as the key to stability of response functions of systems to external perturbations [28]. With the present study in the controlled setting using the synthesized FIDs reminiscent of the corresponding data encoded by way of in vitro MRS, we have the possibility to validate spectra averaging by comparison with the known input data. This stringent trial is hereby passed, notably with PC distinctly identified by shape estimation alone, albeit known beforehand to completely underlie PE. These findings in the controlled setting further benchmark spectra averaging to overcome the estimation sensitivity to changes in model order K.

With the results of the present study alongside those exclusively from the noiseless case [1], it can be confidently recommended to first apply the non-parametric partitioned FPT, by explicitly extracting the analytical expressions for the real and imaginary parts of complex spectrum $P_{K}/Q_{K}$ following (9)−(12). This would be the initial “screening” step to find out whether the breast cancer biomarker PC is present or absent. The quantification procedure through the parametric FPT could then be applied to those clinical cases in which PC was advantageously identified stepwise by non-parametric processing applied first. By comparison, the FFT conventionally computed envelopes for in vivo encoded MRS time signals are unable to identify the recognized cancer biomarker, PC. It is mainly for this reason that the efforts of the MR community have shifted to the usage of stronger magnetic fields in attempts to detect breast cancers and distinguish these from benign breast lesions. Such a strategy would not only be extremely costly [46], but its effectiveness has not been demonstrated either [1, 47]. Notably, with higher field scanners, even tCho remained undetected in a number of breast cancers, whereas in several benign breast lesions and in normal breast, tCho was reported as actually being present [47].

5.1.1 Particular importance for molecular imaging of the breast through MRSI

The outlined strategy could be especially useful for MRSI, which is often needed in clinical evaluations regarding breast cancer. Notably, volumetric coverage through MRSI is frequently warranted since a single voxel may not always be sufficiently representative of the status of the imaged breast tissue. Viewing MRI and MRS together brings us to molecular imaging through MRSI. By way of MRSI, the chemical specificity of MRS is combined with the spatial localization techniques provided through MRI to obtain multiple MRS time signals over a volume of interest of scanned tissue [48]. While with single-voxel spectroscopy, which is MRS, three orthogonal slices are selected to encompass a specific volume of interest of about 1.5 to $4\,\hbox {cm}^{2}$, in MRSI for a chosen slice a full spectrum at each point of the selected grid is computed. By processing the data, one can later make images of reconstructed peaks or one can zoom into an image and obtain the local spectra at a specific region of interest. Further details about MRSI can be found in e.g. Ref. [49]. The non-parametric partitioned FPT could be first applied to the voxels in the region of interest. Insofar as PC is detected in any of the voxels, the parametric FPT would be used therein for quantification to extract the metabolite concentrations, their chemical shifts and relaxation times $T_{2}^{*}$.

At this juncture, to acknowledge the fuller scope of the usefulness of envelope partitioning, it is pertinent to recall the obstacles that hamper more vigorous progress of MRSI in the clinic. It is MRSI (multi-voxel), even more than MRS (single-voxel), that is being eagerly awaited by the radiologist for the understandable reasons of volumetric coverage of the scanned tissue by the former diagnostic modality. The radiologist wants first and, if possible, quickly to know the clinically reliable answer to a key question: Is PC present or absent? If PC is present, then the abundance or concentration of this cancer biomarker would be the next sought information on the radiologist’s “to do list”. However, getting this done is itself a formidable task given that from the myriad of encoded FIDs in MRSI, one is faced with thousands of computed spectral envelopes that should all be subjected to accurate quantification.

This challenge is further exacerbated by noise (encoded alongside the signal) which is more abundant with MRSI than is the case for MRS. How to cope with this twofold obstacle: efficient estimation of a huge number of spectra and reasonable noise suppression? It is precisely here that the Padé-based envelope partitioning with spectra averaging comes into play via a “split and conquer” strategy. Just like the radiologist’s stepwise approach, so would we in signal processing apply envelope partitioning first to pinpoint the presence of PC, and then to perform the local quantification to reconstruct the concentration of this cancer biomarker as the next task. Such a fast and accurate procedure within the FPT, would simultaneously surmount both mentioned obstacles. And, therefore, it is indeed MRSI where the concept of Padé-devised envelope partitioning would find its match—where it matters most, as per the radiologist’s “wish list”. Needless to say, this practical methodology would also be advantageous for the patient and, due to its efficiency, to the health care system for the obvious reason of this newly designed MRSI. Overall, we see that it is actually MRSI where spectra partitioning is expected to find its most prominent applications of paramount clinical importance.

For each of the specific clinical issues that will now be reviewed, applying MRSI in this way could be particularly helpful. On the basis of the present study and the previous one on the envelope partitioning sub-topic [1], it can now be fully recommended that this stepwise, multi-faceted Padé-based approach be applied to in vivo MRS and MRSI, aimed at validation on clinical MR scanners for breast cancer diagnostics, and beyond. As elaborated in the next subsection, it can be envisioned that this approach could contribute in many ways to a more individualized approach to breast cancer, via Padé-optimized MRS and MRSI.

5.2 Potential contributions to a more personalized approach to breast cancer

In Ref. [11], we elaborated extensively on the implications of Padé-optimized MRS for a more personalized approach to various aspects of the diagnostics and treatment of breast cancer. This was viewed within the framework of molecular imaging for personalized cancer medicine, PCM. Saliently, it has been stated: “Molecular imaging is rapidly gaining recognition as a tool that has the capacity to improve every facet of cancer care. The growing demands among physicians, patients and society for personalized care are increasing the importance of molecular imaging and shaping the development of biomedical imaging as a whole” [34] (p. 182). Our findings in Ref. [11] focused upon the capability of the parametric FPT to identify and quantify PC, as well as a substantial number of other diagnostically important metabolites. We concluded that in addition to aiding initial detection of breast cancer, this information could also be used to better adapt therapy to the individual patient and to monitor the patient thereafter. Herein, we direct our attention on the partitioning strategy within the non-parametric FPT for visualizing PC as a cancer biomarker. Bolstered by our experience with in vivo MRS and the findings of the present paper, the application of spectra averaging is strongly recommended to reduce the sensitivity to model order K.

5.2.1 Monitoring the presence versus absence of PC reflecting response to neoadjuvant therapy

Neoadjuvant chemotherapy (NCT) is currently used to treat patients with locally advanced or inflammatory breast cancers, and is often chosen for aggressive tumors. With NCT, more direct evaluation of breast cancer treatment effectiveness becomes possible [50]. Response to NCT for breast cancer is most often evaluated according to change in tumor size as manifested on MRI together with the pathological characteristics at subsequent surgery, both of which are relatively late indicators. Regarding NCT, it has been stated: “no imaging-based, early-response biomarker has been suitably validated to become incorporated either as standard of care … or as a routine component of clinical trials” [51] (p. 140). There is some, albeit mixed evidence, that changes in total choline, tCho, as assessed via Fourier-based MRS is reflective of early response to NCT [52,53,54]. Insofar the Padé-based partitioning strategy for identifying the presence of phosphocholine, PC, is validated in vivo, this might contribute to more efficient assessment of early response to NCT, as well as for repeated monitoring during and after the therapy. This strategy could be especially helpful for younger women with aggressive breast cancer, since these patients often have poor outcomes, despite intensive treatment [55].

5.2.2 Detection of PC to better identify hypoxic regions for tailoring radiation therapy

Identifying the most significant regions to improve target definition for radiotherapy can be aided by MR-modalities [56]. Localization of tumor hypoxia is of key importance, since it reflects radiotherapy resistance, driving genomic instability and progression to invasive/metastatic breast cancer [57]. In breast cancer models, tumor hypoxia has been found to be associated with regions of higher tCho as assessed via MRSI. The elevated tCho was primarily associated with elevated PC, as shown by in vitro MRS [58]. These findings suggest a potential role for the Padé-based partitioning strategy to help identify hypoxic regions associated with PC, so as to aid tailoring radiation therapy for patients with breast cancer.

5.2.3 Post-therapeutic surveillance for the presence versus absence of PC

The number of women who have survived breast cancer is rapidly growing as a consequence of early detection together with more effective treatment. The optimal strategies for post-therapeutic imaging surveillance remain to be determined [59]. Since a critical concern for survivors is that breast cancer might recur [56], false positive results can be devastating. On the other hand, insofar as the patient does not receive adequate follow-up, a feeling of abandonment often arises [60], and this can also heighten fear of recurrence. As the patient herself becomes more actively involved in decisions about the surveillance strategy, and is better informed about the results, her sense of control through self-management of aftercare is found to improve with beneficial effects [61, 62]. With this aim in mind, greater confidence is needed about the findings, for the patient as well as for the clinician. An efficient strategy for monitoring the presence versus absence of PC, could well contribute to improving post-therapeutic surveillance strategies to better meet the needs of individual patients who have survived breast cancer.

5.2.4 Detection of PC for surveillance of women at high breast cancer risk

For women at high breast cancer risk, imaging surveillance, especially with the excellent sensitivity of MRI, clearly contributes to timely breast cancer detection [63,64,65,66,67]. Being non-invasive, having no interference with child bearing and few adverse long-term effects [68], intensive imaging surveillance appears to be preferred by women at high breast cancer risk [69]. Women cite the feeling of security and reassurance that breast cancer will be detected early as a key benefit of imaging surveillance [70, 71]. However, false positive findings discourage participation in imaging surveillance [72, 73]. Since the breast is a radiosensitive organ and surveillance for women at high risk should begin at an early age, and possibly with greater frequency, radiation exposure associated with mammography is also of concern, especially for women with BRCA1/2 mutations, Li Fraumeni syndrome or ataxia-telangiectasia, for whom there is heightened vulnerability to radiation damage [64, 74]. Further, mammography is generally less effective for women with a genetic predisposition to breast cancer due to increased mammographic density [65], although mammography sometimes detects early breast cancers missed on MRI [75]. Moreover, BRCA1-associated breast cancers may appear on MRI as benign, with relatively smooth margins, and/or with enhancement kinetics typical of fibroadenomas [64]. Consequently, biopsy has been recommended for MRI-detected fibroadenomas or even cyst-like masses in young patients with BRCA1/2 mutations or otherwise strong family breast cancer history [64]. However, this recommendation can lead to more false-positive findings, and thereafter to further worsening in the specificity of MRI after excisional biopsy [76]. It should be noted that women at high breast cancer risk often turn down invitations to participate in MRI surveillance programs due to fear of being sent for biopsy or other testing [77].

There is a pressing need to improve breast cancer imaging strategies to be better adapted to women’s needs, especially their risk profile [78]. Improved diagnostic accuracy, especially specificity, is vital. As noted, molecular imaging with the FFT for in vivo MRS and MRSI has contributed to somewhat higher specificity, mainly through assessment of total choline, tCho. However, reliance upon tCho assessments does not provide sufficiently trustworthy standards to diagnose a breast lesion as cancerous versus benign [1]. In addition to the possibilities of applying the parametric FPT-MRS and FPT-MRSI to quantify phosphocholine, PC, as a breast cancer marker, as elaborated in detail in Refs. [11, 16, 19,20,21, 45], the results presented in Ref. [1] and in the present study, suggest a potential role of the associated non-parametric estimation of partitioned envelopes for initial detection of PC on in vivo data. This strategy is anticipated to be a particularly helpful contribution to MR-based imaging surveillance via MRSI.

The stringent, multi-faceted testing carried out in the present study indicates that through partitioning, the non-parametric fast Padé transform, FPT, can reliably detect phosphocholine, PC, the otherwise hidden component of MRS and MRSI envelopes or total shape spectra. It is now fully warranted to apply the non-parametric, partitioned FPT to in vivo MRS time signals encoded from the breast. Partitioning should be applied together with spectra averaging to diminish the sensitivity to model order K and suppress noise. The expectation is justified that this approach will contribute in practical improvements for diagnosis and management of breast cancer, concordant with the aims of personalized cancer medicine, PCM.

6 Conclusion and future outlooks

All sciences ultimately deal with signals. They may come from very different sources and, yet, they all have one generic feature in common: being carriers of the unknown causes that produced the detected effects (in the sense of inverse problems that are in ordinary life known as reverse engineering). There is no blueprint for success in unlocking the system’s secrets so as to determine its structure. Nevertheless, the abundant past experience in different fields (resulting in 9 fascinating Nobel prizes) indicates that nuclear magnetic resonance, NMR, is universally the most powerful strategy for revealing the structure of general matter, including the anatomy and function of human tissues. The method of NMR spectroscopy is renamed to magnetic resonance spectroscopy, MRS, in medicine to avoid the patients’ fear that the word “nuclear” might be associated with nuclear radiation (otherwise non-existent in MR phenomena). The present theoretical study deals with times signals typical of those from MRS in cancer diagnostics.

Signal processing is a huge field. Traditionally, for decades now, it has been a stand-alone branch, focused almost exclusively upon a variety of applied research of primary interest to engineering, technologies and industries. Not so long ago, however, researchers in this community began to take a closer look at the other disciplines also dealing with time signals. Eventually, among the many remarkable achievements, they encountered the recent advances in quantum chemistry within the realm of the so-named “quantum-mechanical signal processing” [2, 79,80,81,82,83,84]. Therein, their main attention was caught by signal processors based upon rational functions, notably those from the vast family belonging to the Padé approximant [85,86,87,88,89,90,91,92,93]. This is hardly coincidental. Namely, although not always referred to explicitly by its proper name: the Padé approximant, it occurs that rational functions, through the ratio of two polynomials, are the most frequently used response functions for modeling the response of generic systems to disturbances of any kind. To take a broadly familiar example, typing the words “Padé approximant” on Google, way back in September 1987, would return some, at that time, impressive 2000 items. Three decades later, we note a nearly 100-fold increase to formidable 190000 items in September 2017.

Thus, if the Padé approximant has become omnipresent, with its road so remarkably often traveled, which novelties of such modeling under the umbrella of quantum-mechanical signal processing could merit the heightened attention of investigators from the standard signal processing arena? It is, in fact, the judicious and practical combination of new insights in signal processing provided by basic science, and a new prospect for advanced applications of signal processing including those in medical diagnostics. Basic science, through its most prestigious proponent, quantum physics, improved the fundamental notions, concepts and outlooks of standard signal processing, and elevated its status to that of a complete theory. This was done by reformulating spectral analysis (quantification), which is the main mathematical tool of standard signal processing, in the powerful and general language of quantum-mechanical stationary (time-independent) and non-stationary (time-dependent) formalisms, advantageously coupled to the Schrödinger eigenvalue problems, with time signals equivalently conceived as quantum auto-correlations functions. This approach, in turn, stimulated a vigorous reformulation of signal processing in biomedicine by placing rational polynomials as the response functions of choice, front-and-center, where they are needed most in early cancer diagnostics by MRS.

The reason for the versatile applicability of rational polynomials in the role of response functions is in their unprecedentedly useful mathematical expression which, from the onset, provides the sought polar representation of the spectrum. Single polynomials, such as those in the fast Fourier transform, FFT, are usually suitable for describing smooth, regular functions. However, they are utterly inadequate for modeling functions with various types of discontinuities, as well as singularities like poles, branch points and branch cuts. For such irregular functions, rational polynomials, as polynomial quotients, e.g. the Padé approximant, are most suitable. It goes without saying that the same polynomial ratios are also adequate for regular, continuous functions. Rational polynomials describe poles, by default. Poles of rational polynomials, as zeros of denominators, automatically yield the peaks in spectra, as the fingerprint of the system’s structure.

In signal processing, the Padé approximant is equivalently called the fast Padé transform, FPT, for a twofold rationale. First, one of its multiple computational codes, called the Euclid algorithm, scales quasi-linearly as ${N({\mathrm {log}}_{2}N)}^{2}$ with the total signal length N . Second, the Padé frequency spectrum can be inverted to retrieve the input time signals to any desired degree of accuracy. The word “transform” is reserved only to those mappings that are able to preserve the full information (despite corruption by e.g. noise) when passing from one to the other equivalent domain (time to frequency and vice versa, in the present context). Thus, renaming the Padé approximant as the fast Padé transform is amply justified. It coheres with the generalized Weierstrass theorem, which asserts that every function can be represented arbitrarily accurately by expanding it in terms of rational polynomials.

As stated, signal processing encompasses two general methodologies, shape and parameter estimations of spectra, with the former and the latter yielding the envelopes and components, respectively. Usually, the two strategies do not mix. However, we show within the FPT that the line separating these two tools need not be so sharp. Rather, the two approaches can be intertwined by means of the so-called partitioned envelopes whereby shape estimations alone can qualitatively uncover the component spectra. This is especially valuable for the abundant cases where e.g. two or more closely packed (overlapped) peaks appear as a single compound lineshape.

Precisely such a challenging situation is addressed in the present illustration for MRS time signals associated with breast cancer. Therein, the FPT-based partitioned total shape spectra clearly separated phosphocholine, PC, a cancer biomarker, from phosphoethanolamine, PE. These two peaks fuse together into a single ideally symmetric peak in the corresponding conventional, non-partitioned envelopes. This line of investigation first started with synthesized noiseless MRS time signals for a single model order of the FPT [1], as an initial proof-of-concept study. The present examination is a feasibility study in a more realistic situation, involving noisy synthesized MRS time signals for a sequence of model orders subsequently accompanied by spectra averaging. The imposed stringent tests and benchmarking of the FPT have successfully been completed in this general setting, as well.

The expounded strategy is a motivation to further apply the partitioned envelopes from the FPT to in vivo MRS times signals encoded from breast and other organs, as well. The idea is to proceed in two successive steps. First, to exhaust the partitioned envelopes, as an additional degree of freedom in non-parametric or shape estimations, determining whether cancer biomarkers are present at all in composite peaks in which they are customarily invisible. Second, if present, the cancer biomarker abundance is to be quantified by the subsequent local quantification via parametric processing in the relevant, narrow frequency window.

Such a fast, stepwise processing is expected to be particularly critical for advances in magnetic resonance spectroscopic imaging, MRSI. In MRSI, its volumetric coverage of the scanned tissue from multiple voxels (as opposed to a single voxel in MRS) leads to thousands of spectra that need to be reconstructed and every saving of computational time as well as human resources to analyze the results would be most welcome. Such savings, not at the expense of accuracy and clinical reliability, should be feasible in the FPT with its shape estimation of partitioned spectra in the pre-quantification stage, followed by parameter estimation via efficient local quantification. This is one of the prospects deemed to be well worth exploring in the near future.