# Estimating fiber orientation in ultrasound imaging through correlation and multivariate Granger causality

Mahshid Dodel, Adrian Basarab, François Varray

► To cite this version:

Mahshid Dodel, Adrian Basarab, François Varray. Estimating fiber orientation in ultrasound imaging through correlation and multivariate Granger causality. Computers in Biology and Medicine, 2025, 198, pp.111223. 10.1016/j.compbiomed.2025.111223. hal-05373311

HAL Id: hal-05373311

https://hal.science/hal-05373311v1

Submitted on 10 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Estimating Fiber Orientation in Ultrasound Imaging
through Correlation and Multivariate Granger Causality

Mahshid Dodel<sup>a</sup>, Adrian Basarab<sup>a</sup>, François Varray<sup>a</sup>

<sup>a</sup>INSA-Lyon, Université Claude Bernard Lyon 1, Inserm, CNRS, CREATIS UMR 5220, U1294, F-69100 Villeurbanne, France

**Abstract**

Estimating myocardial fiber orientation is essential for understanding the heart's microstructure and diagnosing conditions like arrhythmias or heart failure. While Diffusion Tensor Imaging (DTI) is the current gold standard, it is often unavailable in many clinical settings due to cost and accessibility. Ultrasound offers a promising alternative, especially through methods that analyze the spatial coherence of RF signals, such as Backscatter Tensor Imaging (BTI). However, most existing techniques rely on correlation, which can miss important directional dependencies in the signals.

In this work, we investigate the use of multivariate Granger causality (GC) as a complementary tool for estimating ultrasound-based fiber orientation. We first develop a simulation framework based on vector autoregressive models that allows independent control of correlation and causality. This framework enables systematic evaluation of correlation, time-domain GC, and frequency-domain GC under different signal regimes, including the effect of preprocessing such as bandpass filtering. We show that correlation performs well in correlation-dominant conditions, and GC is more effective when causal dependencies dominate. Frequency-domain GC in particular demonstrated robustness to filtering, whereas time-domain GC was more sensitive to preprocessing.

The methodology was further validated experimentally using an *in vitro* nylon-fiber phantom, imaged with a 1024-element matrix array probe and a Verasonics Vantage 256 system. Coherent plane-wave compounding with 25 tilted transmissions was applied, and fiber orientations were varied from -60° to 60° using a rotation stage. Analysis of the compounded RF data confirmed the simulation findings: correlation provided stable estimates in regions dominated by instantaneous dependencies, while GC better captured directional
patterns and was more sensitive to local structure. Balanced cases, observed both in simulation and experiment, suggested that real myocardial signals may arise from a mixture of correlated and causal interactions.

Finally, fusing correlation and GC outputs through consensus and averag-
ing improved robustness and accuracy of orientation estimation, especially in
noisy or heterogeneous regions. While simulations relied on linear Gaussian
VAR models, the phantom was more homogeneous than real tissue, and only
in-plane orientations were studied, the results provide a foundation for future
extensions toward 3D analysis, *ex vivo* validation, and real-time implemen-
tation. Overall, this study establishes multivariate Granger causality as a
valuable complement to correlation for ultrasound-based fiber orientation es-
timation, demonstrating feasibility and promise as a cost-effective alternative
to DTI with potential impact in cardiac diagnosis and monitoring.

Keywords: Fiber tracking, Ultrasound coherence imaging, Fractional Anisotropy, Granger Causality, Correlation

1. Introduction

Cardiovascular disease remains the leading cause of death in Europe, ac-
counting for over 3 million deaths annually [1]. The burden is particularly
high in middle-income countries, where reductions in mortality have been sig-
nificantly slower than in high-income nations. One contributing factor is the
limited availability of advanced diagnostic tools. Diffusion Tensor Imaging
(DTI) using magnetic resonance imaging is currently the standard method
for estimating how heart muscle fibers are arranged. This information is im-
portant for understanding cardiac function and for detecting problems like
irregular heartbeats or heart failure [2]. However, DTI is expensive, time-
consuming, and often unavailable in resource-limited settings.

Ultrasound imaging offers a practical and accessible alternative. It oper-
ates by transmitting high-frequency sound waves into the body and analyzing
the reflected echoes to reconstruct tissue structures in real time. Given its
broad availability, low cost, and ability to capture dynamic cardiac motion,
ultrasound offers strong potential for the non-invasive estimation of myocar-
dial fiber orientation, aiming to improve diagnostic capabilities and reduce
global health disparities in cardiovascular care [3].

Different ultrasound-based techniques have been developed to estimate
fiber orientation. These include measuring backscatter characteristics [4, 5],
analyzing shear wave speed [6], and computing spatial coherence between
signals [7]. Among these, Backscatter Tensor Imaging (BTI) based on spa-
tial coherence has shown promising results in both synthetic and biological
tissues [8]. However, most of these approaches rely on correlation, which
only captures linear relationships between signals. This can limit their per-
formance in more complex tissues where signal interactions are nonlinear or
directional [9].

Alternative statistical methods have been explored to overcome the limi-
tations of standard correlation-based approaches. Mutual Information is one
measure that can capture both linear and nonlinear dependencies between
signals, making it more suitable for analyzing the complex interactions in bio-
logical tissues. Another powerful technique is Granger Causality (GC), which
is a statistical framework developed by Clive Granger in 1969 [10] to assess
whether one time series can predict another, based on past information. Un-
like direct cause-effect inference, GC identifies directional dependencies—if
the history of variable X_t improves the prediction of variable Y_t beyond what
Y_t's own history can achieve, then X_t is said to "Granger-cause" Y_t. This
approach, often implemented using vector autoregressive (VAR) models [11],
has been widely adopted in disciplines ranging from neuroscience to eco-
nomics [12, 13].

Extensions of GC to multivariate and frequency domains [14, 15] have
enabled more analysis of interactions in complex systems. However, standard
GC relies on the assumptions that signals are linear and stable over time.
In practice, these conditions may not always hold. To address this, new
versions of GC have been developed using tools like machine learning and
sparse modeling to better handle complex and high-dimensional data [16, 12].

In cardiac applications, GC has demonstrated strong potential for ana-
lyzing directional dependencies in cardiac signals. It has been applied in
atrial fibrillation studies to uncover causal interactions between intracar-
diac electrograms and latent sources of electrical activity [17]. It has also
gained attention in biomedical imaging and RF signal processing, where
it enables the detection of directional dependencies that correlation-based
methods may miss. Nonetheless, challenges persist. Signal preprocessing
techniques such as filtering can alter time-series characteristics and distort
GC estimates [18, 19]. Frequency-domain GC partially addresses this issue
by examining interactions across frequency bands [20].

Our recent preliminary work evaluated four statistical methods—Correlation, Mutual Information, Kraskov-based Mutual Information, and GC—for esti-
mating myocardial fiber orientation from ultrasound data [21]. While correlation struggled with complex patterns, we found that MI performed better in simpler regions. GC was most effective in capturing highly anisotropic features, such as sharp directional changes. Applying a filtering strategy based on fractional anisotropy further improved performance. Combining these approaches through fused and consensus mapping provided a robust framework for characterizing fiber orientation, critical for advancing cardiac microstructure assessment [21].

Building on this foundation, our current study focuses on using GC for fiber orientation estimation. To better understand its behavior, we simulate RF signals with controllable levels of causality and correlation. These simulations allow us to assess how preprocessing affects causal inference and to identify conditions where strong correlations may confound results. This controlled analysis supports the reliable application of GC in cardiac imaging.

## 1.1. *Research Objectives*

This study aims to advance myocardial fiber orientation estimation from ultrasound RF data by analyzing and integrating both correlation- and causality-based measures. The specific objectives are:

1. Proposing a framework that integrates correlation- and causality-based measures for improved coherence imaging and fiber orientation mapping.

2. Developing a simulation framework that generates RF signals with controllable levels of correlation and causality to evaluate GC performance under known conditions.

3. Investigating the effect of preprocessing steps, such as filtering, on the reliability of causality and correlation metrics.

By combining simulation-driven analysis with methodological development, this work provides a robust foundation for using GC in cardiac ultrasound imaging and supports its application in real-world diagnostic settings.

The paper is organized as follows. Section 2 introduces Backscatter Tensor Imaging (BTI) and the use of spatial coherence for fiber orientation estimation. Section 3 presents the theoretical foundations of Granger causality in both time and frequency domains. Section 4 describes the proposed simulation framework for generating synthetic RF signals with tunable levels of correlation and causality. Section 5 and Section 6 report and analyze the results from both simulations and experimental data.
## 2. Theoretical Background: Backscatter Tensor Imaging

BTI is a high-resolution ultrasound technique that enables estimation of local tissue fiber orientation by analyzing the spatial coherence of backscattered RF signals. Originally proposed by Papadacci et al. [7], BTI leverages plane wave coherent compounding and spatial correlation analysis to visualize anisotropic structures in biological tissue, such as myocardial fibers.

BTI operates by transmitting multiple tilted plane waves and coherently compounding the received backscattered signals to form focal zones at each voxel [22]. The spatial coherence between pairs of transducer elements within these zones is analyzed to infer local structural organization. The orientation of fibers can then be estimated either by fitting ellipses to the coherence maps or by using Principal Component Analysis (PCA) to extract the dominant direction of anisotropy [8]. An overview of this process is illustrated in Figure 1, which shows the wave transmission, coherence computation, and resulting tractography map.

### 2.1. Probe Geometry and Signal Representation

The ultrasound probe is structured as a two-dimensional matrix of transducer elements, with $N_x/N_y$ elements along the x- and y-axis, respectively. Each element can be indexed as a pair $(i, j)$, where $i \in [1, N_x]$ and $j \in [1, N_y]$.

$$ \begin{pmatrix} (1,1) & (1,2) & \cdots & (1, N_y) \\ (2,1) & (2,2) & \cdots & (2, N_y) \\ \vdots & \vdots & \ddots & \vdots \\ (N_x, 1) & (N_x, 2) & \cdots & (N_x, N_y) \end{pmatrix}. \qquad (1) $$

To facilitate further developments, the spatiotemporal RF data is reshaped into a 2D matrix $S \in \mathbb{R}^{N_t \times (N_x \times N_y)}$, where $N_t$ is the number of time samples and $N_x \times N_y$ is the total number of probe elements.

Figure 2 illustrates the transformation of raw 3D ultrasound RF data into a 2D matrix form suitable for subsequent signal processing.

Each element captures a time-domain RF signal, denoted as $s_n \in \mathbb{R}^{1 \times N_t}$, where $n = 1, 2, \dots, N_x N_y$ represents the flattened index of the element. To construct the data matrix, these temporal signals are concatenated column-wise into a matrix $S \in \mathbb{R}^{N_t \times (N_x \times N_y)}$, where each column corresponds to the time series recorded by a unique spatial location on the probe. The flattening is performed using column-major order, meaning the probe elements are traversed first along the x-axis and then along the y-axis.
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic of the BTI workflow (a) Multiple tilted plane waves are transmitted into the tissue to insonify fiber structures from different angles, and the corresponding backscattered echoes are recorded across the transducer aperture. (b) Coherent compounding is applied to synthetically generate focused signals at each voxel. (c) The 2D spatial coherence function is computed for each voxel to quantify local anisotropy. Fiber angles are then estimated using a PCA-based approach. The resulting fiber angles are assembled into a tractography map, which provides a visual representation of the spatial distribution and continuity of fiber orientation throughout the imaged region.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: The diagram illustrates the probe structure, with arrows indicating how the signals are mapped into the matrix S.</figcaption></figure>

## 2.2. Spatial Coherence Estimation

The normalized signal matrix $\hat{S}$ is computed by subtracting the mean value from each RF signal and dividing by its standard deviation. From this, the sample covariance matrix is calculated:

$$ \Sigma_S = \frac{1}{N_t - 1} \times \hat{S} \times \hat{S}^T \in \mathbb{R}^{(N_x \times N_y) \times (N_x \times N_y)}. \quad (2) $$

To evaluate the coherence between element pairs separated by a spatial shift $(\Delta x, \Delta y)$, we consider $\Delta x \in [-N_x+1, N_x-1]$ and $\Delta y \in [-N_y+1, N_y-1]$. A mask $\Psi_{\Delta x, \Delta y}$ is applied to extract the relevant entries of the covariance matrix $\Sigma_S$, and the spatial coherence function is then given by:

$$ \mathbf{R}(\Delta x, \Delta y) = \frac{1}{|\Psi_{\Delta x, \Delta y}|} \sum_{i} \sum_{j} \Sigma_S(i, j), \quad (3) $$

where $|\Psi_{\Delta x, \Delta y}|$ is the number of valid index pairs for the given shift. The index set $\Psi_{\Delta x, \Delta y}$ is defined as:
$$ \Psi_{\Delta x, \Delta y} = \left\{ (i, j) \left| \begin{array}{l} \left\lfloor \frac{i}{N_y} \right\rfloor - \left\lfloor \frac{j}{N_y} \right\rfloor = \Delta x, \\ (i \text{ mod } N_y) - (j \text{ mod } N_y) = \Delta y \end{array} \right. \right\} \quad (4) $$

In this formulation, each index i or j refers to a specific transducer element in the 2D probe, reshaped into a 1D vector. The ceiling operator [·] and the modulo operator ensure proper mapping between the 1D and 2D representations. For instance, $\Delta x = 1, \Delta y = 0$ corresponds to averaging covariances between vertically adjacent elements.

The resulting coherence map $\mathbf{R}$ has dimensions $(2N_x - 1) \times (2N_y - 1)$ and provides a localized measure of signal similarity across the probe aperture. This coherence metric captures the spatial structure of RF signals and is sensitive to fiber orientation. Examples of coherence maps illustrating this structure can be observed in Figure 1(c).

## 2.3. Inferring Fiber Orientation and Anisotropy

Different strategies can be used to estimate the orientation angle, such as ellipse fitting or PCA [8]. In this work, we use PCA, where the direction of the first principal component is extracted as the estimated fiber angle. In 2D, the coherence function typically forms an elliptical shape reflecting local tissue anisotropy. To quantify this anisotropy, we compute the Fractional Anisotropy (FA), a normalized index ranging from 0 (isotropic) to 1 (highly anisotropic), based on the relative lengths of the short and long axes of the ellipse. A higher FA value indicates a stronger directional preference in the coherence structure. To improve the robustness of orientation estimation, we apply a threshold on FA to exclude points with low anisotropy, as they correspond to nearly isotropic structures, where the coherence function approaches a circular or spherical shape, and thus provide unreliable orientation information [7, 23].

## 3. Granger Causality: Foundations and Methodological Considerations

Time series analysis frequently employs VAR models to study dynamic relationships between multiple time series. It extends the univariate autoregressive model to a multivariate setting, capturing linear dependencies among variables over time [11]. This section introduces VAR models and
their connection to GC, followed by a detailed discussion of pairwise GC in multivariate settings.

### 3.1. Definition of Vector Autoregressive Model

The VAR model extends the autoregressive model to multiple time series [24]. It captures the linear interdependencies between multiple variables by modeling each variable as a function of its past values and the past values of all other variables in the system.

A VAR model is characterized by its order, which refers to the number of past periods (lags) used in the model. The order of the model, denoted as $p$, indicates that the model incorporates information from the last $p$ time steps. A VAR model of order $p$ is denoted as VAR($p$).

Mathematically, a VAR($p$) model for an $N$-dimensional multivariate time series $\mathbf{U}_t = (u_{1,t}, u_{2,t}, \dots, u_{N,t})^\top$ is defined as:

$$ \mathbf{U}_t = \nu + \sum_{i=1}^{p} \mathbf{A}_i \mathbf{U}_{t-i} + \boldsymbol{\varepsilon}_t, \qquad (5) $$

where $\nu$ is an $N \times 1$ vector of intercepts, $\mathbf{A}_i$ are $N \times N$ coefficient matrices, and $\boldsymbol{\varepsilon}_t$ is an $N$-dimensional vector of error terms, assumed to be white noise with zero mean and covariance matrix $\Sigma_\varepsilon$. Each equation in this system models one variable as a function of the past values of all variables, allowing for interactions between them.

### 3.2. Univariate Case: Autoregressive Model

A univariate time series $u_t$ follows an autoregressive model of order $p$ if:

$$ u_t = \mathbf{A}_1 u_{t-1} + \mathbf{A}_2 u_{t-2} + \dots + \mathbf{A}_p u_{t-p} + \boldsymbol{\varepsilon}_t, \qquad (6) $$

where $\boldsymbol{\varepsilon}_t \sim \mathcal{N}(0, \sigma^2)$ is white Gaussian noise. This model does not include interactions with other variables, meaning that GC is not applicable in a univariate setting. Two principal approaches are defining GC: a bivariate and a multivariate (partial) one.

### 3.3. Bivariate Case: Granger Causality in VAR(1)

To analyze interactions between two time series, we consider a bivariate process represented as:

$$ \mathbf{U}_t = \begin{bmatrix} X_t \\ Y_t \end{bmatrix} \in \mathbb{R}^2, \qquad (7) $$
where two variables $X_t$ and $Y_t$ evolve according to a first-order vector autoregressive (VAR(1)) process:

$$ 
\begin{aligned}
X_t &= a_{11}X_{t-1} + a_{12}Y_{t-1} + \varepsilon_{1,t}, \\
Y_t &= a_{21}X_{t-1} + a_{22}Y_{t-1} + \varepsilon_{2,t},
\end{aligned} 
\quad (8) $$

where $a_{ij}$ are scalar coefficients and $\varepsilon_{1,t}$, $\varepsilon_{2,t}$ are zero-mean white noise processes with covariances $\sigma_1^2$ and $\sigma_2^2$.

To determine whether $Y_t$ Granger-causes $X_t$, we compare the following two models:

* **Full model:** Includes both $X_{t-1}$ and $Y_{t-1}$ to predict $X_t$.

* **Restricted model:** Excludes $Y_{t-1}$ by setting $a_{12} = 0$.

For the restricted model:

$$ X_t = a_{11}X_{t-1} + \varepsilon'_{1,t}. \quad (9) $$

If excluding $Y_{t-1}$ significantly increases the residual variance, then $Y_t$ is said to Granger-cause $X_t$. The strength of this causality can be quantified using the Granger Causality Index (GCI):

$$ GCI_{Y \to X} = \ln \left( \frac{\sigma_r^2}{\sigma_f^2} \right), \quad (10) $$

where $\sigma_r^2$ and $\sigma_f^2$ are the variances of residuals from the restricted and full models, respectively. A higher value indicates stronger GC. The same approach applies symmetrically to test $X_t \to Y_t$ by setting $a_{21} = 0$.

This framework establishes the basic methodology for GC: evaluating predictive contributions by comparing residual variances between full and restricted models.

### 3.4. Multivariate Granger Causality (Unconditional)

The multivariate extension of GC follows the same principle as in the bivariate case. In a VAR($p$) model, each variable is regressed on the past values of all variables in the system.

The key difference lies in the presence of additional variables: when testing whether $Y_t$ Granger-causes $X_t$, the analysis must account for the influence of the rest of the system. This is done by comparing:
* A full model that includes lagged terms of all variables (including $Y_t$),

* A restricted model that excludes only the lagged terms of $Y_t$.

GC is established if excluding $Y_t$ significantly increases the residual variance of $X_t$. In this setting, the effect of $Y_t$ is assessed while controlling for all other variables in the system.

### 3.5. Multivariate Pairwise Granger Causality (Conditional)

Conditional GC addresses the limitations of the unconditional method by controlling the influence of other variables, enabling the detection of direct rather than indirect causal effects in multivariate systems.

Let the system be partitioned into three subsets: $X_t$ (target), $Y_t$ (source), and $Z_t$ (conditioning set). The conditional test evaluates whether $Y_t$ provides unique predictive information about $X_t$ that is not already contained in $X_t$'s own past and in $Z_t$.

The full model includes all past influences:

$$X_t = \sum_{k=1}^{p} A_{xx,k}X_{t-k} + \sum_{k=1}^{p} A_{xy,k}Y_{t-k} + \sum_{k=1}^{p} A_{xz,k}Z_{t-k} + \varepsilon_{x,t}, \quad (11)$$

while the restricted model excludes the contribution of $Y_t$:

$$X_t = \sum_{k=1}^{p} A'_{xx,k}X_{t-k} + \sum_{k=1}^{p} A'_{xz,k}Z_{t-k} + \varepsilon'_{x,t}. \quad (12)$$

If removing $Y_t$ significantly increases the residual variance of $X_t$, even after accounting for $Z_t$, then $Y_t$ is said to conditionally Granger-cause $X_t$.

The conditional GCI is defined as:

$$GCI_{Y \to X|Z} = \ln \left( \frac{\det(\Sigma'_{xx})}{\det(\Sigma_{xx})} \right), \quad (13)$$

where $\Sigma'_{xx}$ and $\Sigma_{xx}$ are the residual covariance matrices from the restricted and full models, respectively.
### 3.6. Computation of Time-Domain and Frequency-Domain Granger Causality

GC can be analyzed in both the time and frequency domains. The time-domain formulation, previously described, evaluates whether the past of one variable improves the prediction of another, while the frequency-domain approach decomposes causality across different spectral components. Frequency-domain GC extends the time-domain approach by decomposing causality into distinct frequency components, offering insights into how information flows across specific spectral bands [14].

Let the multivariate process $U_t$ (defined in (7)) be modeled by a stable VAR($p$) process. The transition to the frequency domain begins by rewriting the time-domain VAR model in operator form:

$$A(L)U_t = \varepsilon_t, \qquad (14)$$

where $A(L) = I - \sum_{k=1}^{p} A_k L^k$ is the lag polynomial and $L$ is the lag operator such that $L^k U_t = U_{t-k}$. Solving for $U_t$ gives:

$$U_t = A(L)^{-1} \varepsilon_t. \qquad (15)$$

Applying the Fourier transform, each lag operator $L^k$ is replaced with $e^{-ik\lambda}$, resulting in:

$$U(\lambda) = H(\lambda)\varepsilon(\lambda), \quad \text{with } H(\lambda) = A(\lambda)^{-1}, \qquad (16)$$

with $H(\lambda)$ representing the frequency-domain transfer function and $\varepsilon(\lambda)$ the corresponding innovation vector.

The spectral density matrix of $U(\lambda)$ is:

$$S(\lambda) = H(\lambda)\Sigma_{\varepsilon}H(\lambda)^{*}, \qquad (17)$$

where $\Sigma_{\varepsilon}$ is the covariance matrix of the innovations. We now partition $H(\lambda)$ and $\Sigma_{\varepsilon}$ as:

$$H(\lambda) = \begin{bmatrix} H_{xx}(\lambda) & H_{xy}(\lambda) \\ H_{yx}(\lambda) & H_{yy}(\lambda) \end{bmatrix}, \quad \Sigma_{\varepsilon} = \begin{bmatrix} \Sigma_{xx} & \Sigma_{xy} \\ \Sigma_{yx} & \Sigma_{yy} \end{bmatrix}. \qquad (18)$$

Then, the spectral density of $X_t$ (i.e., the first component of $U_t$) is given by the upper-left entry of $S(\lambda)$:
$$ S_{xx}^{\text{full}}(\lambda) = H_{xx}(\lambda)\Sigma_{xx}H_{xx}^{*}(\lambda) + H_{xy}(\lambda)\Sigma_{yy}H_{xy}^{*}(\lambda) \\ + H_{xx}(\lambda)\Sigma_{xy}H_{xy}^{*}(\lambda) + H_{xy}(\lambda)\Sigma_{yx}H_{yx}^{*}(\lambda). \quad (19) $$

To isolate the influence of $Y_t$ on $X_t$, we define a restricted model in which the lagged terms of $Y_t$ are omitted from the VAR model. This results in a modified transfer function $H^{\text{res}}(\lambda)$, where $H_{xy}(\lambda) = 0$. The restricted spectral density of $X_t$ is therefore:

$$ S_{xx}^{\text{res}}(\lambda) = H_{xx}(\lambda)\Sigma_{xx}H_{xx}^{*}(\lambda). \quad (20) $$

The frequency-domain GC from $Y_t$ to $X_t$ is then computed as:

$$ f_{Y \to X}(\lambda) = \ln \left( \frac{\det S_{xx}^{\text{full}}(\lambda)}{\det S_{xx}^{\text{res}}(\lambda)} \right). \quad (21) $$

This expression quantifies the additional spectral power in $X_t$ that is explained by including $Y_t$. Integrating over the full frequency range yields the corresponding time-domain GC:

$$ GCI_{Y \to X} = \frac{1}{2\pi} \int_{0}^{2\pi} f_{Y \to X}(\lambda) d\lambda. \quad (22) $$

This spectral decomposition enables the analysis of how causal influences vary across different frequency bands, which is essential in applications where oscillatory components play a role, such as EEG analysis and economic cycles [13, 11].

### 3.7. Impact of Filtering on Granger Causality

Filtering is commonly applied in time series analysis for denoising or isolating specific frequency bands. However, it can distort GC estimates by altering signal dependencies and spectral properties. While symmetric filtering does not theoretically affect GC, practical issues such as increased model order, phase distortion, and boundary effects can introduce bias [18].

In the time domain, filtering transforms a time series $U_t$ through convolution:

$$ U'_{t} = \sum_{m=-M}^{M} h_{m} U_{t-m}, \quad (23) $$
which can introduce artificial dependencies or suppress true causal relationships [19].

Filtering in the frequency domain modifies the VAR representation through its transfer function. Starting from the time-domain VAR(p) model defined in (5), and omitting the intercept term $\nu$ (which does not influence GC) the system can be written using the lag operator $L^k U_t = U_{t-k}$ as

$$A(L)U_t = \epsilon_t, \quad A(L) = I - \sum_{k=1}^{p} A_k L^k. \qquad (24)$$

By applying the Fourier transform, the lag operator $L^k$ is replaced with $e^{-i2\pi\lambda^k}$, yielding the frequency-domain transfer function

$$A(\lambda) = I - \sum_{k=1}^{p} A_k e^{-i2\pi\lambda^k}. \qquad (25)$$

Filtering alters $A(\lambda)$, impacting spectral GC calculations and potentially leading to misleading frequency-specific causal interpretations [18]. To minimize these issues, causal filters should be used to avoid information leakage, apply filtering symmetrically to all variables to prevent asymmetric distortions, and compare GC estimates before and after filtering to assess their impact. While filtering remains a valuable tool, it has to be carefully implemented to ensure valid causal inferences.

## 4. Proposed methodology

In the ultrasound experimental data, both directional dependencies (causality) and instantaneous correlations were observed between RF signals [21]. To better understand how these two forms of dependency influence coherence estimation, we leverage the theoretical framework described above to simulate multivariate RF signals where causality and correlation can be independently controlled. This simulation-based approach enables us to evaluate the robustness of GC under realistic conditions, including preprocessing effects such as bandpass filtering.

An overview of the complete framework is provided in Fig. 3. Simulated or experimental RF signals are first coherently compounded. From these signals, either the covariance matrix (2) or the GC matrix (13) is computed. Substituting these matrices into (3) yields the spatial coherence function
$R(\Delta x, \Delta y)$. A complementary map is then constructed by combining information based on correlation and causality. Finally, PCA with FA thresholding is applied to the coherence maps to extract local fiber orientations, enabling tractography visualization.

### 4.1. Controlling Causality and Correlation in VAR Models

VAR models offer a flexible framework for analyzing multivariate time series data, enabling the examination of both causal and correlational relationships. In many applications, controlling and manipulating these relationships is crucial for better understanding the underlying dynamics and evaluating the robustness of causal inference methods. This section discusses the approaches to control causality and correlation in VAR models, supported by foundational literature.

#### 4.1.1. Controlling Causality in VAR Models

In a VAR($p$) model, GC is governed by the coefficient matrices $A_k$ structure. Directed influences between variables are introduced or removed by modifying the off-diagonal elements of these matrices [10, 11].

We consider the bivariate process $U_t$ defined in Equation (7). A general VAR($p$) process for a multivariate system is given by:

$$U_t = \sum_{i=1}^{p} A_i U_{t-i} + \varepsilon_t. \qquad (26)$$

To eliminate GC from variable $Y_t$ to $X_t$, we enforce:

$$A_{xy,k} = 0, \quad \forall k = 1, \dots, p. \qquad (27)$$

Conversely, to introduce a specific causal influence, we set nonzero values for $A_{xy,k}$ while maintaining the stability condition of the VAR model. Stability is typically maintained by verifying that the spectral radius of the companion matrix constructed from the VAR coefficient matrices {$A_1, \dots, A_p$} is less than one [11, 13].

#### 4.1.2. Controlling Correlation in VAR Models

Correlation in VAR models primarily arises from the noise covariance structure $\Sigma$:

$$\Sigma = \mathbb{E}[\varepsilon_t \varepsilon_t^T]. \qquad (28)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Workflow of the proposed framework. Simulated or experimental RF signals are coherently compounded, either the covariance matrix (2) or the Granger causality matrix (13) is computed. By substituting these matrices into (3), a coherence function $R(\Delta x, \Delta y)$ is obtained, and then combined into a complementary representation. PCA with FA thresholding provides local fiber orientation estimates</figcaption></figure>
To control correlation between $X_t$ and $Y_t$ without introducing causality (27), the off-diagonal elements of $\Sigma$ can be adjusted, while keeping $\mathbf{A}_k$ fixed. Specifically, we set:

$$ \Sigma_{xy} = \rho \sqrt{\Sigma_{xx} \Sigma_{yy}}, \qquad (29) $$

where $\Sigma_{xx}$ and $\Sigma_{yy}$ are the variances (diagonal entries of $\Sigma$) corresponding to $X_t$ and $Y_t$, respectively, and $\rho$ is the desired correlation coefficient. By tuning $\rho$, we can impose varying degrees of correlation independently of the directional dependencies encoded in $\mathbf{A}_k$ [11, 25].

### 4.1.3. Trade-offs Between Causality and Correlation

It is essential to distinguish between correlation and causality, as high correlation does not necessarily imply a causal relationship. A well-conditioned VAR model should be designed to ensure that:

1. Causal relationships are controlled via $\mathbf{A}_k$.

2. Correlation structures are manipulated independently via $\Sigma_\varepsilon$.

By carefully adjusting these components, one can generate synthetic data that isolates the effects of causality and correlation, thereby aiding in validating causality estimation methods, such as GC. This framework is particularly useful for assessing the impact of confounding variables, noise interactions, and external influences in time series analysis [14, 13].

Table 1 summarizes the roles of $\mathbf{A}_k$, $\Sigma_\varepsilon$, and $H(\lambda)$ in capturing the system's causal and correlational dynamics.

<table><caption>Table 1: Key components in VAR modeling and their roles.</caption><thead><tr><td>Component</td><td>Role</td><td>Description</td></tr></thead><tbody><tr><td>A<sub>k</sub></td><td>Causality</td><td>Coefficient matrices; off-diagonal elements define directional influences between variables.</td></tr><tr><td>S<sub>ε</sub></td><td>Correlation</td><td>Residual covariance; off-diagonal terms encode instantaneous, non-causal dependencies.</td></tr><tr><td>H(λ)</td><td>Freq. Response</td><td>Transfer function from VAR coefficients; decomposes GC across spectral components.</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 4: Overview of the simulation pipeline used to generate synthetic RF signals with known causality and correlation. The process begins by defining signal parameters and generating a reference coherence map based on a 2D Gaussian model, shaped by eccentricity and orientation. This reference is then expanded into a coherence matrix representing all pairwise relationships. Using the MVGC toolbox, multivariate time series are simulated from VAR models defined by their coefficient matrix, A, representing causality, and their covariance matrix, Σ, representing correlation. The resulting signals are used to analyze the effects of varying causality and correlation.</figcaption></figure>

## 4.2. Simulation Framework

A key challenge in simulating VAR signals for causality analysis lies in generating coefficient matrices that both (i) encode a desired strength of causality and correlation, and (ii) preserve a specific directional orientation. This aspect is particularly important in fiber imaging, where orientation carries physiological significance. To meet these goals, the proposed method begins by constructing a reference coherence map that reflects the intended spatial structure. The process consists of the following components:

**Spatial Coherence Prior Construction.** The reference coherence map is constructed by modeling spatial coherence as a two-dimensional anisotropic Gaussian distribution with a defined orientation. Lag coordinates along the lateral and axial axes are centered around zero to represent spatial delays. The Gaussian variances are scaled globally and shaped by an eccentricity parameter to control anisotropy. A rotation matrix is applied to align the principal axis of the Gaussian with the desired orientation, and the distribution is evaluated over the lag grid. The resulting map is normalized to have a unit value at zero lag for consistency. This continuous coherence map is then discretized into a square matrix representing pairwise coherence be-
tween transducer elements, where each entry corresponds to the expected coherence based on spatial lag.

*Mapping Coherence to VAR Parameters.* The reference coherence matrix guides the construction of the VAR model parameters. The coefficient matrix $A$ is shaped to reflect directional dependencies aligned with the principal axis of the coherence structure. The covariance matrix $\Sigma$ of the innovation process controls baseline correlation between elements.

This setup ensures that GC arises from structured interactions encoded in $A$, while marginal correlation is governed by $\Sigma$, as introduced in (26) and discussed in Sections 3 and 4.1.

*Signal Simulation.* Time series data are simulated using the MVGC tool-box [13], which generates multi-trial, multivariate VAR($p$) data as described in (26). The model takes as input:

• a coefficient matrix set ($A_1, \dots, A_p$), defining the system dynamics;

• a noise covariance matrix $\Sigma$, defining the variance and correlation of innovation terms;

and produces time series driven by independent and identically distributed Gaussian residuals $\varepsilon_t \sim N(0, \Sigma)$.

*Parameter Selection.* The simulation parameters include:

• **Orientation angle $\theta$, eccentricity $e$, and variance scale**, which define the reference coherence map structure (see Section 2);

• **VAR order $p$;**

• **Trial number and number of time points per trial**, for statistical reliability.

*Output and Evaluation.* The simulated signals are processed to evaluate Time-domain GC, frequency-domain GC, and correlation metrics. By systematically varying $A$ and $\Sigma$, we generate scenarios that range from high causality / low correlation to the reverse. This enables the analysis of how signal structure influences causality estimation.

The entire simulation pipeline is illustrated in Figure 4.
### 4.3. Processing the RF Signals

Once synthetic RF signals are generated, using the VAR model defined in (26), the next step is to quantify directional dependencies among signal components using GC analysis, conducted in both the temporal and spectral domains.

In the time domain, as detailed in Section 3.5, GC is computed by comparing the prediction errors of the unrestricted VAR model (11) and the restricted model (12). A statistically significant reduction in prediction error upon including an additional signal indicates a causal relationship, as formalized by the GC index in (13). In the frequency domain, the VAR model is transformed via a Fourier-based decomposition, following Geweke's formulation [14] and the framework outlined in Section 3.6, enabling the detection of frequency-specific directional interactions.

To examine how bandpass filtering affects GC estimation, we applied filtering to the synthetic VAR-generated time series before computing GC in the time or frequency domain. A symmetric FIR bandpass filter with pass-band $[0.7, 1.3] \cdot f_0 / (f_s/2)$ was used, where $f_0$ is the central transmit frequency and $f_s$ is the sampling rate, reflecting the normalized bandwidth of typical ultrasound RF signals.

Our simulations revealed a notable difference in how the two domains respond to filtering: time-domain GC estimates were often attenuated or suppressed, even when the original data encoded strong causal structure. In contrast, frequency-domain GC remained consistent after filtering, suggesting greater robustness to preprocessing effects. These findings highlight that filtering can disrupt the temporal structure necessary for time-domain GC, whereas spectral GC retains directional information distributed across frequency bands.

Both single-trial and multi-trial GC analyses are performed to mimic the structure of experimental data. In the single-trial analysis, GC is computed at a single spatial position, corresponding to one realization of the simulated signal. This reflects a localized estimation and serves as a baseline for variability. In contrast, the multi-trial analysis computes GC across multiple spatial positions, similar to how experimental acquisitions produce RF data at various locations within the target tissue. This enables the estimation of spatially resolved causality, which is essential for reconstructing directional patterns such as fiber orientation.

This distinction allows us to evaluate how causality behaves locally versus across a spatial ensemble, and assess the robustness of fiber direction
estimation under varying conditions.

For each causality or correlation matrix, the principal direction was ex-
tracted using PCA, which was consistent with the approach used in BTI-
based angle estimation. This ensures comparability across methods and fa-
cilitates fiber orientation mapping.

4.4. *Complementary Maps*

To make the most of the strengths of each metric, two combination strategies were developed: the **Fused Map** and the **Consensus Map**. The fused map is created by averaging the normalized outputs from correlation and GC. This integration tends to produce a smoother, more balanced representation of fiber orientation, potentially reflecting a more coherent estimate across metrics.

In contrast, the consensus map aims to emphasize areas of strong agree-
ment. Each metric is first thresholded to retain only the values above 0.6,
producing binary masks. These masks are then added together, and regions
where both correlation and GC agree are highlighted more clearly.

**5. Materials and Results**

5.1. *Data Acquisition and Preprocessing*

To simulate the structure of an oriented muscle sheet, an *in vitro* exper-
iment was set up to replicate the structure of aligned muscle fibers using a
nylon thread. A 0.1 mm nylon thread was wrapped around a disk to form a
synthetic fiber sheet, similar to the method described in [8]. The final sheet
measured 40 × 45 mm, which is larger than the aperture of the 1024-element
matrix array probe used in this study (Vermon, 3.47 MHz, 0.3 mm pitch,
10.5 × 9.6 mm footprint). The disk was mounted on a rotation stage to
control the in-plane fiber orientation angle θ<sub>d</sub>, while keeping the out-of-plane
angle close to zero (ϕ<sub>d</sub> ≈ 0°). Photos of the setup are shown in Fig. 5.

Ultrasound data were acquired using a Verasonics Vantage 256 system.
A total of 25 tilted plane waves were transmitted between -20° and 20°.
Because the system could only handle 256 channels simultaneously, the full
1024-element aperture was divided into four sub-apertures and combined syn-
thetically. Beamforming was done using coherent plane-wave compounding.

To study how orientation affects the signal, acquisitions were performed
at seven different thread angles: θ<sub>d</sub> ∈ {-60°, -45°, -30°, 0°, 30°, 45°, 60°}.
<figure><img src="image_6.png" /><figcaption>Figure 5: Top (a) and front (b) views of the *in vitro* setup for fiber orientation measurements.</figcaption></figure>

Computations for each method were performed at a fixed depth of 37 mm to allow for consistent comparisons.

The coherent plane-wave compounded RF data were subsequently pro-
cessed and analyzed using both correlation and causality-based frameworks.
Specifically, Eq. (2) was applied to compute the covariance matrix for corre-
lation analysis, while Eq. (13) was used to estimate the GC matrix. These
matrices were then inserted into Eq. (3), replacing Σ_s, to obtain the spatial
coherence function R(Δx, Δy). Based on the coherence maps, fiber orienta-
tion was extracted using a PCA-based approach, combined with FA thresh-
olding for robustness, and tractography visualization was performed. The
overall workflow is summarized in Fig. 3. All processing steps were imple-
mented in MATLAB using standard signal-processing routines, together with
the MVGC toolbox [13] for multivariate, time- and frequency-domain GC es-
timation. As ultrasound RF signals are inherently bandpass in nature due
to transducer characteristics, additional filtering was not required.

5.2. *Simulation Results*

This section compares orientation estimation performance across syn-
thetic simulations using correlation, unfiltered single-trial time-domain GC,
and filtered spectral GC methods. Simulated RF signals were generated us-
ing a VAR model following (26). Causal and correlated dependencies were
independently controlled through the coefficient matrix *A* and noise covari-
ance matrix Σ. Seven synthetic configurations were defined, ranging from
purely correlated to purely causal systems, by varying the structure of *A* and
Σ, while maintaining a consistent ground truth orientation (θ) from −75° to
75°.

We applied a banding strategy based on diagonal support to generate
structured coefficient and noise matrices for simulation. Specifically, we used
<figure><img src="image_7.png" /><figcaption>Figure 6: Reference covariance matrices used to generate the coefficient matrix A (top left) and the innovation covariance matrix Σ (bottom left), along with their corresponding banding masks (top and bottom right). The masks are derived using 95% energy-based bandwidth estimation</figcaption></figure>

a tapered estimator of the sample covariance matrix, where each entry is
weighted according to its distance from the diagonal. A binary **banding**
**mask** $W \in \mathbb{R}^{N \times N}$, with $N = N_x \times N_y$, was defined such that $w_{ij} = 1$
if $|i-j| \le k$ and zero otherwise, where $k$ is the **bandwidth parameter**
controlling the number of retained diagonals. This process enforces a spars-
sity structure that reflects local dependencies and is particularly suitable for
modeling spatial coherence in RF signals [26].

In our simulation framework, this approach was used to separately band
two matrices: the VAR coefficient matrix **A**, the noise covariance matrix **Σ**.
The resulting masks W₁ and W₂ were derived using a 95% energy-based cri-
terion, where the smallest bandwidth *k* that captures 95% of the total Frobe-
nius norm was selected. These masks define the effective diagonal support
of **A** and **Σ**, controlling the degree of causality and correlation, respectively,
<figure><img src="image_8.png" /><figcaption>Figure 7: Bandwidth and diagonal support ratios across simulation cases. Left: ratio of estimated bandwidths ($K_1/K_2$). Right: ratio of diagonal spans $(2K_1 + 1)/(2K_2 + 1)$.</figcaption></figure>

introduced into the system.

Fig. 6 illustrates this process for Case 4. The top-left matrix represents the coherence prior used to structure **A**, while the bottom-left matrix shapes **Σ**. The corresponding banding masks $W_1$ and $W_2$ appear on the right. Fig. 7 further summarizes the relative bandwidths across all cases, using two metrics: the ratio of estimated bandwidths $K_1/K_2$, and the ratio of diagonal spans $(2K_1+1)/(2K_2+1)$. These ratios serve as interpretable control parameters for tuning the relative dominance of Granger-causal versus correlation-based dependencies in the simulated signals.

Fig. 8 shows correlation, time-domain GC, and frequency-domain GC maps for Cases 1, 4, and 7 at $\theta = 60^{\circ}$. These configurations represent distinct interaction patterns, characterized by the ratio of Granger-to-correlation support, quantified as the estimated bandwidth ratio $K_1/K_2$.

* **Case 1 (K-ratio = 0.029):** Strong correlation, weak causality; structure primarily in **Σ**.

* **Case 4 (K-ratio = 2):** Balanced contribution from both **A** and **Σ**.

* **Case 7 (K-ratio = 34):** Strong causality, minimal correlation; structure concentrated in **A**.

Fig. 9 presents each method's estimated orientation angles' mean and standard deviation as a function of the ground truth angle $\theta$. Each subplot
<figure><img src="image_9.png" /><figcaption>Figure 8: Directional analysis results for Cases 1, 4, and 7 at $\theta = 60^\circ$, showing (from left to right) correlation, time-domain GC, and frequency-domain GC. The three cases illustrate distinct interaction patterns: Case 1 is correlation-dominant, Case 4 reflects balanced contributions, and Case 7 is causality-dominant.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Comparison of the estimated orientation angles' mean and standard deviation across simulations with varying signal structures and experimental *in vitro* data. Each subplot corresponds to a specific simulation case with increasing Granger-to-correlation dominance, quantified by the estimated bandwidth ratio *K*₁/*K*₂. The final subplot shows experimental results. The dashed line represents perfect agreement with the ground truth angle.</figcaption></figure>

corresponds to a different simulation case (RF 1–7), with the final panel show-
ing results from the experimental *in vitro* data. The dashed line represents
perfect agreement with the ground truth.

In the correlation-dominated cases (RF 1–3), the correlation method
closely matches the ground truth with low variance. Time-domain GC con-
sistently underestimates angles and shows higher variability. Spectral GC
remains near zero in RF 1, but gradually improves in RF 2 and RF 3, though
still less accurate than correlation. These results confirm that correlation is
more reliable than GC methods when the causal structure is weak.

In the balanced case (RF 4, *K*-ratio = 2), all methods achieve simi-
lar performance with small errors and low variance across the full angular
range. This configuration closely mirrors the behavior seen in the experi-
mental dataset, suggesting that the experimental signals may result from a
combination of moderate correlation and strong causality. Notably, in RF 4,
the three methods yield consistent estimates with minimal bias, indicating
mutual reinforcement between correlation and directional information.

GC methods dramatically outperform correlation in the causality-dominated
cases (RF 5–7). The filtered spectral GC and time-domain GC track the
ground truth accurately. Meanwhile, the correlation method breaks down
completely—its estimates are flat, misaligned, or near zero, especially in
RF 6 and RF 7, where the coherence structure is absent from the noise com-
ponent. This confirms that in highly directional systems, causality-based
methods are essential to recover fiber orientation, while correlation fails to
reflect the underlying structure.

5.3. *Experimental Results*

<figure><img src="image_11.png" /><figcaption>Figure 10: Estimated fiber orientation using correlation, GC, and consensus approaches for a ground truth angle of 60°. Top row: angular estimation from one position. Bottom row: spatial fiber maps across 1120 positions.</figcaption></figure>

Fig. 10 presents the orientation estimation results for different methods when the ground truth fiber angle is 60°, with an interpolation factor of 2 applied for visualization clarity.

In the single-position analysis (top row), the correlation estimated an angle of 59.4°, which is close to the reference, but the response is less directional. GC provided a slightly lower estimate of 52.9°, but with a sharper orientation
profile, suggesting better sensitivity to directional structure. The consensus
method yielded an angle of 58.3°, striking a balance between sharpness and
stability by integrating agreement across methods.

The bottom row displays fiber tracking based on the orientation esti-
mated across 1120 spatial positions. GC was able to detect fiber directions
even in edge regions where correlation struggled, indicating its advantage in
capturing causal dependencies in areas with lower coherence. In contrast,
correlation appeared more stable in regions with uniform structure. These
observations suggest that the methods are complementary: GC may perform
better in detecting complex structures (e.g., corners or noisy areas), while
correlation provides robustness in simpler regions.

Finally, the fused and consensus maps combine the strengths of the indi-
vidual methods. These composite maps enhance overall performance by inte-
grating both stability and structural sensitivity, improving coherence analysis
in ultrasound-based fiber orientation estimation.

6. Discussion

Previous studies in ultrasound imaging have mainly relied on correlation-
based techniques, such as Backscatter Tensor Imaging (BTI) [7], coherence-
based anisotropy mapping [23]. While these approaches have shown strong
potential for characterizing myocardial microstructure, they are fundamen-
tally limited to capturing linear relationships between RF signals. Alter-
native strategies include mutual information approaches [9], which extend
beyond linear assumptions but remain less explored in cardiac ultrasound.

In parallel, GC has been successfully applied in other biomedical do-
mains, such as neuroscience for analyzing directional dependencies in EEG
and fMRI [18, 12], and in cardiac electrophysiology for uncovering causal
interactions during atrial fibrillation [17]. These works confirm the value of
causality-based analysis for detecting directional patterns that correlation
alone may miss. To the best of our knowledge, our study is the first to
apply multivariate and frequency-domain GC for myocardial fiber orienta-
tion in ultrasound, thus extending the scope of causality methods into tissue
characterization.

This study evaluated the performance of correlation and GC methods
for estimating fiber orientation using both simulated and *in vitro* ultrasound
data. While the standard approach to fiber orientation estimation is based
on spatial correlation, directional interactions observed in experimental RF
signals suggested the presence of temporal dependencies not captured by
correlation alone. This motivated us to explore whether causality analysis,
particularly GC, could provide complementary or improved information in
these cases. To evaluate this idea under controlled conditions, we developed a
flexible VAR model where the coefficient and noise covariance matrices were
adjusted to generate signals with varying levels of correlation and causality,
from purely correlated to strongly causal.

The results showed that correlation performs reliably in correlation dom-
inant cases, providing stable and accurate angle estimates. However, its
performance decreases in situations where the signal contains strong causal
interactions, as it does not account for temporal dependencies. In contrast,
GC-based methods, especially in the frequency domain, were more effective in
detecting directional patterns in causal regimes. Time-domain GC was more
sensitive to preprocessing, particularly filtering, which reduced its ability to
detect causal influence. All methods produced similar orientation estimates
in the balanced configuration (Case 4), where both correlation and causal-
ity were present. Interestingly, this case also matched most closely with the
experimental results, suggesting that real tissue signals may reflect a com-
bination of spatial and temporal dependencies. This highlights the benefit
of using simulation as a tool to better understand the behavior of different
estimation methods in practical conditions.

We also observed that averaging GC results across multiple spatial posi-
tions improved stability, compared to single-trial analysis. This is consistent
with the way the experimental acquisitions benefit from spatial redundancy.

Several limitations of this study should be acknowledged. First, the sim-
ulations were based on linear VAR models with Gaussian noise, which may
not fully capture the complexity of real biological signals. Moreover, both
correlation and GC assume stationarity, whereas myocardial RF signals can
be influenced by motion and heterogeneity that challenge this assumption.

The computational cost of GC is also higher compared to correlation.
Multivariate and frequency domain GC require autoregressive modeling, mak-
ing them computationally intensive. To evaluate their feasibility, we bench-
marked execution time on a representative RF block. Correlation completed
in approximately 0.06 s and Granger Causality in 1.00 s. These results high-
light the trade-off between speed and robustness, motivating future GPU or
parallel optimization for real-time use. Another limitation concerns the ex-
perimental validation, which was conducted on a synthetic fiber phantom.
While this model is controlled and reproducible, it does not fully replicate
the structural and compositional complexity of myocardial tissue.

Additionally, US RF signals are inherently bandpass due to the charac-
teristics of the US probe. To ensure fair comparisons, we applied similar
filtering on synthetic data. This revealed that filtering can attenuate or
suppress time-domain GC, consistent with previous observations in neuro-
science and electrophysiology studies [18, 19]. Testing different filter designs
confirmed that this sensitivity is not specific to one filter type, but a general
property of causal inference.

The present work focused on in-plane orientation estimation. Extending
these methods to 3D and out-of-plane fiber structures is a necessary step. In
this context, Dumas et al.[8] expanded the 2D BTI framework of Papadacci
et al.[7] to 3D, enabling anisotropy measurement in all directions. To the
best of our knowledge, there is no equivalent 3D extension for GC, which
highlights both the novelty and the challenge of our ongoing work. Our
current efforts aim to extend the GC framework to 3D, inspired by these
approaches, but with the ability to capture both temporal dependencies and
spatial coherence.

Future work will therefore focus on extending the framework to *ex vivo*
cardiac data, exploring alternative causality measures that can capture non-
linear or non-Gaussian effects, and developing fusion strategies that combine
the stability of correlation with the directional sensitivity of GC. Efforts will
also be directed toward improving computational efficiency and finalizing
the 3D GC framework for robust estimation of out-of-plane myocardial fiber
structure.

Beyond methodological development, this work has a clear perspective to-
ward practical application in cardiac imaging. By providing a low-cost and
widely available alternative to diffusion tensor imaging, ultrasound-based co-
herence and causality analysis can enable routine assessment of myocardial
fiber orientation in clinical environments where MRI is not feasible. The
complementary use of correlation for speed and robustness, and GC for di-
rectional sensitivity, offers a balance between computational feasibility and
diagnostic accuracy. With further optimization and validation on *ex vivo* car-
diac data, the proposed framework could support real-time visualization of
fiber orientation during ultrasound examinations. Such capabilities may ulti-
mately assist in improved diagnosis and monitoring of myocardial disorders,
including arrhythmias and heart failure.

Overall, this work highlights how correlation and GC capture different
aspects of signal structure. Their combination can improve fiber orientation
estimation, especially in complex or mixed signal conditions.

## Acknowledgments

This work was supported by the LABEX PRIMES (ANR-11-LABX-0063) of the Université de Lyon within the program "Investissements d'Avenir" (ANR-11-IDEX-0007), operated by the French National Research Agency (ANR), the GDR IASIS, and the ANR DETLA project ANR-24-CE45-7060. This material is based upon work done on the ISO 9001:2015 PILOT facility. The RF Verasonics generator was cofounded by the FEDER program, Saint-Etienne Metropole (SME), and Conseil General de la Loire (CG42) within the framework of the SonoCardio-Protection Project leaded by Pr Pierre Croisille.

## References

[1] A. Timmis, V. Aboyans, P. Vardas, N. Townsend, A. Torbica, M. Kavousi, G. Boriani, R. Huculeci, D. Kazakiewicz, D. Scherr et al., "European society of cardiology: the 2023 atlas of cardiovascular disease statistics," *European Heart Journal*, vol. 45, no. 38, pp. 4019-4062, 2024.

[2] J. Cleutjens, "The infarcted myocardium Simply dead tissue, or a lively target for therapeutic interventions," *Cardiovascular Research*, vol. 44, no. 2, pp. 232-241, 1999.

[3] V. Chan and A. Perlas, "Basics of ultrasound imaging," *Atlas of ultrasound-guided procedures in interventional pain management*, pp. 13-19, 2011.

[4] Y. Song, C. M. Kube, J. Zhang, and X. Li, "Higher-order spatial correlation coefficients of ultrasonic backscattering signals using partial cross-correlation analysis," *The Journal of the Acoustical Society of America*, vol. 147, no. 2, pp. 757-768, 2020.

[5] M. Yang, T. M. Krueger, M. R. Holland, and J. G. Miller, "Anisotropy of the backscatter coefficient of formalin-fixed ovine myocardium," *The Journal of the Acoustical Society of America*, vol. 122, no. 1, pp. 581-586, 2007.
[6] W.-N. Lee, M. Pernot, M. Couade, E. Messas, P. Bruneval, A. Bel, A. A. Hagege, M. Fink, and M. Tanter, "Mapping myocardial fiber orientation using echocardiography-based shear wave imaging," *IEEE transactions on medical imaging*, vol. 31, no. 3, pp. 554-562, 2011.

[7] C. Papadacci, V. Finel, J. Provost, O. Villemain, P. Bruneval, J.-L. Gennisson, M. Tanter, M. Fink, and M. Pernot, "Imaging the dynamics of cardiac fiber orientation in vivo using 3d ultrasound backscatter tensor imaging," *Scientific reports*, vol. 7, no. 1, p. 830, 2017.

[8] R. Dumas, B. Pialot, and F. Varray, "Complete 3d anisotropy measurement with coherence imaging," in *2023 IEEE International Ultrasonics Symposium (IUS)*. IEEE, 2023, pp. 1-4.

[9] A. Papana, "Connectivity analysis for multivariate time series: Correlation vs. causality," *Entropy*, vol. 23, no. 12, p. 1570, 2021.

[10] C. W. Granger, "Investigating causal relations by econometric models and cross-spectral methods," *Econometrica: journal of the Econometric Society*, pp. 424-438, 1969.

[11] H. Lütkepohl, *New introduction to multiple time series analysis*. Springer Science & Business Media, 2005.

[12] A. Shojaie and E. B. Fox, "Granger causality: A review and recent advances," *Annual Review of Statistics and Its Application*, vol. 9, no. 1, pp. 289-319, 2022.

[13] L. Barnett and A. K. Seth, "The mvgc multivariate granger causality toolbox: a new approach to granger-causal inference," *Journal of neuroscience methods*, vol. 223, pp. 50-68, 2014.

[14] J. Geweke, "Measurement of Linear Dependence and Feedback between Multiple Time Series," *Journal of the American Statistical Association*, vol. 77, no. 378, pp. 304-313, Jun. 1982.

[15] M. Eichler, "Granger causality and path diagrams for multivariate time series," *Journal of Econometrics*, vol. 137, no. 2, pp. 334-353, 2007.

[16] S. Basu, A. Shojaie, and G. Michailidis, "Network granger causality with inherent grouping structure," *The Journal of Machine Learning Research*, vol. 16, no. 1, pp. 417-453, 2015.
[17] D. Luengo and V. Elvira, "Latent variable analysis of causal interactions in atrial fibrillation electrograms," in *2016 Computing in Cardiology Conference (CinC)*. IEEE, 2016, pp. 981-984.

[18] L. Barnett and A. K. Seth, "Behaviour of granger causality under filtering: theoretical invariance and practical application," *Journal of neuroscience methods*, vol. 201, no. 2, pp. 404-419, 2011.

[19] E. Florin, J. Gross, J. Pfeifer, G. R. Fink, and L. Timmermann, "The effect of filtering on granger causality based multivariate causality measures," *Neuroimage*, vol. 50, no. 2, pp. 577-588, 2010.

[20] J. F. Geweke, "Measures of conditional linear dependence and feedback between time series," *Journal of the American Statistical Association*, vol. 79, no. 388, pp. 907-915, 1984.

[21] M. Dodel, R. Dumas, A. Basarab, and F. Varray, "On statistical methods beyond correlation in coherence analysis based on ultrasound," in *Proceedings of the 33rd European Signal Processing Conference (EUSIPCO)*, Sicily, Italy, 2025, to appear.

[22] A. Derode and M. Fink, "Spatial coherence of ultrasonic speckle in composites," *IEEE Transactions on ultrasonics, ferroelectrics, and frequency control*, vol. 40, no. 6, pp. 666-675, 1993.

[23] E. Turquin, L. Petrusca, M. Viallon, H. Liebgott, and F. Varray, "Full 3d anisotropic estimation of tissue in ultrasound imaging," in *2019 IEEE International Ultrasonics Symposium (IUS)*. IEEE, 2019, pp. 56-59.

[24] E. Zivot and J. Wang, "Vector autoregressive models for multivariate time series," *Modeling financial time series with S-PLUS®*, pp. 385-429, 2006.

[25] C. A. Sims, "Macroeconomics and reality," *Econometrica: journal of the Econometric Society*, pp. 1-48, 1980.

[26] E. Ollila and A. Breloy, "Regularized tapered sample covariance matrix," *IEEE Transactions on Signal Processing*, vol. 70, pp. 2306-2320, 2022.