Fourier decomposition (FD)¹ MRI provides regional ventilation and perfusion information from a single free-brething acquisition series. In the current FD MRI implementation, the ventilation and perfusion images are produced by integration of frequencies corresponding to respiratory and cardiac cycles in a power spectrum. For proper quantification of FD MRI, a precise estimation of the underlying frequencies is mandatory; however, the peak search in the Fourier spectrum is usually very unstable, inaccurate and slow. Furthermore, the analysis of oscillatory truncated signals using simple fast Fourier transform often results in spectral leakage and Gibbs ringing. In this work we propose to use matrix pencil (MP) method² to get a proper estimation of the respiratory and cardiac frequency components in combination with a linear least squares approach to find the corresponding amplitudes.

As for the contemporary FD, the MP analysis is performed voxel-wise on registered image series i.e. with lung structures aligned along the time axis. MP method has been already suggested as an efficient tool to estimate the spectral components for a sum of $$$M$$$ damped exponentials of form:

$$s_n=\sum^{M}_{p=1}a_{p}e^{(i\omega_p-R_p)}t_n+\eta_n=\sum^{M}_{p=1}a_{p}z^n_p+\eta_n$$

where: $$$s_n$$$- NMR signal in a voxel, $$$n=0...N-1$$$- equidistant data points, $$$\eta_n$$$- noise-term, $$$a_p$$$- complex amplitude, $$$\omega_p$$$- frequency, and $$$R_p$$$-damping factor of $$$M$$$ distinct exponentials with signal poles $$$z_p{\equiv}e^{(i\omega_p-R_p)\Delta{t}}$$$ (see Ref.2).

For MP the pulmonary signal is first summed over the whole image and subsequently low and high-pass filtered to separate cardiac and respiratory cycles. The poles corresponding to the respiratory and cardiac cycles are then retrieved from MP decomposition of low-pass filtered $$$S_r(n)$$$ and high-pass filtered $$$S_c(n)$$$ signals:

$$S_r(n)=\sum^{u}_{p=1}a_{r,p}z^n_{r,p},{\quad}S_c(n)=\sum^{v}_{p=1}a_{c,p}z^n_{c,p}$$

where $$$u$$$ and $$$v$$$ indicate number of detected respiratory and cardiac frequency components. Note that generally, one pole for the low-pass filtered signal will have $$$\omega\sim{0}$$$ and will thus retrieve the local signal offset. From the corresponding respiratory and cardiac poles, two matrices are defined:

$$\textbf{Z}_r=\begin{bmatrix}z_{r,1}&z_{r_2}&{\cdots}&z_{r,u}\\z_{r,1}^2&z_{r_2}^2&{\cdots}&z_{r,u}^2\\{\vdots}&\vdots&{\ddots}&\vdots\\z_{r,1}^N&z_{r_2}^N&{\cdots}&z_{r,u}^N\end{bmatrix},{\quad}\textbf{Z}_c=\begin{bmatrix}z_{c,1}&z_{c_2}&{\cdots}&z_{c,u}\\z_{c,1}^2&z_{c_2}^2&{\cdots}&z_{c,u}^2\\{\vdots}&\vdots&{\ddots}&\vdots\\z_{c,1}^N&z_{c_2}^N&{\cdots}&z_{c,u}^N\end{bmatrix}$$

Subsequently, the complex respiratory and cardiac amplitudes can now be obtained voxel-wise by the Moore-Pensore pseude-inverse of the matrix $$$Z$$$ and the measured local signal time-course $$$s_n(x,y)$$$:

$$\textbf{a}_{r,c}(x,y) = \textbf{Z}^{+}s(x,y)$$

Finally, summation yields the total respiratory and cardiac amplitudes:

$$A_r(x,y)=\sum_{p=1}^{u}a_{r,p}(x,y),{\quad}A_c(x,y)=\sum_{p=1}^{v}a_{c,p}(x,y)$$

MR data acquisition

All measurements were performed on a 1.5T MR-scanner (Avanto, Siemens Healthcare). Four healthy volunteers were scanned during free-breathing with a time-resolved 2D ultra-fast steady-state free precession (ufSSFP) sequence³ using the following parameters: field-of-view=450x450mm², slice thickness=12mm, TE/TR/TA=0.67ms/1.46ms/119ms, flip angle α=65°, bandwidth=2056Hz/pixel, matrix=128x128 (interpolated to 256x256), 200 coronal images, 3.33 images/s, parallel imaging GRAPPA factor=2.

Image post-processing

Respiratory motion in ufSSFP data was compensated with a non-rigid image registration algorithm⁴. Subsequently, FD and MP decomposition of the registered time-series were performed voxel-wise to estimate the respiratory $$$(A_r)$$$ and cardiac $$$(A_c)$$$ amplitudes and create qualitative ventilation- and perfusion-weighted images. Eventually, quantitative analysis of fractional ventilation⁵, perfusion⁶, and blood arrival time⁷ was performed. The analysis was accomplished for time-series of different length $$$N{\in}{140...190}$$$.

Illustrative qualitative ventilation- and perfusion-weighted images obtained in a volunteer using MP and the contemporary FD for time series containing 190 images are shown in Figs.1a-1d. The variability in $$$A_r$$$ and $$$A_c$$$ as a function of time series length measured in segmented lung area (see Figure 1e) is shown in Figure 1f. The difference in maximal and minimal amplitudes $$$A_r$$$ and $$$A_c$$$ can reach up to 32% and 24% with FD, respectively; while for the MP method the difference was 5% and 12% for $$$A_r$$$ and $$$A_c$$$, respectively. Figure 2 shows exemplar quantitative maps of fractional ventilation, perfusion and blood arrival time obtained in a healthy volunteer and generated using either MP or FD for time-series length $$$N=154$$$. Furthermore, maps of differences between the MP and FD quantitative images of fractional ventilation, perfusion and blood arrival time were shown in the bottom row in Figure 2.

In this work, we presented an improved method of regional pulmonary function assessment using contrast-agent free time-resolved ufSSFP imaging. The proposed method employs a spectral analysis based on a matrix pencil approach and linearized least squares fitting. The method is shown to allow for an fully automatic and improved estimation of amplitudes of respiratory and cardiac signal modulation in the lung parenchyma. We have shown that the standard implementation of FD MRI based on the straightforward FFT is due to its increased variability in amplitude estimation not well-suited for quantitative analysis. However, the improved signal analysis based on the MP decomposition method may result in increased accuracy and reproducibility of quantitative parameters derived from the time-resolved chest acquisitions such as regional ventilation, perfusion and blood arrival time.