Processing Beamformed IQ Data¶

Introduction¶

After loading your beamformed IQ data, the next step in most fUSI workflows is processing it to extract hemodynamic signals. This guide explains how to use ConfUSIus to transform beamformed IQ data into measures like B-mode, power Doppler, and axial velocity.

Why Process IQ Data?¶

Beamformed IQ data contains signals from all moving and stationary structures in the imaging field:

• Tissue signals: Strong echoes from slow-moving tissue (~40-60 dB stronger than blood).
• Blood signals: Weak echoes from fast-moving red blood cells.

Beamformed IQ signals can be processed to extract different types of information depending on your application:

• Structural imaging (B-mode): Use amplitude information to visualize tissue anatomy.
• Blood volume (power Doppler): Separate blood from tissue to quantify cerebral blood volume changes.
• Blood flow velocimetry (color Doppler): Estimate blood flow speed and direction.
• Tissue motion: Track tissue displacement (e.g., cardiac pulsatility, respiration).

This guide covers the IQ processing functions in ConfUSIus, which focus primarily on extracting blood volume information for fUSI applications. The typical workflow involves:

1. Filtering to isolate signals of interest (e.g., blood).
2. Computing derived measures from the filtered IQ signals (e.g., power Doppler, velocity).

For fUSI specifically, this means using clutter filtering to separate blood from tissue, then computing measures like power Doppler or axial velocity. The result is a series of volumes showing changes in cerebral blood volume or flow.

Understanding Your IQ Data¶

Before processing, it's important to understand what beamformed IQ data represents and how it's structured.

Quick Primer on Beamformed IQ Signals¶

Most fUSI users work with beamformed IQ signals rather than raw radio-frequency (RF) signals. Indeed, ultrasound systems emit pulses at high frequencies (typically 1–30 MHz). The received echoes oscillate at these same high frequencies. However, for imaging we only care about:

• Amplitude: How strong is the reflection (related to tissue properties)?
• Phase: How much has the phase shifted (related to motion via Doppler effect)?

IQ demodulation extracts this amplitude and phase information while discarding the high-frequency carrier, allowing for a significant reduction in sampling rate and thus drastically reducing data volume. A beamformed IQ sample can be represented as

where \(A\) is amplitude, \(\phi\) is phase, \((x, y, z)\) are spatial coordinates, and \(k\) is the slow-time index.

For single plane-wave imaging, the sampling rate of the beamformed IQ signals equals the pulse repetition frequency (PRF), typically kHz rates. With compound imaging (common in fUSI), multiple plane waves are emitted at different angles and their echoes are combined to form a single compound volume¹. In this case, the effective sampling rate, the compound frequency, is the PRF divided by the number of angles, typically hundreds of Hz.

Expected Data Structure¶

Beamformed IQ data in ConfUSIus should be a complex-valued Xarray DataArray with dimensions (time, z, y, x), where time is the slow-time dimension and (z, y, x) are the spatial dimensions. ConfUSIus typically adopts the following convention for the spatial dimensions:

• z: Elevation dimension (out-of-plane dimension for linear arrays).
• y: Axial dimension (depth dimension).
• x: Lateral dimension (dimension along the transducer array).

Additionally, ConfUSIus expects certain metadata attributes to be present:

• compound_sampling_frequency: Effective sampling frequency of the IQ data (Hz).
• transmit_frequency: Frequency of the transmitted ultrasound pulse (Hz).
• sound_velocity: Speed of sound used for beamforming (m/s).

import xarray as xr

iq = cf.load("sub-01_task-rest_iq.zarr")

print(iq)

<xarray.DataArray 'iq' (time: 50000, z: 1, y: 118, x: 52)> Size: 24GB
dask.array<...>
Coordinates:
  * time     (time) float64 400kB 0.0 0.002 0.004 ... 99.994 99.996 99.998
  * z        (z) float64 8B 0.0
  * y        (y) float64 944B 4.656 4.705 4.753 ... 10.231 10.279 10.328
  * x        (x) float64 416B -2.671 -2.57 -2.469 ... 2.268 2.369 2.469
Attributes:
    transmit_frequency:           15625000.0
    compound_sampling_frequency:  500.0
    pulse_repetition_frequency:   15000.0
    sound_velocity:               1540.0

In this example, the PRF is 15 kHz, but the data uses compound imaging with a compound_sampling_frequency of 500 Hz. This means each IQ sample (compound volume) results from combining 30 individual pulses, giving an effective temporal sampling period of 2 ms rather than 67 μs.

How ConfUSIus Processes IQ Data¶

ConfUSIus offers a flexible framework for processing IQ data using sliding windows. At its core is the process_iq_blocks function, which allows you to run any user-defined function on sliding windows of IQ signals along the temporal dimension. This base function is then wrapped by specialized processing functions for common fUSI applications.

Generic IQ Processing¶

The process_iq_blocks function is the foundation of IQ processing in ConfUSIus. It provides a generic interface for applying any user-defined function to "blocks" of IQ data using sliding temporal windows. You will generally not call this function directly, but it is worth knowing it exists if you ever need to implement a custom processing workflow (for example, when implementing new clutter filters).

---
config:
  layout: elk
---

flowchart LR
    IQ -->|Window 1| F1{{"User Function"}}
    IQ --> ellipsis{{"..."}}
    IQ -->|Window <em>n</em>| F2{{"User Function"}}

    F1 -->|"Output Volume(s)"| output
    F2 -->|"Output Volume(s)"| output

    IQ@{ shape: processes, label: "Input IQ Volumes" }
    output@{ shape: processes, label: "Stacked Output Volumes" }
    ellipsis@{ shape: text}

Convenience Wrappers¶

The specialized functions process_iq_to_power_doppler and process_iq_to_axial_velocity wrap process_iq_blocks to implement a nested sliding window approach. process_iq_to_bmode is similar but uses a single window level since no clutter filtering is involved:

1. Outer windows (clutter filtering): A first sliding window is used for clutter filtering. Within each window, the clutter filter is applied to separate blood from tissue signals.

2. Inner windows (derived measure computation): Inside each clutter-filtered window, a second sliding window runs to compute the desired measure (power Doppler or axial velocity). Multiple output volumes can be extracted from each clutter window.

---
config:
  layout: elk
---

flowchart LR
    IQ -->|"Outer Window 1"| CF1{{"Clutter Filter"}}
    IQ --> ellipsis{{"..."}}
    IQ -->|"Outer Window <em>n</em>"| CFn{{"Clutter Filter"}}

    CF1 -->|"Inner Window 1"| C1{{"Compute"}}
    CF1 --> e1{{"..."}}
    CF1 -->|"Inner Window <em>m</em>"| C2{{"Compute"}}

    CFn -->|"Inner Window 1"| C3{{"Compute"}}
    CFn --> e2{{"..."}}
    CFn -->|"Inner Window <em>m</em>"| C4{{"Compute"}}

    C1 -->|"Output Volume"| output
    C2 -->|"Output Volume"| output
    C3 -->|"Output Volume"| output
    C4 -->|"Output Volume"| output

    IQ@{ shape: processes, label: "Input IQ Volumes" }
    output@{ shape: processes, label: "Stacked Output Volumes" }
    ellipsis@{ shape: text}
    e1@{ shape: text}
    e2@{ shape: text}

This nested approach allows independent control over:

• Clutter filtering temporal resolution: Controlled by clutter_window_width and clutter_window_stride.
• Measure computation temporal resolution: Controlled by doppler_window_width / velocity_window_width and their corresponding strides.

Memory and Performance Optimization¶

Beamformed IQ datasets are typically very large (tens to hundreds of GB). ConfUSIus relies on the Dask backend for Xarray to enable out-of-core processing, allowing you to work with datasets larger than available RAM. Dask uses lazy evaluation, meaning data stays on disk until you explicitly compute results. This allows you to build complex processing pipelines that only execute when needed.

For optimal performance, always start a Dask distributed Client before running computations. Even for local processing (single machine), the distributed scheduler provides significantly better performance than the default scheduler, along with task monitoring and a web-based dashboard to track progress:

from dask.distributed import Client
import xarray as xr
import confusius

# Start a local Dask cluster.
client = Client()

# The client prints a dashboard link for monitoring progress.
# Example: http://127.0.0.1:8787/status

# Load IQ data (metadata only, no data in memory yet).
iq = cf.load("large_iq_data.zarr")

# Build computation graph (still no execution).
pwd = iq.fusi.iq.process_to_power_doppler(low_cutoff=40)

# Execute and save results (data processed in chunks, written to disk).
pwd.to_zarr("power_doppler.zarr")

# Limit workers and memory.
client = Client(
    n_workers=4,              # Number of worker processes.
    threads_per_worker=2,     # Threads per worker.
    memory_limit="8GB"        # Memory limit per worker.
)

import dask
import confusius

pwd = iq.fusi.iq.process_to_power_doppler(...)
velocity = iq.fusi.iq.process_to_axial_velocity(...)
bmode = iq.fusi.iq.process_to_bmode(...)

# Compute all three with a single pass over the IQ data.
pwd, velocity, bmode = dask.compute(pwd, velocity, bmode)

Clutter Filtering¶

Clutter filtering is an essential step in fUSI processing to separate weak blood signals from strong tissue signals. Indeed, blood signals from red blood cells are typically 40-60 dB (100-1000×) weaker than tissue signals, making effective clutter rejection essential for detecting hemodynamic changes.

How Clutter Filtering Works¶

Clutter filters exploit the different spatiotemporal characteristics of tissue and blood signals:

• Tissue signals: Spatially coherent across large regions, with low temporal frequencies from slow tissue motion.
• Blood signals: Spatially localized to vascular structures, with higher temporal frequencies from faster blood flow.

The challenge is that tissue motion (especially in awake subjects) can create Doppler frequencies that overlap with blood flow signals, making simple frequency-based separation insufficient.

Available Filter Methods¶

ConfUSIus provides several clutter filtering methods, each with different trade-offs:

Filter Parameters¶

The interpretation of cutoff parameters depends on the filter method:

Computing the Cumulative Energy Threshold¶

Setting high_cutoff for svd_cumulative_energy is non-trivial because cumulative energy values depend on the amount of tissue motion and change from window to window. compute_svd_cumulative_energy_threshold automates this: it iterates over sliding temporal windows, computes the cumulative eigenvalue spectrum for each window (accumulating from low-energy to high-energy components), and returns the minimum cumulative energy at a given singular-value index across all windows. Using the minimum ensures that the chosen cutoff does not over-filter any window.

import confusius as cf

iq = cf.load("recording.zarr")["iq"]
brain_mask = cf.load("brain_mask.zarr")["mask"]

# Compute the minimum cumulative energy at singular value index 40 across all windows.
threshold = cf.iq.compute_svd_cumulative_energy_threshold(
    iq,
    singular_value_index=40,
    clutter_mask=brain_mask,
    window_width=200,
).compute()

# Apply the cumulative energy filter using the computed threshold.
pwd = iq.fusi.iq.process_to_power_doppler(
    filter_method="svd_cumulative_energy",
    clutter_mask=brain_mask,
    high_cutoff=threshold,
    clutter_window_width=200,
)

Clutter Masks for SVD Filters¶

For SVD-based filters, you can provide a spatial mask identifying tissue regions. This improves filter performance by computing tissue subspaces only from masked voxels⁴:

import xarray as xr

# Load a brain mask (e.g., from an anatomical atlas or segmentation).
brain_mask = cf.load("brain_mask.zarr")

# Use the mask for clutter filtering.
pwd = iq.fusi.iq.process_to_power_doppler(
    filter_method="svd_indices",
    clutter_mask=brain_mask,
    low_cutoff=50,
)

Computing Derived Measures¶

ConfUSIus can compute several derived measures from beamformed IQ data. In fUSI, the primary measures are B-mode (tissue anatomy), power Doppler (cerebral blood volume changes), and axial velocity (blood flow speed and direction).

Power Doppler¶

Power Doppler allows estimation of cerebral blood volume by integrating the Doppler power spectrum over time. It is the most common fUSI readout and is ideally suited for mapping functional brain activity.

Computing Power Doppler¶

To compute power Doppler signals from beamformed IQ data, you may use either the confusius.iq.process_iq_to_power_doppler function or the corresponding Xarray accessor method.

from dask.distributed import Client
import xarray as xr
import confusius as cf

# Start a local Dask cluster for optimal performance.
client = Client()

# Load IQ data.
iq = cf.load("sub-01_task-rest_iq.zarr")

# Process to power Doppler using SVD filtering.
pwd = cf.iq.process_iq_to_power_doppler(
    iq,
    clutter_window_width=200,       # Use 200 volumes per clutter filter window.
    doppler_window_width=100,       # Integrate power in nested windows of 100 volumes.
    filter_method="svd_indices",    # Use the static SVD filter.
    low_cutoff=50,                  # Remove 50 strongest singular vectors.
)

from dask.distributed import Client
import xarray as xr
import confusius

# Start a local Dask cluster for optimal performance.
client = Client()

# Load IQ data.
iq = cf.load("sub-01_task-rest_iq.zarr")

# Process to power Doppler using SVD filtering.
pwd = iq.fusi.iq.process_to_power_doppler(
    clutter_window_width=200,       # Use 200 volumes per clutter filter window.
    doppler_window_width=100,       # Integrate power over 100 volumes.
    filter_method="svd_indices",    # Use static SVD filter.
    low_cutoff=50,                  # Remove 50 strongest singular vectors.
)

The resulting power Doppler DataArray will be backed by a Dask array (even if the input wasn't Dask-backed) and rely on lazy evaluation.

print(pwd)

<xarray.DataArray 'power_doppler' (time: 860, z: 1, y: 128, x: 86)> Size: 76MB
dask.array<...>
Coordinates:
  * time     (time) float64 7kB 0.299 0.899 1.499 2.099 ... 514.5 515.1 515.7
  * z        (z) float64 8B 0.0
  * y        (y) float64 1kB 5.664 5.713 5.762 5.811 ... 11.72 11.77 11.82 11.87
  * x        (x) float64 688B -3.583 -3.492 -3.402 -3.311 ... 3.946 4.037 4.127
Attributes: (11/16)
    transmit_frequency:           15625000.0
    compound_sampling_frequency:  500.0
    pulse_repetition_frequency:   15000.0
    sound_velocity:               1540.0
    units:                        'a.u.'
    long_name:                    'Power Doppler intensity'
    clutter_filter_method:        svd_indices
    clutter_window_width:         300
    clutter_window_stride:        300
    doppler_window_width:         300
    doppler_window_stride:        300
    clutter_low_cutoff:           40

To actually compute the power Doppler values and load them into memory, you must call .compute():

pwd = pwd.compute()
print(pwd)

<xarray.DataArray 'power_doppler' (time: 860, z: 1, y: 128, x: 86)> Size: 76MB
array([[[[74.26474299, 76.99585806, 84.46316397, ..., 79.11284661,
          77.6835962 , 69.19274666],
         [70.16886299, 83.84617566, 76.62624364, ..., 74.06344271,
          74.82030741, 74.48822094],
         ...
...

Axial Velocity¶

Axial velocity estimation quantifies the speed and direction of blood flow along the ultrasound beam axis using autocorrelation-based methods.

Computing Axial Velocity¶

To compute axial velocity from beamformed IQ data, you may use either the confusius.iq.process_iq_to_axial_velocity function or the corresponding Xarray accessor method.

from dask.distributed import Client
import xarray as xr
import confusius as cf

# Start a local Dask cluster for optimal performance.
client = Client()

# Load IQ data.
iq = cf.load("sub-01_task-rest_iq.zarr")

# Process to axial velocity using SVD filtering.
velocity = cf.iq.process_iq_to_axial_velocity(
    iq,
    clutter_window_width=50,
    velocity_window_width=50,
    filter_method="svd_indices",
    low_cutoff=50,
    lag=1,                          # Autocorrelation lag (volumes).
    estimation_method="average_angle",  # Velocity estimation method.
)

from dask.distributed import Client
import xarray as xr
import confusius

# Start a local Dask cluster for optimal performance.
client = Client()

# Load IQ data.
iq = cf.load("sub-01_task-rest_iq.zarr")

# Process to axial velocity using SVD filtering.
velocity = iq.fusi.iq.process_to_axial_velocity(
    clutter_window_width=50,
    velocity_window_width=50,
    filter_method="svd_indices",
    low_cutoff=50,
    lag=1,                          # Autocorrelation lag (volumes).
    estimation_method="average_angle",  # Velocity estimation method.
)

The resulting axial velocity DataArray will be backed by a Dask array (even if the input wasn't Dask-backed) and rely on lazy evaluation.

print(velocity)

<xarray.DataArray (time: 1000, z: 1, y: 118, x: 52)> Size: 24MB
dask.array<...>
Coordinates:
  * time     (time) float64 8kB 0.024 0.074 0.124 ... 99.874 99.924 99.974
  * z        (z) float64 8B 0.0
  * y        (y) float64 944B 4.656 4.705 4.753 ... 10.231 10.279 10.328
  * x        (x) float64 416B -2.671 -2.57 -2.469 ... 2.268 2.369 2.469
Attributes:
    units:    m/s

Interpreting Velocity¶

Velocity values are returned in meters per second with sign indicating direction:

• Positive values: Flow toward the transducer.
• Negative values: Flow away from the transducer.
• Magnitude: Speed of flow.

B-mode¶

B-mode imaging visualizes tissue anatomy by averaging the IQ signal magnitude over time. Unlike power Doppler, no clutter filtering is applied and the mean magnitude (not squared magnitude) is computed within each temporal window.

Computing B-mode¶

To compute B-mode volumes from beamformed IQ data, you may use either the confusius.iq.process_iq_to_bmode function or the corresponding Xarray accessor method.

from dask.distributed import Client
import xarray as xr
import confusius as cf

# Start a local Dask cluster for optimal performance.
client = Client()

# Load IQ data.
iq = cf.load("sub-01_task-rest_iq.zarr")

# Process to B-mode using a sliding window of 200 volumes.
bmode = cf.iq.process_iq_to_bmode(iq, bmode_window_width=200)

from dask.distributed import Client
import xarray as xr
import confusius

# Start a local Dask cluster for optimal performance.
client = Client()

# Load IQ data.
iq = cf.load("sub-01_task-rest_iq.zarr")

# Process to B-mode using a sliding window of 200 volumes.
bmode = iq.fusi.iq.process_to_bmode(bmode_window_width=200)

Next Steps¶

Now that you've processed your beamformed IQ data into power Doppler or velocity volumes, you're ready to:

1. Visualization: explore your data interactively with napari and generate static figures with Matplotlib.
2. Quality Control: assess data quality and identify artifacts (coming soon)
3. Registration: correct for motion and align to anatomical templates (coming soon)
4. Signal Processing: extract regional signals and apply denoising (coming soon)