Encoding workflows

Slepian/DPSS sequences

Discrete prolate spheroidal sequences concentrate energy in a frequency band while minimizing spectral leakage. Use them for resting-state data or any analysis focused on a specific frequency range.

spec_slep <- spec_time_slepian(tr = 2, bandwidth = 0.08, k = 4)
lat_slep <- encode(X, spec_slep, mask = mask_vol, materialize = "matrix")

# 4 basis vectors spanning the 0-0.08 Hz band
dim(basis(lat_slep))
#> [1] 20  4

# Reconstruction error
sqrt(mean((as.matrix(lat_slep) - X)^2))
#> [1] 0.8974354

The bandwidth parameter (in Hz) and tr (in seconds) together determine the normalized half-bandwidth W = bandwidth * tr. The number of well-concentrated tapers is approximately 2 * n_time * W - 1; setting k lower than this gives heavier compression.

DCT

The discrete cosine transform provides orthogonal frequency components. It is fast and parameter-free beyond choosing k — a good default when you have no strong prior about the signal bandwidth.

lat_dct <- encode(X, spec_time_dct(k = 8), mask = mask_vol, materialize = "matrix")
dim(basis(lat_dct))
#> [1] 20  8
sqrt(mean((as.matrix(lat_dct) - X)^2))
#> [1] 0.7783513

B-splines

B-splines produce smooth, localized basis functions with controllable flexibility. They are useful for modeling slow drifts or the hemodynamic response, where you want smooth temporal structure without frequency-domain assumptions.

lat_bs <- encode(X, spec_time_bspline(k = 6, degree = 3), mask = mask_vol,
                 materialize = "matrix")
dim(basis(lat_bs))
#> [1] 20  6

The degree parameter controls smoothness (3 = cubic, the default). Set orthonormalize = TRUE (the default) for orthonormal columns, which simplifies downstream linear algebra.

Comparing temporal families

All three temporal encodings follow the same pattern. The trade-off is between spectral precision (Slepian), simplicity (DCT), and localization (B-spline):

rmse <- function(lat) sqrt(mean((as.matrix(lat) - X)^2))

data.frame(
  family = c("Slepian (k=4)", "DCT (k=8)", "B-spline (k=6)"),
  k      = c(ncol(basis(lat_slep)), ncol(basis(lat_dct)), ncol(basis(lat_bs))),
  rmse   = round(c(rmse(lat_slep), rmse(lat_dct), rmse(lat_bs)), 4)
)
#>           family k   rmse
#> 1  Slepian (k=4) 4 0.8974
#> 2      DCT (k=8) 8 0.7784
#> 3 B-spline (k=6) 6 0.8239

Graph Slepian

Graph Slepian bases are eigenvectors of the graph Laplacian, spectrally concentrated on spatial clusters. They are the spatial analog of temporal DPSS sequences.

lat_sp <- encode(X, spec_space_slepian(k = 3, k_neighbors = 6),
                 mask = mask_vol)

The k parameter sets the number of components per cluster; k_neighbors controls the k-NN graph used for spatial adjacency.

Heat wavelet

Heat wavelets approximate diffusion on the voxel graph at multiple time scales, capturing both local and global spatial structure.

lat_heat <- encode(X,
  spec_space_heat(scales = c(1, 2, 4, 8), order = 30, threshold = 1e-6),
  mask = mask_vol
)

The scales parameter controls diffusion time scales (larger = smoother); order sets the Chebyshev polynomial approximation order; threshold truncates small coefficients for sparsity.

HRBF (hierarchical radial basis functions)

HRBF places smooth Gaussian atoms at data-driven centers. This family gives you explicit, interpretable spatial basis functions.

lat_hrbf <- encode(X,
  spec_space_hrbf(params = list(
    sigma0 = 2, levels = 1,
    radius_factor = 2.5, kernel_type = "gaussian"
  )),
  mask = mask_vol
)

The sigma0 parameter controls atom width; levels controls hierarchy depth (more levels = finer resolution); radius_factor sets atom spacing relative to sigma.

Wavelet active (CDF 5/3 lifting)

The wavelet active encoder applies CDF 5/3 lifting wavelets along each spatial axis, touching only in-mask voxels. It is exact (biorthogonal) and very fast — the best choice for speed-critical pipelines.

lat_wav <- encode(X,
  spec_space_wavelet_active(levels_space = 2, levels_time = 0, threshold = 0),
  mask = mask_vol
)

# Wavelet active returns an ImplicitLatent; reconstruct with predict()
recon_wav <- predict(lat_wav)
dim(recon_wav)

Setting threshold > 0 zeroes out small wavelet coefficients for lossy compression. With threshold = 0 the transform is exact.

PCA (cluster-local)

Parcel-local PCA computes per-cluster SVD to build a block-sparse spatial dictionary. Each parcel contributes k components. This is a data-driven approach — the loadings are subject-specific.

lat_pca <- encode(X, spec_space_pca(k = 3, center = TRUE), mask = mask_vol,
                  reduction = red)

# 4 parcels x 3 components = 12 total spatial atoms
dim(loadings(lat_pca))
#> [1] 64 12

The center parameter subtracts voxel means before PCA (stored in offset()); whiten = TRUE rescales scores to unit variance.

Shared parcel templates

When you have a parcellation (atlas) and want the same spatial dictionary across subjects, use parcel_basis_template(). The default basis is Laplacian eigenvectors — smooth functions from the voxel adjacency graph within each parcel. No training data is needed.

# Build template once from the shared reduction.
tmpl <- parcel_basis_template(red, basis_pca(k = 1), data = X, center = TRUE)
tmpl
#> ParcelBasisTemplate
#>   Basis family: spec_pca 
#>   Atoms: 4 
#>   Voxels: 64 
#>   Parcels: 4 
#>   Center at encode: TRUE

# Project each subject onto the shared dictionary
lvec_s1 <- encode(X, spec_space_parcel(tmpl), mask = mask_vol)

# Loadings are identical across subjects; scores and offset vary
dim(basis(lvec_s1))     # time x atoms
#> [1] 20  4
dim(loadings(lvec_s1))  # voxels x atoms
#> [1] 64  4

The example uses one PCA component per parcel so it is fast and fully reproducible. Geometric templates use the same encoding surface:

atlas <- neuroim2::ClusteredNeuroVol(mask_vol, parcel_map)
tmpl_slepian <- parcel_basis_template(atlas, basis_slepian(k = 5))
lvec <- encode(X, spec_space_parcel(tmpl_slepian), mask = mask_vol)

Any lift()-compatible basis spec works: basis_slepian() (default), basis_pca(), or basis_heat_wavelet().