Cardiac and Respiratory Self-Gating in Radial MRI using an Adapted Singular Spectrum Analysis (SSA-FARY)

TL;DR

This paper introduces SSA-FARY, a novel self-gating method for radial MRI that uses singular spectrum analysis to reliably detect cardiac and respiratory motion, improving image reconstruction without external gating devices.

Contribution

The paper presents SSA-FARY, a new self-gating technique based on singular spectrum analysis for radial MRI, enhancing motion detection and image quality without external synchronization.

Findings

SSA-FARY reliably detects cardiac and respiratory signals.

It improves image reconstruction quality using parallel imaging and compressed sensing.

The method is validated with simulations and in-vivo experiments.

Abstract

Cardiac Magnetic Resonance Imaging (MRI) is time-consuming and error-prone. To ease the patient's burden and to increase the efficiency and robustness of cardiac exams, interest in methods based on continuous steady-state acquisition and self-gating has been growing in recent years. Self-gating methods extract the cardiac and respiratory signals from the measurement data and then retrospectively sort the data into cardiac and respiratory phases. Repeated breathholds and synchronization with the heart beat using some external device as required in conventional MRI are then not necessary. In this work, we introduce a novel self-gating method for radially acquired data based on a dimensionality reduction technique for time-series analysis (SSA-FARY). Building on Singular Spectrum Analysis, a zero-padded, time-delayed embedding of the auto-calibration data is analyzed using Principle Component Analysis. We demonstrate the basic functionality of SSA-FARY using numerical simulations and apply it to in-vivo cardiac radial single-slice bSSFP and Simultaneous Multi-Slice radiofrequency-spoiled gradient-echo measurements, as well as to Stack-of-Stars bSSFP measurements. SSA-FARY reliably detects the cardiac and respiratory motion and separates it from noise. We utilize the generated signals for high-dimensional image reconstruction using parallel imaging and compressed sensing with in-plane wavelet and (spatio-)temporal total-variation regularization.

Tables

Table 1. Table 1: Average heart rate f ¯ heart subscript ¯ 𝑓 heart \bar{f}_{\text{heart}} and corresponding standard deviation σ heart subscript 𝜎 heart \sigma_{\text{heart}} , standard deviation of SSA-FARY from ECG trigger σ trig subscript 𝜎 trig \sigma_{\text{trig}} . End-systolic and end-diastolic left-ventricular blood-pool area of a mid-ventricular slice for ECG-CINE and SSA-FARY and corresponding error. ∗ denotes the volunteer with highly erratic respiration, ∗∗ denotes the volunteer who yawned.

Volunteer	${\overline{𝒇}}_{\text{heart}}$ [$\frac{\mathrm{Hz}}{}$]	${𝝈}_{\text{heart}}$ [$\frac{\mathrm{Hz}}{}$]	${𝝈}_{\text{trig}}$ [$\frac{\mathrm{ms}}{}$]	Diastole [$\frac{{\mathrm{mm}}^2}{}$]			Systole [$\frac{{\mathrm{mm}}^2}{}$]
				ECG	SSA-FARY	Err. [$\%$]	ECG	SSA-FARY	Err. [$\%$]
V1	0.910(4)	0.057(3)	23(1)	2254	2222	1.4	1068	1017	4.8
V2	1.380(3)	0.0541(2)	14.7(7)	1536	1647	7.2	573	706	23.2
V3	0.979(2)	0.0245(1)	28(2)	1911	2161	13.1	1072	1168	11.0
V4	1.245(3)	0.045(2)	23(1)	1998	1956	2.1	866	961	11.0
V5	0.955(4)	0.056(3)	19(1)	2348	2380	1.4	1005	1179	17.3
V6∗	1.195(7)	0.103(5)	34(2)	2244	2220	1.1	1351	1407	4.1
V7∗∗	1.08(1)	0.129(7)	27(1)	3398	3314	2.5	2060	1930	6.3

Table 2. Table 2: Variables and values for the numerical simulation.

Variable	Value
$N_c$	$30$
$A$	$3$
${\varphi }_a$	$0$
${\omega }_a$	$\frac{2\pi }{80}$
$B$	$1$
${\varphi }_b$	$0.5$
${\omega }_b$	$\frac{2\pi }{30}$
${\sigma }_{\text{noise}}$	$2$
$C$	$1$
${\varphi }_c$	$0$
${\omega }_c$	$\frac{2\pi }{10}$
$[t_1,t_2]$	$[220,300]$
$D_1$	$0.15$
$D_2$	$1.5$
$\xi$	$3.75\cdot {10}^{-3}$
$W$	$101$
$\Phi$	$2$

Equations

$${\phi }^t=t\cdot {\phi }_0$$

$$n^{t}:=\left(\begin{array}[]{c}{e^{i\varphi^{t}}}\\ {e^{-i\varphi^{t}}}\\ {e^{i\cdot 2\varphi^{t}}}\\ {e^{-i\cdot 2\varphi^{t}}}\\ \vdots\\ {e^{i\cdot N_{H}\varphi^{t}}}\\ {e^{-i\cdot N_{H}\varphi^{t}}}\\ \end{array}\right),$$

$${\mathbf{X}}_{\text{cor}}={\mathbf{X}}_{\text{raw}}-\mathbf{n}({\mathbf{n}}^{†}{\mathbf{X}}_{\text{raw}}),$$

$${\mathbf{X}}^T=\mathbf{U}\mathbf{S}{\mathbf{V}}^H.$$

$$(X^T{)}_t=\sum _{k=1}^{N_t}{U}_t^k(S{V}^H{)}_k,$$

$$(S{V}^H{)}_k:=\sum _{l=1}^{N_c\cdot N_p}{S}_k^l(V^H{)}_l={\lambda }_k(V^H{)}_k.$$

$$\tilde{\mathbf{X}}=\mathcal{Z}\mathbf{X}.$$

$$\mathbf{A}=\mathcal{H}\tilde{\mathbf{X}}$$

$$\mathbf{A}=\mathbf{U}\mathbf{S}{\mathbf{V}}^H,$$

$$(S{V}^H{)}_k:=\sum _{l=1}^{W\cdot N_c\cdot N_p}{S}_k^l(V^H{)}_l={\lambda }_k(V^H{)}_k.$$

$${\tilde{X}}_c^{t+j}=A_t^{cj}=\sum _{k=1}^{N_t}{U}_t^k(S{V}^H{)}_k^{cj},$$

$$\mathbf{U}\mathbf{S}=\mathbf{A}\mathbf{V},$$

$$U_t^k=\frac{1}{{\lambda }_k}\sum _{c=1}^{N_p\cdot N_c}\sum _{j=1}^W{\tilde{X}}_c^{t+j}{V}_{cj}^k.$$

$$\delta f_{\text{B}}=\frac{f_{\text{s}}}{W},$$

$$\underset{\alpha }{\mathrm{argmin}}{|}_{t,k}\sum _{c=1}^{N_p\cdot N_c}\sum _{j=1}^W{\Vert {\tilde{X}}_c^{t+j}-\alpha (S{V}^H{)}_k^{cj}\Vert }^2.$$

$$\mathbf{A}{|}_{W=1}=\mathcal{H}\tilde{\mathbf{X}}=\mathcal{H}\mathcal{Z}\mathbf{X}={\mathbf{X}}^T.$$

$$\mathbf{C}=\mathbf{A}{\mathbf{A}}^H.$$

$$a(t)=A\sin ({\phi }_a+{\omega }_{a}t+{\varphi }_a(t)),$$

$$b(t)=B\sin ({\phi }_b+{\omega }_{b}t+{\varphi }_b(t)),$$

$${\varphi }_a(t)=\Phi \sin (2\pi t/2T),$$

$${\varphi }_b(t)=\Phi \sin (2\pi t/T),$$

$$x_i(t)=\frac{i}{N_c}a(t)+\frac{(N_c+1-i)}{N_c}b(t),1\le i\le N_c.$$

$$h_{\text{spell}}(t)=C\sin ({\phi }_c+{\omega }_{c}t),t_1\le t\le t_2,$$

$$h_{\text{trend}}(t)=D_1{e}^{\xi t}-D_2,$$

$$X_i(t)=x_i(t)+h(t).$$

Full text

Keywords: Self-Gating, MRI, dimensionality reduction, PCA, Singular Spectrum Analysis

1 Introduction

Magnetic Resonance Imaging (MRI) is an intrinsically slow imaging technique, which makes imaging of moving organs particularly challenging. Still, from the early years of MRI, researchers recognized the great chances and implications of monitoring the beating heart without the use of ionizing radiation and with the superior tissue contrast of MRI. Here, the respiratory and cardiac motion pose high demands on the acquisition and reconstruction.

One solution is real-time imaging, which resolves the true dynamics of the heart but is limited in terms of temporal and spatial resolution and restricted to two-dimensional imaging [1, 2, 3, 4]. In clinical practice, pro- or retrospective gating is typically used, which exploits the quasi-periodicity of the respiratory and cardiac motion to compose a single synthetic heartbeat from data acquired during several actual beats. To synchronize data-acquisition with breathing motion, external devices like a respiratory belt or adapted sequences with navigator readouts are commonly used [5, 6]. However, these devices have to be placed and adjusted individually for each patient. Furthermore, the resulting respiratory signal is not always directly correlated to the motion of the heart [7]. Sequences with additional interleaved navigator acquisitions prolong the measurement substantially and complicate the use of steady-state sequences. Otherwise, breath-hold commands can be used to avoid the need for respiratory gating completely, which, however, can be exhausting, time-consuming and not expedient for sick or non-compliant patients and children. For cardiac gating, the standard in clinical practice is the use of an electrocardiogram (ECG) [8], but the ECG signals can experience signal distortion when MRI sequences with fast gradient switching are utilized [9].

To avoid these drawbacks and to gain more flexibility, techniques have been developed to extract cardiac motion from the data itself, which is known as retrospective self-gating [10]. Similar approaches can also be used to extract respiratory motion [11]. A large number of different strategies for cardiac, respiratory or combined self-gating with Cartesian or Non-Cartesian acquisition were proposed in the past, e.g. [10, 12, 13, 14, 15, 16]. Still, the fundamental idea in most approaches is similar: Either a 1D signal is extracted from certain receive channels using a band-pass filter and specific properties of the acquired auto-calibration (AC) data, or a (sliding window) low-spatial high-temporal resolution reconstruction of a specific region of interest (ROI) is analyzed. A more sophisticated yet simple idea was proposed by Pang et al. [17]: The general concept of dimensionality reduction [18] is applied to the AC data by using a Principle Component Analysis (PCA) to extract the required motion signals. However, the resulting signals are often spoiled by noise or trajectory-dependent oscillations, which makes additional filtering necessary [19, 20, 21]. Moreover, cardiac and respiratory motion are not always clearly separated [22], which complicates data binning into the respective breathing and heart phases and requires the use of further post-processing steps such as coil clustering [23].

To overcome these limitations, we propose the use of an adapted Singular Spectrum Analysis (SSA), which can be thought of as a temporally localized PCA or equivalently as a PCA applied to time-delay embedded coordinates. SSA is an application of the general Karhunen-Loève theorem [24] and a powerful tool for the analysis of dynamical systems, incorporating elements of classical time-series analysis, multivariate statistics, multivariate geometry and signal processing [25]. Broomhead and King derived SSA from Takens’ theorem for the analysis of chaotic dynamical systems, and applied it to the problems of dynamical systems theory [26, 27]. Further development was promoted by Vautard et al. [28, 29]. Since the birth of SSA in 1986 [27, 30] it has found wide-spread application in various fields [31, 32, 33, 34, 35, 36]. SSA can be used for noise reduction, detrending and the identification of oscillatory components [29], hence it is ideally suited for the extraction of vital motion signals such as respiratory and cardiac motion in self-gated MRI.

Nevertheless, conventional univariate SSA can only be applied to single-channel time series, whereas in parallel MRI multiple receive channels (phased array coils) are available. Channels located closer to the heart tend to capture cardiac motion, while coils placed near the diaphragm rather monitor respiratory motion. Manual coil selection [37] can enable the use of univariate SSA for vital motion extraction, but correlated information from other coils is then lost. Moreover, for routine clinical use a fully automated technique is preferred.

Fortunately, univariate SSA has a natural extension for the analysis of a multi-channel time series [30]. However, this multivariate SSA is not a dimensionality reduction technique, but recovers the specific oscillations for each channel rather than to extract a single signal that describes the temporal evolution of the principle motion components.

Here, we adapt the Singular Spectrum Analysis For Advanced Reduction of dimensionalitY, which we dub SSA-FARY [38, 39]. In its original form (multivariate) SSA consists of four steps [25]: I) Hankelization, II) Decomposition, III) Grouping, IV) Backprojection. In SSA-FARY, we remove steps III) and IV) and instead perform a zero-padding operation at the start. We will demonstrate the basic functionality of SSA-FARY in numerical simulations and show reconstructions of in vivo cardiac measurements acquired with single-slice bSSFP and Simultaneous Multi-Slice (SMS) radio-frequency (RF)-spoiled gradient-echo (FLASH) sequences, as well as with a Stack-of-Stars (SOS) bSSFP sequence.

2 Theory

In radial single-slice, SMS or Stack-of-Stars imaging the central k-space point or the central line along the slice-dimension $k_{z}$ ($k_{x}=0$, $k_{y}=0$), respectively, have proved to be ideally suited for self-gating [10, 20]. We extract this AC data from the measurement data and stack all coils and partitions into a single dimension. This yields a multi-channel time-series ${X}_{c}^{\;t}$ of size $[(N_{p}\times N_{c})\times N_{t}]$, with $1\leq t\leq N_{t}$ and $1\leq c\leq(N_{p}\cdot N_{c})$, which contains information about the respiratory and cardiac motion. $N_{t}$ is the total number of central k-space points or lines used for auto-calibration, $N_{p}$ is the number of partitions and $N_{c}$ is the number of receive coils. Each channel $c$ is normalized to have zero mean.

2.1 Correction of the AC Data

System imperfections such as gradient delays and off-resonances usually cause a corruption of the AC data $\bm{X}$, which manifests an oscillation of a trajectory-dependent frequency in radial imaging [19]. This signal fluctuation is often misinterpreted by dimensionality reduction methods as a major signal contribution. This contribution can mostly be removed by a simple orthogonal projection based on its known frequency. Here, we extend this method to also include higher-order harmonics which yields a method that almost completely removes the unwanted signal. For simplicity, we assume a golden angle acquisition scheme. Let $\varphi_{0}$ be the incremental projection angle, then

[TABLE]

is the projection angle used for the acquisition at time step $t$. We define the vector

[TABLE]

containing the oscillations up to the $N_{H}$-th harmonic as a basis for the perturbing oscillation and constrain $\bm{X}$ to be orthogonal to $\bm{n}$,

[TABLE]

with † denoting the pseudo-inverse. This procedure cleans the corrupted signal $\bm{X}_{\text{raw}}$ and yields a corrected time series $\bm{X}_{\text{cor}}$. We use this AC correction method in all presented in vivo experiments.

2.2 Dimensionality Reduction Methods

Principle Component Analysis

PCA can be understood as the rotation of the original coordinate system to a new one with orthogonal axes that coincide with the directions of maximum variable variance [40]. The PCA of a time series $\bm{X}$ can be performed using the Singular Value Decomposition (SVD).

[TABLE]

Here, the diagonal matrix $\bm{S}$ contains the real eigenvalues $\lambda_{1}\geq\dots\geq\lambda_{{N}_{t}}\geq 0$ in decreasing order of magnitude. PCA provides the expansion of $\bm{X}^{T}$ onto the orthonormal [$N_{t}\times N_{t}$] basis $\bm{U}$,

[TABLE]

where the principle components $(SV^{H})_{k}$ are given by

[TABLE]

Since the cardiac and respiratory motion signals contribute as main sources of variation to the time series $\bm{X}$, their temporal behavior should be captured by one of the first basis vectors $U^{k}$, respectively [17].

Singular Spectrum Analysis For Advanced Dimensionality

Reduction (SSA-FARY)

A schematic of the SSA-FARY procedure is depicted in Fig. 1a. In contrast to conventional (multivariate) SSA, we first zero-pad ($\mathcal{Z}$) the second dimension of the AC data $\bm{X}$ to obtain matrix $\tilde{\bm{X}}$ of size [$(N_{p}\times N_{c})\times({N}_{t}+W-1)$],

[TABLE]

Next, we construct a Block-Hankel calibration matrix

[TABLE]

of size [$N_{t}\times((N_{p}\times N_{c})\times W)$]. Here, the Hankelization operator $\mathcal{H}$ slides a window of size [$1\times W$] through channel $\tilde{X}_{c}$ of the zero-padded AC data and takes each block to be a row in the $c$-th column of the calibration matrix. This operation is similar to the construction of the calibration matrix in ESPIRiT [41].

We decompose $\bm{A}$ using an SVD

[TABLE]

and consider $\bm{U}$ of size [${N}_{t}\times{N}_{t}$] as the orthonormal basis that consists of the Empirical Orthonormal Functions (EOFs) $U^{k}$, $1\leq k\leq{N}_{t}$. The principle components $(SV^{H})_{k}$ are given by

[TABLE]

The expansion of $\bm{A}$, or $\tilde{\bm{X}}$, in the basis $\bm{U}$ then reads

[TABLE]

where $1\leq{t}\leq{N}_{t}$ iterates through the temporal samples, $1\leq c\leq(N_{c}\cdot N_{p})$ through the channels and $0\leq j<W$ is the index inside the sliding window. The EOFs can be considered as data-adaptive weighted moving averages of the original time series $\tilde{\bm{X}}$, with $\bm{V}$ being the data-adaptive filters [29, 42],

[TABLE]

[TABLE]

In fact, the columns of $\bm{V}$ can bee seen as a complete eigenfilter decomposition of the original time series [43]. These filters $V^{k}$ act as data-adaptive band-pass filters with a frequency bandwidth $\delta f_{\text{B}}$ given by

[TABLE]

where $f_{\text{s}}$ is the sampling rate [44, 45]. Harris and Yuan showed in [42] for the univariate case that a periodic oscillation contained in the data lead to an even and odd filter. The application of these filters to the original time series constitutes for each oscillation one EOF which is in phase and one which is in quadrature to the original oscillation, respectively.

Vautard et al. [29] proposed another interpretation of the EOFs considering the minimization problem

[TABLE]

The solution of eq. (15) is $\alpha=U_{t}^{\;k}$, thus the EOFs can be obtained by a local least-squares fit of the k-th principle component to the original time series. This locality, determined by the window size $W$, distinguishes SSA-FARY from classical PCA, which does not take the temporal past and future of a sample into account. In fact, PCA is a special case of SSA-FARY with $W=1$,

[TABLE]

In conventional SSA, i.e. where no zero-padding is applied, the EOFs are of reduced length $N_{t}-W+1$. Thus, the exact correspondence in time is lost, which inhibits their use as self-gating signal. In contrast, the EOFs $U^{k}$ in SSA-FARY preserve the length $N_{t}$ of the original time series and can directly be used for self-gating, similar to the eigenvectors $U^{k}$ in PCA. In distinction to PCA, the EOFs in SSA(-FARY) capture temporal oscillations via oscillatory pairs [28]. In particular, if two consecutive eigenvalues are nearly equal, the two corresponding EOFs are nearly periodic with the same period and in quadrature [46], which is a consequence of the filtering property of SSA-FARY [42]. To ensure a proper separation the singular values of different EOF pairs should be distinct, which is called the strong separability condition [47] and usually fulfilled for our application.

These pairs can be seen as the data-adaptive equivalent to the sine-cosine pairs of Fourier analysis [48]. A single EOF-pair might suffice for the analysis of nonlinear and inharmonic oscillations, as it automatically locates intermittent oscillatory regions. In contrast, classical spectral analysis would require a large amount of harmonics or subharmonics of the fundamental period [29, 49].

Comments

The EOF $U^{k}$ is also the k-th left eigenvector of the [${N}_{t}\times{N}_{t}$] real-symmetric cross-covariance matrix

[TABLE]

Depending on the number of acquired spokes and partitions, computing the eigen-decomposition of $\bm{C}$ is usually more efficient than computing the SVD of $\bm{A}$.

In contrast to the parameter-free PCA, for SSA we must define a window size $W$. To long-range correlations in time, $W$ should be large, which - as a trade-off - results in a lower degree of statistical confidence [49]. Vautard et al. [29] showed that SSA can resolve oscillations best when the periods are shorter than the window size $W$. In our study, the window size $W\approx$3\text{,}\frac{\mathrm{s}}{}$$ proofed to be a robust choice for most measurements, independently of the utilized sequence, and this value was chosen as the default. More information on the choice of the window size is given in the Methods and Discussion section.

The fundamental concept behind the use of a temporal window $W$ is Taken’s delay embedding theorem [26], one of the backbones of chaotic dynamical system analysis. Instead of considering each temporal sample individually and isolated from other time points, so called time-delay coordinates are constructed by embedding the samples in a higher-dimensional space with embedding dimension $W$. Consequently, each time point is represented by a time-delay coordinate vector, which comprises not only the sample of the respective time but also its temporal past and future. It is therefore a natural choice to pick an odd value for $W$ in order to incorporate the same amount of past and future information.

Towards the beginning and the end of a time series this embedding can no longer be constructed due to lack of past or future samples, respectively. There are two strategies to overcome this limitation: 1.) Time-delay coordinates are constructed for the central $N_{t}-(W-1)$ samples only, which means that $W-1$ samples would be discarded from future processing. 2.) $(W-1)/2$ samples are zero-padded on both ends of the time series, which comes at the expense of increased inaccuracy for the marginal samples of the time series. However, for the second approach no samples have to be discarded and through the symmetric zero-padding the time-delay coordinates remain in sync with the actual temporal evolution of the signal. In this manuscript, the second approach is used.

2.3 Binning

The myocardium shows different behavior for contraction and expansion, so usually the entire cardiac cycle is divided into multiple distinct bins to accurately resolve the temporal motion. For respiratory gating it is usually assumed that inspiration and expiration do not have to be distinguished [50, 20, 16]. However, various studies reveal that respiratory motion is heavily subject-dependent and exhibits a strong variability as well as hysteresis, which affects the global position of the myocardium [51, 52, 53, 54]. Hence, inspiration and expiration should be distinguished to properly resolve the effects of breathing motion on the heart.

Since SSA-FARY yields EOF quadrature pairs that capture the phase information of periodic oscillations, binning is straight-forward for both cardiac and respiratory motion: The phase portrait, i.e. the amplitude-amplitude scatter plot, of an EOF pair is divided into $N$ circular sectors with central angle $\phi=$360\text{,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$/{N}$. The samples are then binned according to their respective circular sector, see Fig. 1b.

3 Methods

3.1 Numerical Simulations

We compare the capability of PCA and SSA-FARY in extracting and separating oscillatory signals in simple numerical simulations.

The signals we want to extract are two frequency-modulated sinusoids

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

which account for frequency variations. To simulate various channels $i$, we use a weighted sum

[TABLE]

To spoil the composite signals $x_{i}(t)$, we add Gaussian white noise $h_{\text{noise}}(t)$ with standard deviation $\sigma_{\text{noise}}$, or an oscillatory spell from time $t_{1}$ to time $t_{2}$,

[TABLE]

or an exponential trend

[TABLE]

to all channels, which yields

[TABLE]

We analyze the time-series $\bm{X}^{\text{noise}}$, $\bm{X}^{\text{spell}}$ and $\bm{X}^{\text{trend}}$ using PCA and SSA-FARY with window size $W$. Note, that the aim of this numerical experiment is to demonstrate the general benefits of SSA-FARY over PCA for the analysis of time series, and not to simulate cardiac and respiratory motion in the most accurate way. We therefore did not include the modeling of a more complex frequency variability or motion signal shapes. More details on the simulation are provided in the appendix.

3.2 In Vivo experiments

All measurements were performed on a Skyra 3T scanner (Siemens Healthcare GmbH, Erlangen, Germany) using 30 channels of a thorax and spine coil. Gradient delay correction was performed using RING [55, 56]. The AC data was corrected using the orthogonal projection with $N_{H}=5$. In the following, the AC data’s real and imaginary part are treated as individual channels. The field of view in all experiments was $256\times$256\text{,}\frac{\mathrm{mm}}{}$^{2}$ at base resolution $192$. All presented experiments were performed on volunteers with no known diseases, who gave written informed consent. All SOS measurements were performed on different volunteers. The study had received approval from the local ethics committee.

Sequence Design, Auto-Calibration and Reconstruction

We utilize a radial bSSFP sequence for the single-slice measurement, an RF-spoiled gradient-echo sequence with randomized RF spoiling [57] for the SMS measurement and a radial bSSFP sequence with undersampling in $k_{z}$ direction [58] for the SOS measurement. To obtain maximum k-space coverage and thus improved image quality [59, 60], the projection angle $\varphi$ is increased in each shot about the seventh tiny golden angle $\varphi_{0}\approx$23.6\text{,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$$ [61]. For the SMS and SOS measurements, the partitions are acquired in an interleaved fashion, i.e. one spoke is recorded for each partition before the next in-plane spoke of a partition is acquired.

For single-slice imaging, the central sample of all spokes is used for auto-calibration. For SMS, more AC data is available as not only a single sample but the central line along the $k_{z}$ ($k_{x}=0,\;k_{y}=0$) direction can be utilized.

Inspired by [58], we make use of variable-density $k_{z}$-undersampling for the SOS acquisition. From the total number of 14 partitions, the central 6 partitions are always acquired and the corresponding central line is used for auto-calibration. The remaining 8 partitions are undersampled by a factor of 4. This center-dense sampling scheme does not only increase the temporal resolution of the AC data but also improves the image quality [62].

For self-gating with PCA and SSA-FARY as well as for imaging reconstruction we use BART [63]. Image reconstruction is performed using combined parallel imaging and compressed sensing (PICS) [64] applying the alternating direction method of multipliers (ADMM) [65] with in-plane wavelet-regularization on the spatial dimensions and total variation (TV) on the cardiac and respiratory dimension [20, 66]. For SOS imaging, we additionally apply TV regularization in slice direction. The coil sensitivities for the single-slice and SMS measurements are generated using radial ENLIVE allowing two maps [3, 60, 67, 41]. To reduce the memory demand we allow only one map in the SOS reconstruction. We apply coil compression [68, 69] to reduce the number of coils to 13 for single-slice and SMS, and 10 for SOS imaging and perform the calibration of the sensitivities using a lower resolution.

The SSA-FARY gating signals is distributed into 25 cardiac and 9 respiratory bins for image reconstruction. Although not always necessary, we standardly perform an additional detrending of the EOFs using a moving average filter of length $L_{\text{avg}}\approx 3W$. This further improves the binning accuracy of SSA-FARY by removing a possibly remaining residual trend.

In the spirit of reproducible research, code and data to reproduce the experiments are made available on Github.111https://github.com/mrirecon/SSA-FARY

Single-slice imaging

We perform a $72$ second free-breathing bSSFP scan (TE/TR = $1.90/$3.80\text{,}\frac{\mathrm{ms}}{}, flip angle $40\text{\,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$) with slice-thickness $7\text{\,}\frac{\mathrm{mm}}{}$ of the human heart in short-axis view and use the first 30 seconds of data for further analysis and image reconstruction. The full $72$ second scan is used for a gridding reconstruction in Supplementary Material chapter II. We furthermore conduct an ECG-triggered CINE bSSFP breath-hold scan of the same slice (TE/TR = $1.52/$3.04\text{\,}\frac{\mathrm{ms}}{}, flip angle $57\text{\,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$, slice-thickness $7\text{\,}\frac{\mathrm{mm}}{}$).

We compare the principal motion signals using PCA and SSA-FARY with window size $W=751$, which covers a period of about $3\text{\,}\frac{\mathrm{s}}{}$.

SMS imaging

We perform a 60 second free-breathing SMS RF-spoiled gradient-echo scan (TE/TR = $1.79/$2.90\text{,}\frac{\mathrm{ms}}{}$$, flip angle $12\text{\,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$) with three simultaneously acquired slices in short-axis view. The slice thickness is $5\text{\,}\frac{\mathrm{mm}}{}$ and the slice gap $10\text{\,}\frac{\mathrm{mm}}{}$. We use SSA-FARY with window size $W=345$, which covers a period of about $3\text{\,}\frac{\mathrm{s}}{}$. We perform a joint reconstruction of all slices using binning based on SSA-FARY. To evaluate the accuracy of SSA-FARY, we compare the SSA-FARY respiration quadrature-signals of the complete time series with the breathing pattern extracted from a pneumatic respiratory belt (Siemens Healthcare GmbH, Erlangen, Germany) and a real-time reconstruction [60, 70].

SOS imaging

We measure 8 volunteers and on each we perform a free-breathing three-minute radial SOS bSSFP scan (TE/TR = $1.90/$3.80\text{,}\frac{\mathrm{ms}}{}$$, flip angle $35\text{\,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$) with fourteen partitions in short-axis view and slice thickness $7\text{\,}\frac{\mathrm{mm}}{}$. By default we use SSA-FARY with window size $W=91$, which covers a period of about $2.8\text{\,}\frac{\mathrm{s}}{}$. For one volunteer the respiratory EOF pair revealed highly irregular breathing with occasional breath-holds of up to $12\text{\,}\frac{\mathrm{s}}{}$ and for another volunteer the cardiac EOF pair showed a pronounced frequency variation. To improve the gating accuracy in these two cases, we determined another respiratory EOF pair using window size $W=51$ and another cardiac EOF pair using window size $W=21$, respectively.

To evaluate the precision of the cardiac gating with SSA-FARY, we analyze the SSA-FARY signal in comparison to the simultaneously acquired ECG trigger using a Python 3 script. Therefore, we define a synthetic trigger point when the phase related to the orthogonal cardiac SSA-FARY quadrature pair experiences a zero-phase crossing. Since the global phase-offset of the quadrature pair is arbitrary, we correct the synthetic SSA-FARY trigger by a constant shift using the average distance to the ECG trigger. We then compute the standard deviation $\sigma_{\text{trig}}$ and standard error of the corrected SSA-FARY trigger to the ECG trigger.

Moreover, we acquire the same slices using a conventional ECG-triggered breath-hold CINE single-slice stack bSSFP measurement with Cartesian read-out (TE/TR = $1.52/$3.04\text{,}\frac{\mathrm{ms}}{}, flip angles depending on the Specific Absorbtion Rate (SAR) limits between $47\text{\,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$ and $62\text{\,}\frac{\mathrm{\SIUnitSymbolDegree}}{}$) and cardiac bin size $T_{\text{bin}}\approx$48\text{\,}\frac{\mathrm{ms}}{}. The patient dependent measurement time is around $6-$8\text{,}\frac{\mathrm{min}}{}$$. We compare the end-diastolic and end-systolic left-ventricular blood-pool area of a mid-ventricular slice to the SSA-FARY based reconstruction using ImageJ.

4 Results

Numerical simulations

Fig. 2 depicts the results of the numerical simulations. PCA is able to extract, at least essentially, the shape of the two oscillations $a(t)$ and $b(t)$ from $\bm{X}^{\text{noise}}$ and $\bm{X}^{\text{spell}}$. However, both the noise and the oscillation spell are still evident in the resulting eigenvectors and corrupt the signal. PCA fails to produce a useful result for $\bm{X}^{\text{trend}}$. While oscillation $a(t)$ is present in two principle eigenvectors, oscillation $b(t)$ is not distinctly separated in any of the eigenvectors.

In contrast, SSA-FARY extracts the oscillation signals with almost no spoiling residuals in all three investigated cases. Only at the borders deviations from the ideal signal can be observed. Notably, SSA-FARY does not only yield a one-dimensional signal of the temporal evolution of an oscillation, but preserves the phase information by quadrature pairs, as can be appreciated from the amplitude-amplitude plots. The two EOFs corresponding to the same pair have very similar singular values. In the scree plot of Fig. 2 b) and c), the first plateau corresponds to EOF 1 and 2, and the second plateau belongs to EOF 3 and 4. In Fig. 2 d), the first two plateaus correspond to EOF 1 and 2, and to EOF 4 and 5. Hence, SSA-FARY does not mix the trend into the oscillations but creates an additional EOF to account for it.222In Supplementary Material Fig. 1-3 we provide results for simulations with different noise variance and trend and spell amplitudes, as well as a different frequency variation. More information on these figures is provided in chapter I of the supplementary document.

Single-slice imaging

For each coil, Fig. 3 depicts the DC component of 100 consecutive spokes before and after the data correction using the orthogonal projection. Before correction some coils exhibit pronounced oscillations with approximately 15 samples per period. As we have used the seventh tiny golden angle ($\varphi_{0}\approx 23.6^{\circ}$) the oscillations period of 15 samples corresponds to $\varphi^{15}=15\cdot\varphi_{0}\approx 354^{\circ}$ (see eq. (1)). Hence, the oscillation period in the AC data is linked to the period of the projection angle. By removing this frequency and the higher-order harmonics these oscillations can be completely eliminated.

Fig. 4 shows the self-gating signals generated with PCA and SSA-FARY. In SSA-FARY, the first two EOFs represent respiratory motion and the third and fourth EOF cardiac motion. The EOFs of the pairs are in quadrature, respectively. Both the cardiac and the respiratory phases are well separated. In contrast, PCA cannot fully extract and separate the signals as respiratory and cardiac motion are superposed and heavily spoiled by noise.

Fig. 4c shows six representative images of the SSA-FARY-gated reconstruction. Depicted are from bottom to top the end-systolic, an intermittent and the end-diastolic frame for end-expiration and end-inspiration, respectively. For comparison, we also show the result of the ECG-triggered CINE breath-hold scan.

In Supplementary Material chapter II we present the results of a conventional gridding reconstruction of the full 72 second measurement gated with SSA-FARY.333Supplementary Material Fig. 4 shows the SSA-FARY and PCA gating signals of the 72 second scan. Moreover, a gridding reconstruction for three different cardiac phases is depicted and the results of a breath-hold and an ECG-gated CINE scan are shown for comparison.

In Supplementary Material chapter III we present a similar experiment with a RF-spoiled gradient-echo sequence.444Supplementary Material Fig. 5 shows the effect of the proposed correction on the AC region. Supplementary Material Fig. 6 shows the self-gating signals determined with SSA-FARY and PCA. Furthermore, six representative frames of the SSA-FARY-gated reconstruction for different respiratory and cardiac motion states are depicted.

SMS imaging

Fig. 5a shows 5 out of 25 cardiac phases for all three slices in end-expiration. The different systolic and diastolic phases are well resolved.

The background of Fig. 5b shows the temporal evolution of a line extracted from an SMS-NLINV real-time reconstruction of the full time series. The line was placed in slice two in vertical direction such that the actual motion of the diaphragm can be observed. This diaphragmatic motion is linearly related to the translation of the heart during breathing [71]. On top of the background image, we have plotted the respiratory EOF quadrature pair obtained from SSA-FARY as well as the signal provided by the respiratory belt. For the purpose of comparison, the motion signals were scaled according to the amplitude of the diaphragm motion in the background image.

One of the EOFs (light blue) coincides very well both in amplitude and phase with the temporal evolution of the diaphragm. This is in line with the filtering interpretation of SSA-FARY, which states that one EOF of the quadrature pair is in phase with the related underlying motion, whereas the other EOF (dark blue) constitutes the quadrature signal. The in-phase EOF is furthermore in good agreement with the motion signal provided by the respiratory belt.

The presented volunteer shows a largely periodic breathing pattern. In Supplementary Material chapter IV we present the results of the same experiment on another volunteer, which happend to exhibit a highly irregular respiration pattern during measurement.555 In Supplementary Material Fig. 7 we provide results of the same experiment conducted on a different volunteer. We depict 5 cardiac phases for all three slices, as well as a comparison between the SSA-FARY respiratory self-gating signal, the respiratory belt and the diaphragm motion extracted from a real-time reconstruction.

SOS imaging

Fig. 6a depicts a zoomed view on EOFs representing respiratory and cardiac motion for different window sizes. For window size $W=31$, which - considering the undersampling scheme - covers a period of $\approx$1.0\text{,}\frac{\mathrm{s}}{}, respiratory and cardiac motion are not fully separated and appear superposed in one EOF for respiratory and one EOF for cardiac motion. In contrast, for the proposed window size $W=91$ ($\approx$2.8\text{\,}\frac{\mathrm{s}}{}) the motion signals are well separated. Then again, for $W=151$ ($\approx$4.5\text{,}\frac{\mathrm{s}}{}$$) a signal loss can be observed in the cardiac EOFs. Fig. 6b shows one respiratory EOF and one cardiac EOF for the windows $W=81$, $W=91$ and $W=101$, covering periods from $2.5\text{\,}\frac{\mathrm{s}}{}$ to $3.0\text{\,}\frac{\mathrm{s}}{}$. All signals are in good agreement. The corresponding pairing-components of the EOFs show similar behavior and are therefore not depicted. Fig. 6c presents end-diastolic and end-systolic frames for three out of fourteen slices in end-expiration. All respiratory and cardiac states are well separated, the cardiac wall and the diaphragm are sharply resolved. Note, however, that some slices at the fringe of the slab have low signal intensity due to an unoptimized excitation profile. The image quality is comparable to the CINE reconstruction Fig. 6d, although the latter tends to be sharper. Due to the higher flip-angle, which is restricted by SAR limitations in volumetric sequences, the CINE images possess a better blood myocardium contrast.

For the different volunteers, Table (1) shows the average heart-rate $\bar{f}_{\text{heart}}$, its standard deviation $\sigma_{\text{heart}}$, the standard deviation $\sigma_{\text{trig}}$ of the synthetic SSA-FARY trigger to the ECG trigger and the end-diastolic and end-systolic left-ventricular blood-pool area of a mid-ventricular slice for CINE and SSA-FARY.

For imaging on a 3T system an insufficient shim can lead to banding artifacts. Measurements with bandings affecting the heart could not immediately be noticed and repeated as the reconstruction was performed offline. Therefore, these measurements were discarded for which the analysis of one volunteer is omitted.

The observed average heart rates range from $0.91\text{\,}\frac{\mathrm{Hz}}{}$ ($54.60\text{\,}\frac{}{}$ bpm) to $1.38\text{\,}\frac{\mathrm{Hz}}{}$ ($82.20\text{\,}\frac{}{}$ bpm) with different heart rate variabilities. For the given SOS acquisition with 14 partitions (6 AC lines and undersampling factor of 4), the temporal resolution of the SSA-FARY trigger is $30.4\text{\,}\frac{\mathrm{ms}}{}$. The standard deviation of the SSA-FARY trigger to the ECG trigger, $\sigma_{\text{trig}}$, is of similar size. Hence, the SSA-FARY trigger is in good agreement with the ECG signal and also matches the temporal resolution of the ECG-CINE acquisition.

The areas of the chosen mid-ventricular slices are comparable for ECG-CINE and SSA-FARY, particularly for end-diastole the difference lies mostly in the lower single-digit percent range, whereas a larger uncertainty can be observed for end-systole.

Volunteer V6 exhibits a highly erratic breathing pattern and volunteer V7 possesses a strongly irregular heartbeat and furthermore unintentionally yawned three times during the measurement. Still, SSA-FARY can provide satisfying results as Table 1 and the figures in the Supplementary Material show.666Supplementary Material Fig. 8 and 9 display the SSA-FARY gating signal, representative frames of the SSA-FARY-based image reconstruction and the corresponding CINE reconstructions for volunteers V6 and V7. More information on the figures is provided in chapter V of the supplementary document.

For all in vivo experiments, we have attached representative movies as Supplementary Material.777The files Mov1-14 of the Supplementary Material show representative movies of all in vivo reconstructions. More information on the movies is provided in chapter VI of the supplementary document.

5 Discussion

We introduced a novel dimensionality reduction method dubbed ”SSA-FARY”, which is based on Singular Spectrum Analysis and showed that the proposed technique can successfully recover the cardiac and respiratory signal from the AC data of single-slice, SMS and SOS MRI measurements. Moreover, we have proposed an extended orthogonal projection to correct for system imperfections in the AC data.

AC correction

The reasons for the oscillations in the AC data are many fold and according to our experience cannot be eliminated completely by techniques that correct for trajectory errors only [72, 73]. Particularly for bSSFP sequences, eddy current related dephasing of spins additionally compromises the data [74, 75], although this effect should be rather small for the tiny golden angle [61].

Zhang et al. [76] find this to be particularly problematic for 3D bSSFP imaging on a 3T system and propose to eliminate these measurement errors by averaging over the central five samples of each spoke. However, we found this to produce even more oscillations in the AC data, which requires additional filtering. In contrast, the orthogonal projection used here can remove most of the signal perturbations directly. For constant increments of the projection angle, the correction can be thought of as a set of sharp band-stop or notch filters corresponding to higher-order harmonics of a base frequency, which can be calculated from the increment of the projection angle.

Strictly speaking, the AC correction is not mandatory for SSA-FARY as without filtering these spurious oscillations appear as additional EOFs. However, as these oscillations usually manifest distinctly in the AC region, the corresponding EOFs often possess the highest singular values. Consequently, they dominate the result of eq. (15) which reduces the significance and accuracy of the components of interest, i.e. cardiac and respiratory motion, and complicates the analysis. This can be easily avoided by using the proposed orthogonal projection, which corrects for first-order system imperfections.

Alternatively, advanced techniques to measure the trajectory error [77, 78] or higher-order system imperfection corrections could be utilized to account for the oscillations in the AC region [79, 80]. Still, these approaches require additional hardware and/or sequence modifications, which limits their accessibility.

SSA-FARY

The presented dimensionality reduction technique for time-series SSA-FARY can be considered a PCA applied to a time-delayed embedding of the AC data to exploit the locally low-rankness of dynamic time series. The Block-Hankel matrix $\bm{A}$ (eq. (8)) consists of shifted segments of the original time series. Its covariance matrix $C_{t}^{\;t^{\prime}}=A_{t}(A^{H})^{t^{\prime}}$ (eq. (17)) is an array of scalar products specifying the correlation of all pairs of multi-channel segments in the embedding space. The considerable redundancy in the correlations of the segments in the presence of (quasi) repetitive oscillations, e.g. cardiac and respiratory motion, causes $\bm{C}$ to have low-rank. By exploiting not only spatial but also these temporal correlations, SSA-FARY performs better in recognizing temporal patterns than classical PCA and allows the separation of trend, oscillations and noise from the signal. This was successfully demonstrated on numerical simulations of superposed and spoiled sinusoidal time series. In the actual in vivo measurements the respiratory and cardiac motion could be detected and clearly separated for all investigated sequence types.

To yield comparable results, methods like classical PCA must be combined with various pre- and post-processing techniques such as coil-selection or coil-clustering [23], (iterative) band-pass filtering [50] and signal smoothing [20], which - especially for small AC regions - may be unstable and demand further manual tuning. By contrast, in SSA-FARY these steps are implicitly integrated and therefore surplus to requirement. Moreover, since SSA-FARY preserves phase information by producing quadrature pairs it allows for direct binning, which also renders the otherwise mandatory peak detection obsolete [10].

We demonstrated the high quality of the SSA-FARY gating signals by comparison with a real-time reconstruction, a respiratory belt and ECG triggers. Moreover, at free breathing and at a considerably lower acquisition time SSA-FARY achieved a reconstruction quality that comes close to the results of the ECG-triggered CINE breath-hold scans. The analyzed left-ventricular blood-pool area of SSA-FARY reconstructions mostly corresponds well to the CINE results for end-diastole, whereas larger uncertainties were found for end-systole. Here, because the left-ventricular blood-pool area is significantly decreased compared to end-diastolic states, any deviation will result in larger relative errors. One reason might be that due to the high acceleration factor the short end-systolic phase is not perfectly resolved. In addition, some discrepancy between a breath-hold scan and a gated free-breathing scan is expected, since it cannot be guaranteed that the selected respiratory bins exactly matches the anatomical state of the breath-hold scan.

As the aim of this manuscript was the introduction of the self-gating technique, we did not fully optimize the sequence and reconstruction parameters, particularly the number of cardiac and respiratory bins, the regularization values of the ADMM and the undersampling scheme of the SOS sequence. The parameter tuning and the setup of a clinically applicable protocol is left for future investigations.

In single-slice imaging we only use a single sample per time step for auto-calibration and despite the proportionally large window size, which reduces statistical significance, SSA-FARY yields reliable results. Still, SSA-FARY tends to be more resilient when multiple partitions are used and thus more AC data is available, as in SMS or SOS experiments, or when more samples relative to the window size are used for auto-calibration. For SMS or SOS imaging a bSSFP sequence is recommended since for RF-spoiled gradient-echo imaging a loss of contrast in systolic phases can occur due to pre-saturated blood flowing in from other slices, see Fig. 5. Note, however, that bSSFP sequences suffer from SAR limitations due to the increased flip-angle of the RF pulse and are prone to banding artifacts when no proper shimming is conducted.

Due to the zero-padding operation, eq. (7), and the subsequent Hankelization, eq. (8), the first and last samples in the SSA-FARY EOFs suffer from slight approximation errors. Still, for all presented experiments and analysis we did not discard these samples.

Occasionally, the trend was not completely separated from the EOFs for cardiac and respiratory motion, which we fixed by standardly using a moving average filter. The reason for the incomplete separation of trend can be understood by considering eq. (14), which relates the window size $W$ to the frequency bandwidth $\delta f_{\text{B}}$ of the eigenfilters $\bm{V}$ generating the EOFs $\bm{U}$. Our default choice of window size and sampling rate corresponds to $\delta f_{\text{B}}\approx$0.35\text{,}\frac{\mathrm{Hz}}{} for all measurements and sequences. If we choose $W$ too small, the EOFs capture a wider range of frequencies, which can result in a mixing of trend and oscillations. Similarly, if the respiration frequency (usually $f_{\text{resp}}\approx$0.3\text{\,}\frac{\mathrm{Hz}}{}) and cardiac frequency (usually $f_{\text{resp}}\approx$1.0\text{,}\frac{\mathrm{Hz}}{}) happen to be spectrally close to one another, a very small window and the corresponding large frequency pass-band of the filters can hinder a proper separation, as it is the case in Fig. ([6](#S4.F6)a), $W=31$). Then again, if we choose a very large $W$, the frequency response $f_{B}$ of the filters can be too narrow. Thus, possible frequency variations in the cardiac and respiratory motion can no longer be captured by a single (band-limited) EOF, which results in signal voids and/or the generation of additional EOFs representing higher order harmonics. Detecting these related EOFs corresponds to the ’Grouping’ step of conventional SSA [[81](#bib.bib81)]. In the cardiac signal of Fig. ([6](#S4.F6)a) such a signal void caused by a cardiac frequency variation of $\approx$0.3\text{\,}\frac{\mathrm{Hz}}{} can be observed for $W=151$, whereas the signal can still be adequately captured with $W=91$. As a summary, we found $\delta f_{\text{B}}\approx$0.35\text{,}\frac{\mathrm{Hz}}{}$$ to be a robust choice and we propose to choose the window size accordingly using eq. (14).

The computational bottle-neck of SSA-FARY is the SVD of $\bm{C}$ of size $[N_{t}\times N_{t}]$. Especially for single-slice imaging, the window size $W$ required to obtain $\delta f_{\text{B}}\approx$0.35\text{,}\frac{\mathrm{Hz}}{}$$ and the corresponding number of AC samples $N_{t}$ to obtain good results is relatively large. Hence, the decomposition of $\bm{C}$ is rather time-consuming. Nevertheless, there are approaches to significantly speed up the decomposition stage [82, 83].

SSA-FARY reliably detected the EOF pairs corresponding to cardiac and respiratory motion, which for RF-spoiled gradient echo measurements usually possess the highest singular values. It is, however, not determined that the two largest components belong to the cardiac or the respiratory motion. For bSSFP sequences, which are generally more prone to system imperfections, we frequently found other components such as trends to have high singular values, too. Although for this study the components used for gating were chosen by visual inspection of the EOFs, it is fairly simple to automatize the assignment using a frequency analysis.

If one finds the separation of SSA-FARY to work insufficiently, a variation of the window size usually helps to recover a suitable gating result. Still, for patients with strong cardiac or respiratory frequency variations, highly non-periodic respiratory motion or arrhythmia the results of the presented SSA-FARY method might still be insufficient. In this case various promising extensions of SSA(-FARY) exist to improve the results, e.g. ’Nonlinear Laplacian Spectral Analysis’, ’Sliding SSA’, ’Oblique SSA’, ’Nested SSA’ and their combinations [84, 85, 47].

Outlook

The property of SSA-FARY to not only reliably separate cardiac and respiratory motion but also to extract the trend in data, suggests further applications. In dynamic contrast-enhanced MRI, SSA-FARY could replace spline-fitting [20] for separating motion and contrast enhancement. Furthermore, preliminary results suggest that the T1 decay in inversion recovery sequences is detected as an individual component and separated from the motion signals, which would enable the simultaneous reconstruction of parameter maps and self-gated anatomical motion [86].

Last but not least we want to mention that the zero-padding approach used in SSA-FARY is not limited to self-gated MRI but turns multi-variate SSA into a general dimensionality reduction method, which might be useful for many other problems related to time-series analysis.

6 Conclusion

We have introduced a novel SSA-based dimensionality reduction method called SSA-FARY. Its intuitive approach, the easy implementation and its capability to separate cardiac and respiratory motion as well as trend makes it a promising approach for MRI self-gating, particularly when it is combined with efficient data acquisition schemes and a state-of-the-art image reconstruction technique.

7 Acknowledgement

The orthogonality projection to correct the AC data was originally conceived by Dr. Anthony G. Christodoulou (Biomedical Imaging Research Institute, Cedars-Sinai Medical Center). We extended it by including the higher-order harmonics.

Supported by the DZHK (German Centre for Cardiovascular Research), the Physics-to-Medicine Initiative Göttingen (LM der Niedersächsischen Vorab), the DFG (German Research Foundation) under grant UE 189/1-1, and the Studienstiftung des deutschen Volkes. We gratefully acknowledge the support of the NVIDIA Corporation with the donation of one NVIDIA TITAN Xp GPU for this research.

Appendix A Details for Numerical Simulations

Here we present the values for the variables used in the numerical simulation (Fig. 2). The total duration $T=800\;\text{[a.u.]}$ consists of $1000$ discrete samples. The other parameters can be found in Table 2.