Reference

CFunction{E, D, A, T, IP, IE} <: IsotropicEstimate{E, D}

Spatial C function estimate derived from spectral estimates.

The C function is the reduced covariance measure of a process evaluated on a punctured ball, of a function of distance. So at radius r it is C(r) = C̆({x : 0 < ||x|| ≤ r}), where C̆ is the reduced covariance measure.

Type Parameters

• E: Estimate trait (e.g., MarginalTrait, PartialTrait)
• D: Spatial dimension
• A: Type of radii array
• T: Type of C function values
• IP: Type of process information
• IE: Type of estimation information

Fields

• radii::A: Distances at which the C function is evaluated
• value::T: C function values
• processinformation::IP: Information about the processes
• estimationinformation::IE: Information about the estimation procedure

Mathematical Background

For a stationary process with power spectral density f(k), the C function is: C(r) = ∫ f(k) W(r,k) dk where W(r,k) is a spatial weighting function depending on the dimension.

Examples

# Compute C function from spectral estimate
cf = c_function(spectrum, radii=0.1:0.1:2.0)

# Direct computation from data
cf = c_function(data, region, radii=0.1:0.1:2.0, nk=(32,32), kmax=(0.5,0.5), tapers=tapers)

CenteredLFunction{E, D, A, T, IP, IE} <: IsotropicEstimate{E, D}

Centered L function estimate for spatial pattern analysis.

The centered L function is defined as $CL(r) = L(r) - r$, which removes the theoretical expectation for a Poisson process. This makes it easier to detect deviations from complete spatial randomness:

• $CL(r) = 0$ indicates Poisson behavior (complete spatial randomness)
• $CL(r) > 0$ indicates clustering at scale r
• $CL(r) < 0$ indicates regularity/inhibition at scale r

Type Parameters

• E: Estimate trait (e.g., MarginalTrait, PartialTrait)
• D: Spatial dimension
• A: Type of radii array
• T: Type of centered L function values
• IP: Type of process information
• IE: Type of estimation information

Mathematical Background

For a Poisson process, $L(r) = r$ theoretically. The centering transformation: $LCL(r) = L(r) - r$ removes this trend, making deviations more apparent and easier to interpret.

Examples

# Compute centered L function from L function
clf = centered_l_function(l_func)

# Direct computation from data
clf = centered_l_function(data, region, radii=0.1:0.1:2.0, nk=(32,32), kmax=(0.5,0.5), tapers=tapers)

Coherence{E, D, N, A, T, IP, IE} <: AbstractEstimate{E, D, N}

A coherence estimate structure containing wavenumber information and coherence values.

Type Parameters

• E: Estimate trait (e.g., MarginalTrait, PartialTrait)
• D: Spatial dimension
• N: Number of wavenumber dimensions
• A: Type of wavenumber argument
• T: Type of coherence values
• IP: Type of process information
• IE: Type of estimation information

Fields

• wavenumber: Wavenumber grid or values
• coherence: Coherence estimates
• processinformation: Information about the processes
• estimationinformation: Information about the estimation procedure

GaussKernel{T <: Real}

Gaussian smoothing kernel for rotational averaging.

Fields

• bw::T: Bandwidth parameter controlling the width of the Gaussian

Mathematical Form

K(x) = exp(-x²/(2σ²)) / σ where σ is the bandwidth.

(f::GaussKernel)(x)

Evaluate the Gaussian kernel at distance x.

KFunction{E, D, A, T, IP, IE} <: IsotropicEstimate{E, D}

Ripley's K function estimate for spatial point pattern analysis.

Ripley's K function $K(r)$ measures the expected number of points within distance r of a typical point, normalized by the intensity. It is fundamental in spatial statistics for detecting clustering or regularity in point patterns.

Type Parameters

• E: Estimate trait (e.g., MarginalTrait, PartialTrait)
• D: Spatial dimension
• A: Type of radii array
• T: Type of K function values
• IP: Type of process information
• IE: Type of estimation information

Fields

• radii::A: Distances at which the K function is evaluated
• value::T: K function values
• processinformation::IP: Information about the processes
• estimationinformation::IE: Information about the estimation procedure

Mathematical Background

For a stationary point process with intensity λ, Ripley's K function is: $Kᵢⱼ(r) = λ⁻¹ E$[number of additional i points within distance r of a typical j point]

The relationship to the C function is: $K(r) = C(r)/λ² + V_d r^d$ where V_d is the volume of a unit ball in d dimensions.

Examples

# Compute K function from spectral estimate
kf = k_function(spectrum, radii=0.1:0.1:2.0)

# Direct computation from data
kf = k_function(data, region, radii=0.1:0.1:2.0, nk=(32,32), kmax=(0.5,0.5), tapers=tapers)

LFunction{E, D, A, T, IP, IE} <: IsotropicEstimate{E, D}

L function estimate derived from Ripley's K function for spatial analysis.

The L function is a variance-stabilizing transformation of Ripley's K function: $L(r) = (K(r)/V_d)^(1/d)$ where V_d is the volume of a unit ball in d dimensions.

For a Poisson process, $L(r) = r$, making deviations easier to interpret.

Type Parameters

• E: Estimate trait (e.g., MarginalTrait, PartialTrait)
• D: Spatial dimension
• A: Type of radii array
• T: Type of L function values
• IP: Type of process information
• IE: Type of estimation information

Mathematical Background

The L function transformation provides:

• Better variance stabilization than K function
• Linear behavior $L(r) = r$ for Poisson processes
• Easier interpretation of clustering/regularity patterns

Examples

# Compute L function from K function
lf = l_function(k_func)

# Direct computation from data
lf = l_function(data, region, radii=0.1:0.1:2.0, nk=(32,32), kmax=(0.5,0.5), tapers=tapers)

MarginallyTransformedEstimate{E, D, N, S, A, T, IP, IE, F} <: AbstractEstimate{E, D, N}

A wrapper for estimates that have undergone element-wise transformations.

This structure represents an estimate where a mathematical function has been applied element-wise to the estimate values while preserving the argument structure and metadata. Common examples include taking the real part, magnitude, phase, or logarithm of estimates.

Type Parameters

• E: Estimate trait (e.g., MarginalTrait, PartialTrait)
• D: Spatial dimension
• N: Number of wavenumber dimensions
• S: Type of the original estimate before transformation
• A: Type of the argument (e.g., wavenumber grid)
• T: Type of the transformed estimate values
• IP: Type of process information
• IE: Type of estimation information
• F: Type of the transformation function

Fields

• argument: The argument (typically wavenumber) of the estimate
• estimate: The transformed estimate values
• processinformation: Information about the processes
• estimationinformation: Information about the estimation procedure

Examples

# Take magnitude of a coherence estimate
mag_coh = abs(coherence_estimate)

# Extract real part of cross-spectrum
real_spec = real(cross_spectrum)

# Compute log-magnitude spectrum
log_spec = log(abs(spectrum))

NoRotational

Indicates that no rotational averaging should be performed.

PartialResampler(
    data::NTuple{P, PointSet},
    region;
    tapers,
    nk,
    kmax,
    mean_method = DefaultMean(),
    shift_method,
    nk_marginal_compute,
    kmax_marginal_compute

) where {P}

Creates a PartialResampler object that can be used to generate partial spectra.

Trait system for handling different process structures in estimates.

• SingleProcessTrait: One process → D-dimensional array of numbers
• MultipleVectorTrait: Multiple processes as vector → (D+2)-dimensional array
• MultipleTupleTrait: Multiple processes as tuple → D-dimensional array of SMatrix

RectKernel{T <: Real}

Rectangular (uniform) smoothing kernel for rotational averaging.

Fields

• bw::T: Bandwidth parameter controlling the width of the rectangular window

Mathematical Form

K(x) = 1/h if |x/h| < 1/2, 0 otherwise, where h is the bandwidth.

(f::RectKernel)(x)

Evaluate the rectangular kernel at distance x.

RotationalEstimate{E, D, S, A, T, IP, IE} <: IsotropicEstimate{E, D}

An estimate that has undergone rotational averaging to produce isotropic results.

Rotational averaging transforms anisotropic estimates (which vary with direction) into isotropic estimates that depend only on radial distance from the origin. This is commonly used in spectral analysis to study scale-dependent properties independent of orientation.

Type Parameters

• E: Estimate trait (e.g., MarginalTrait, PartialTrait)
• D: Spatial dimension
• S: Type of the original anisotropic estimate
• A: Type of radii array
• T: Type of the rotationally averaged estimate values
• IP: Type of process information
• IE: Type of estimation information

Fields

• radii::A: Radial distances at which the estimate is evaluated
• estimate::T: Rotationally averaged estimate values
• processinformation::IP: Information about the processes
• estimationinformation::IE: Information about the estimation procedure

Examples

# Create rotational estimate with default parameters
rot_spec = rotational_estimate(anisotropic_spectrum)

# Create with custom radii and kernel
radii = 0.1:0.1:1.0
kernel = GaussKernel(0.05)
rot_spec = rotational_estimate(spectrum, radii=radii, kernel=kernel)

SpatialData{T,R}

A struct to hold spatial data and its associated region.

All methods in SpatialMultitaper.jl use SpatialData objects to represent spatial datasets. Most will allow the user to pass other raw formats directly and then convert them internally to SpatialData objects.

Fields

• data::T: The spatial observations (PointSet, GeoTable, or collections thereof)
• region::R: The spatial region where data is observed

Usage

You should not use this constructor directly. You should specify data and region using the spatial_data function, which will automatically unpack, validate, and mask the input data to the specified region.

StandardShift(shift)

A shift method that just applies a shift to all points. Does not make any assumptions about the region they are recorded on.

collect(estimate::AbstractEstimate)

Produces a Tuple containing the arguments and the estimate.

Base.getindex(estimate::AbstractEstimate, p, q)
Base.getindex(estimate::AbstractEstimate, p, q, i...)

Get a subset of the estimate corresponding to processes p and q and potentially spatial indices i.

If p and q are integers, the result will be a single process estimate. Otherwise the result will be the same type of estimate but with the specified processes. If one of the indices is an integer and the other is not, the result will still have the original number of dimensions internally, in the sense that if results are stored as Array then the size would be (1, Q, ...) or (P, 1, ...) as appropriate. If they are Array{SMatrix} then the size of the array is the same, but each element is an SMatrix{1, Q} or SMatrix{P, 1} as appropriate.

interp_ip_2(h::AbstractArray{T,D}, g::AbstractArray{T,D}, grid::CartesianGrid) where {T,D}

The L₂ inner product of interpolated sequences g and h recorded on grid grid.

_anisotropic_c_weight(r, u, ::Val{D})

Compute spatial weighting function for D-dimensional c function.

These functions represent the Fourier transform of spherical/circular domains.

_apply_centering(radius, l_value)

Apply centering transformation: $CL(r) = L(r) - r$.

_c_to_k_transform!(radii, c_values, trait, mean_prod, ::Val{D})

Transform C function values to K function values for multiple vector traits.

Handles the case where we have multiple processes stored as arrays.

_c_to_k_transform!(radii, c_values, trait, mean_prod, ::Val{D})

Transform C function values to K function values for single process or tuple traits.

_choose_wavenumbers_1d(nk, kmax)

Returns the wavenumbers for our fft interfaces, with nk wavenumbers, and kmax as the kmax wavenumber.

_compute_coherence_estimate!(spectrum, transform_func, trait_type)

Internal helper to compute coherence estimates with proper error handling.

_compute_k_from_c(radius, c_value, mean_prod, ::Val{2})

Optimized computation for 2D case where unit ball volume is π.

_compute_k_from_c(radius, c_value, mean_prod, ::Val{D})

Core computation for transforming C to K function values.

Implements $K(r) = C(r)/λ² + V_d r^d$ where $V_d$ is the volume of a unit d-ball.

_compute_l_from_k(k_value, ::Val{2})

Optimized L function computation for 2D: $L(r) = sign(K) * sqrt(|K|/π)$.

_compute_l_from_k(k_value, ::Val{D})

Core L function computation: $L(r) = (K(r)/V_d)^(1/d)$.

Handles sign preservation for negative K values by using: $L(r) = sign(K) * (|K|/V_d)^(1/d)$

_compute_spectral_matrix(x::AbstractVector)

Compute the outer product spectral matrix for a vector.

_compute_spectral_matrix(x::AbstractMatrix)

Compute the averaged spectral matrix for a matrix.

_construct_estimate_subset(::Type{<:MarginallyTransformedEstimate{...}}, trait, args...)

Internal constructor for creating estimate subsets with different traits.

_construct_estimate_subset(
    ::Type{<:RotationalEstimate{E, D, S}},
    trait::Type{<:EstimateTrait},
    argument, estimate, processinfo, estimationinfo
) where {E, D, S}

Construct a subset estimate with a different trait from an existing RotationalEstimate.

This is used to create partial or marginal estimates from a rotational estimate.

_count_processes(indices)

Count the number of processes from the index structure.

_dft_to_spectral_matrix!(
    S_mat::Array{<:SMatrix{P, P}}, J_n::AbstractArray{<:SVector{P}, N},
    ::MultipleTupleTrait) where {P, N}

Fill spectral matrix for multiple process case.

Each array in Jn has size n1 × ... × n_D × M.

_dft_to_spectral_matrix(J_n::AbstractArray, ::SingleProcessTrait)
_dft_to_spectral_matrix(J_n::AbstractArray, ::MultipleVectorTrait)
_dft_to_spectral_matrix(J_n::NTuple{P, AbstractArray}, ::MultipleTupleTrait) where {P}

Compute the spectral matrix from DFTs for different process types.

Arguments

• J_n: The tapered DFTs, with shape depending on the process type:
  • For SingleProcessTrait: n₁ × ... × n_D × M array (single process)
  • For MultipleVectorTrait: P × M × n₁ × ... × n_D array (multiple vector processes)
  • For MultipleTupleTrait: Tuple of P arrays, each of size n₁ × ... × n_D × M (multiple processes as tuple)
• Trait type: Specifies the process structure.

Returns

• The spectral matrix or array of spectral matrices, with shape depending on the input.

Notes

• For single process: returns an array of power spectral densities.
• For multiple processes (tuple): returns an array of spectral matrices.
• For multiple vector processes: returns a D+2 array where each slice is a spectral matrix.

_isotropic_c_weight(r, k, spacing, ::Val{D})

Compute weighting function for isotropic C functions.

For isotropic estimates, the weighting accounts for the radial integration that has already been performed.

_k_to_l_transform!(k_values::AbstractArray, ::Val{D})

Transform K function values to L function values for any array structure.

_l_to_centered_l_transform!(radii, l_values, ::MultipleVectorTrait)

Transform L function values to centered L function for multiple vector processes.

_l_to_centered_l_transform!(radii, l_values, ::Union{MultipleTupleTrait, SingleProcessTrait})

Transform L function values to centered L function for single process or tuple traits.

_make_wavenumber_grid(nk, kmax)

Create a wavenumber grid with dimension-specific parameters.

_nk_in_box

Internal function to compute the number of wavenumbers at which the Fourier transform of a taper was evaluated in given box.

_partial_from_marginal_error(operation_name, estimate_type)

Generate a descriptive error for unsupported partial operations on marginal estimates.

_rotational_estimate(a::AbstractEstimate, radii, ::NoRotational)

Pass-through function when no rotational averaging is requested.

_rotational_estimate(a::AbstractEstimate, radii, kernel)

Internal function to compute rotational averaging with a given kernel.

_smoothed_rotational(x::NTuple{D}, y::AbstractArray{T, D}, radii, kernel) where {D, T}

Core rotational averaging computation.

Computes weighted averages over circular annuli using the specified kernel.

_validate_estimate_dimensions(processinformation, est, P, Q)

Validate that estimate dimensions are consistent with process configuration.

apply_marginal_transform(transform, x::AbstractArray{<:Number})

Apply transformation to arrays of numbers.

For arrays of scalar values, applies the transformation element-wise across the array.

apply_marginal_transform(transform::F, x::AbstractArray{<:SMatrix}) where {F}

Apply transformation to arrays of static matrices.

For arrays containing static matrices (common in multi-process spectral estimates), applies the transformation element-wise to each matrix element.

apply_marginal_transform(transform::F, est::AbstractEstimate{E, D, N}) where {F, E, D, N}

Apply a transformation function element-wise to an estimate.

This is the core function for creating transformed estimates. It applies the given transformation function to each element of the estimate while preserving all metadata and argument structure.

Arguments

• transform::F: A function to apply element-wise (e.g., abs, real, log)
• est::AbstractEstimate: The original estimate to transform

Returns

A MarginallyTransformedEstimate wrapping the transformed values.

Examples

# Apply absolute value transformation
abs_est = apply_marginal_transform(abs, complex_estimate)

# Apply logarithm with custom function
log_est = apply_marginal_transform(x -> log(x + 1e-10), spectrum)

box2sides(box::Box)

Converts a box to a Tuple of Tuples, whose j'th entry is min and max of thej`th side.

c_function(data, region; kwargs...)

Compute spatial C function directly from data and region.

c_function(data::SpatialData; radii, nk, kmax, wavenumber_radii, rotational_method, spectra_kwargs...)

Compute spatial C function from spatial data.

First computes the power spectral density, then transforms to the C function via inverse Fourier transform with appropriate spatial weighting.

Arguments

• data::SpatialData: Input spatial data
• radii: Distances at which to evaluate the C function
• wavenumber_radii: Radial wavenumbers for rotational averaging (default: from nk, kmax)
• rotational_method: Kernel for rotational averaging (default: from nk, kmax)
• spectra_kwargs...: Additional arguments passed to spectra

Returns

A CFunction object containing the spatial C function.

c_function(spectrum::Spectra; radii, wavenumber_radii, rotational_method)

Compute spatial C function from a spectral estimate.

Arguments

• spectrum::Spectra: Input power spectral density estimate
• radii: Distances for C function evaluation
• wavenumber_radii: Radial wavenumbers for rotational averaging (default: from spectrum)
• rotational_method: Smoothing kernel for rotational averaging (default: from spectrum)

Returns

A CFunction object with the C function values.

centered_l_function!(c::CFunction)

Compute centered L function from C function via L function transformation.

centered_l_function!(k::KFunction)

Compute centered L function from K function via L function transformation.

centered_l_function!(l::LFunction{E, D}) where {E, D}

Transform L function to centered L function.

centered_l_function(data, region; kwargs...)

Compute centered L function directly from data and region.

centered_l_function(est::AbstractEstimate; kwargs...)

Compute centered L function from some estimate.

centered_l_function(data::SpatialData; kwargs...)

Compute centered L function from spatial data via L function transformation.

centerpad(x::AbstractArray, i)

Pads x with zeros on all sides by i elements. For example, if x = 1:2 and i = 1, then the output is [0, 1, 2, 0].

check_interpolated_tapers_cross(raw_tapers, grid; tol = 1e-2)

Checks how much the interpolation changes the orthogonality of the tapers. If the maximum cross-correlation is greater than 1e-2, it will warn the user.

check_observations(data, region)

Validate that spatial data is compatible with the specified region.

Ensures all data has the same parameter dimension as the region.

check_tapers_for_data(data, tapers, bandwidth; wavenumber_res, tol, min_concentration)

This function checks if the tapers are suitable for the provided data. In particular, it checks that the tapers work for the type of grids provided in the data. This assumes that the base tapers have an L2 norm of 1 and are continuous. We use the term grid here to refer to the sampling mechanism of a process. So it also refer to continuously recorded data (which is the NoGrid type).

If you have a problem, you can use taper_checks to check the tapers directly.

Importantly, these checks are very expensive to compute. In some cases failures in the checks may be due to the quality of the approximation of the Fourier transform of the tapers. Sometimes this can be mitigated by increasing the wavenumber resolution of the tapers here. However, when using interpolated tapers, you may also need to increase the resolution of the grid used to compute the Fourier transform of those tapers.

Arguments

• data: The data to check the tapers for.
• tapers: The tapers to check.
• bandwidth: The bandwidth of the tapers.
• wavenumber_res: The wavenumber resolution to use for the tapers generated from the data grids.
• tol: The tolerance for the L2 norm of the tapers.
• min_concentration: The minimum concentration of the tapers before they are considered invalid and an error is thrown.

checkprocessinformation(processinformation, est)

Validate that the ProcessInformation is consistent with the estimate array dimensions.

Logic

1. Determines process counts from the structure of processindices1 and processindices2
2. Validates estimate array dimensions against expected structure
3. For SMatrix estimates, validates against matrix dimensions

Process Count Rules

• Vector/Tuple indices → multiple processes (count = length)
• Single indices → single process (count = 1)

choose_concentration_resolution(tapers_on_grids, concentration_region)

Chooses the number of wavenumbers to sample at to integrate on the concentration region. Assumes that all tapers in each taper family are the same in some sense (the Fourier transform is sampled in the same way and interpolated). Does not assume this across families. Will return the lowest resolution to mitigate some effects of interpolation that can matter for fine resolutions.

choose_wavenumber_oversample(n, nk; maxoversample=100)

Find the smallest integer oversample such that nk*oversample ≥ n.

This function is needed because if nk ≥ n, we can simply pad the data to size nk. However, if nk < n, we need to pad to a larger size to be able to recover the desired wavenumbers.

coherence(data, region; kwargs...)

Compute coherence estimate directly from data and region.

coherence(x::Number)

Coherence for a single process (always returns 1).

coherence(data::SpatialData; kwargs...)

Compute coherence estimate from spatial data.

coherence(x::AbstractMatrix)

Compute coherence from a spectral matrix by normalizing with diagonal elements.

The coherence γᵢⱼ between processes i and j is computed as: γᵢⱼ = Sᵢⱼ / √(Sᵢᵢ * Sⱼⱼ)

Arguments

• x::AbstractMatrix: Spectral matrix

Returns

Coherence matrix with values bounded between 0 and 1 for the magnitude.

coherence(spectrum::NormalOrRotationalSpectra{E}) where {E}

Compute coherence from a spectral estimate.

Arguments

• spectrum: A Spectra or RotationalSpectra estimate

Returns

A Coherence object containing the coherence estimates.

Throws

• ArgumentError: If the spectrum does not have equal input and output process sets

Examples

# Compute coherence from spectral estimate
coh = coherence(spec)

# Direct computation from data
coh = coherence(data, region, nk=(32, 32), kmax=(0.5, 0.5), tapers=tapers)

covariance_zero_atom(data, region)

Estimate the zero atom of the reduced covariance measure for spatial data.

For point data, computes the density as the number of points divided by the region's measure. For marked point data, computes the sum of squared mark values divided by the region's measure. For gridded data, returns zero.

Returns a scalar for single datasets or a diagonal matrix for multiple datasets.

Arguments

• data: Spatial data (PointSet, GeoTable, CartesianGrid, or collections thereof)
• region: Spatial region for which to compute the estimate

Returns

• Scalar estimate for single datasets
• Diagonal matrix for multiple datasets (Tuple or Vector of datasets)

default_rotational_kernel(est::AbstractEstimate)
default_rotational_kernel(nk, kmax)
default_rotational_kernel(wavenumber::NTuple{D, AbstractVector{<:Real}}) where {D}

Construct a default smoothing kernel for rotational averaging.

The default kernel is a rectangular kernel with bandwidth twice the maximum wavenumber step size across dimensions. This provides reasonable smoothing while preserving wavenumber resolution.

Returns

A RectKernel with bandwidth = 2 * max(wavenumber_steps)

default_rotational_radii(wavenumber::NTuple{D,AbstractVector{<:Real}}) where {D}
default_rotational_radii(s::AbstractEstimate)
default_rotational_radii(nk, kmax)

Construct default radii for rotational averaging based on wavenumber vectors.

The default radii are chosen to span from near zero to the minimum of the maximum wavenumbers across dimensions, with the number of points equal to the maximum wavenumber vector length.

Arguments

• wavenumber: Tuple of wavenumber vectors for each dimension
• s: An estimate from which to extract wavenumber information
• nk, kmax: Number of wavenumbers and maximum wavenumber per dimension

Returns

A range of radii suitable for rotational averaging.

Algorithm

• Maximum radius = min(maxwavenumberper_dimension)
• Number of radii = max(lengthperdimension)
• Radii are offset by half a step to avoid the origin singularity

downsample(x::AbstractArray{T,D}, spacing) where {D,T}

Down sample an AbstractArray. Specify spacing as an Int, for same spacing in all dims. Specify spacing as an NTuple{D,Int} for different spacing. Specify spacing=nothing to just return x.

fft_anydomain(x::Array, grid::Grid, nk, kmax; kwargs...)

Apply an fft to data stored in x given that x is recorded on a grid grid. The output is at wavenumbers determined by nk and kmax. The output is fftshifted so that the zero wavenumber is in the middle of the output array. x can have more dimensions than grid but the first dimensions (up to the ndims(grid)) must be the same as the grid. When x has more dimensions than grid, the fft is only applied over the dimensions corresponding to the grid.

fft_only(points::PointSet, region; nk, kmax)

Does a non-uniform FFT of points, but moves them so the origin of the bounding box of the region is at zero. Makes no adjustments for tapering etc.

get_grid(data)

Gets the grid associated with data. If no grid is associated, it returns NoGrid(), i.e. if data is PointSet or domain(data) isa PointSet.

getargument(f::CenteredLFunction)

Get the radii at which the centered L function is evaluated.

getargument(f::KFunction)

Get the radii at which Ripley's K function is evaluated.

getargument(f::LFunction)

Get the radii at which the L function is evaluated.

getargument(f::RotationalEstimate)

Get the radii at which the rotational estimate is evaluated.

getestimate(f::CenteredLFunction)

Get the centered L function values $CL(r) = L(r) - r$.

getestimate(f::KFunction)

Get the K function values.

getestimate(f::LFunction)

Get the L function values.

getestimate(f::RotationalEstimate)

Get the rotationally averaged estimate values.

getestimatename(::Type{<:RotationalEstimate{E, D, S}}) where {E, D, S}

Get the name of the estimate, reflecting its rotational and original traits.

The name indicates the order of application of rotational and partial traits.

getoriginaltype(::Type{<:MarginallyTransformedEstimate{E, D, N, S}}) where {E, D, N, S}

Extract the original estimate type before transformation.

getoriginaltype(::Type{<:RotationalEstimate{E, D, S}}) where {E, D, S}

Extract the original estimate type before rotational averaging.

gettransformname(est::MarginallyTransformedEstimate)

Get a string representation of the transformation function name.

gettransformtype(est::MarginallyTransformedEstimate)

Get the type of transformation function applied to the estimate.

grid2sides(g::CartesianGrid)

Converts a grid to a tuple of the location of the centers for each side.

interpolated_taper_family(raw_tapers, grid; wavenumber_res=size(grid))

Constructs a multitaper object from multiple taper discrete impulse responses on a grid.

Arguments

• raw_tapers: A Vector of AbstractArrays containing the taper impulse response.
• grid: The grid corresponding to the tapers.
• wavenumber_res: Optional named argument specifying the upsampling for computing the transfer function.

Wavenumber sampling

• The wavenumber grid on which the transfer function will be sampled has:

- length `wavenumber_res`.
- interval `unitless_spacing(grid)/wavenumber_res`.

k_function(spectrum::NormalOrRotationalSpectra; kwargs...)

Compute Ripley's K function from a spectral estimate.

k_function(data, region; kwargs...)

Compute Ripley's K function directly from data and region.

k_function(c::CFunction{E, D}) where {E, D}

Compute Ripley's K function from a C function estimate.

Transforms the C function C(r) to Ripley's K function using: $K(r) = C(r)/λ² + V_d r^d$ where $λ$ is the process intensity and $V_d$ is the volume of a unit d-ball.

Arguments

• c::CFunction: Input C function estimate

Returns

A KFunction object with the corresponding K function values.

k_function(data::SpatialData; kwargs...)

Compute Ripley's K function from spatial data.

First computes the C function, then transforms it to the K function using the relationship $K(r) = C(r)/λ² + V_d r^d$.

Arguments

• data::SpatialData: Input spatial data
• kwargs...: Additional arguments passed to C function computation

Returns

A KFunction object containing Ripley's K function estimates.

l_function!(spec::NormalOrRotationalSpectra; kwargs...)

Compute L function from spectral estimates.

l_function!(k::KFunction{E, D}) where {E, D}

Transform Ripley's K function to L function.

l_function(data, region; kwargs...)

Compute L function directly from data and region.

l_function(est::AbstractEstimate)

Compute L function from some estimate.

l_function(data::SpatialData; kwargs...)

Compute L function from spatial data.

magnitude_coherence(args...; kwargs...)

Compute magnitude coherence directly from data arguments.

magnitude_coherence(coh::Coherence)

Compute magnitude of existing coherence estimate.

magnitude_coherence(spectrum::Spectra)

Compute magnitude coherence |γᵢⱼ| from a spectral estimate.

magnitude_squared_coherence(spectrum::Spectra)

Compute magnitude-squared coherence |γᵢⱼ|² from a spectral estimate.

Also known as the coherency or coherence function.

make_tapers(region; bandwidth, threshold = 0.99, space_res = 100, wavenumber_res = 500, wavenumber_generate_res = 500, ntapers_max = 30)

A function to create tapers for a given region.

Arguments

• region: The region to create the tapers on.
• bandwidth: The bandwidth of the tapers
• threshold: The threshold for the number of tapers to use.
• space_res: The resolution of the grid to create the tapers on.
• wavenumber_res: The resolution of the wavenumber grid to compute the transfer function on.
• wavenumber_generate_res: The resolution of the wavenumber grid to generate the tapers on.
• ntapers_max: The maximum number of tapers to generate.

As a rough rule of thumb, setting the bandwidth to an integer multiple (say x) of a side length of the bounding square of a region will produce x^d tapers (note this is a very crude approximation, but you will find out how many good tapers you can make when you create them).

mask(pp::PointSet, region)
mask(geotable::GeoTable, region)

Masks a PointSet or GeoTable to only include points within the specified region. If the data is already fully contained within the region, a copy is still made and returned. If the data is a GeoTable on a grid, values outside the region are set to NaN.

newweight(h, wavenumber_downsampling)

Necessary because eigs gives normalised vector under the downsampling.

nufft1d1_anydomain!(output_storage, input_data, nk, kmax, iflag, eps; kwargs...)

In place version of nufft1d1_anydomain. Use nufft1d1_anydomain_precomp to precompute the input data and memory required.

nufft1d1_anydomain(interval, nk, kmax, xj, cj, iflag, eps; kwargs...)

Computes the approximate dft of a 1d function using the nufft algorithm, with points in any interval, at wavenumbers defined by kmax and nk. xj are the points coordinates and cj are the function values. Set iflag<0 for -i in exponent.

nufft2d1_anydomain!(output_storage, input_data, nk, kmax, iflag, eps; kwargs...)

In place version of nufft2d1_anydomain. Use nufft2d1_anydomain_precomp to precompute the input data and memory required.

nufft2d1_anydomain(box, nk, kmax, xj, yj, cj, iflag, eps; kwargs...)

Computes the approximate dft of a 2d function using the nufft algorithm, with points in any box, at wavenumbers defined by kmax and nk. box here should be a tuple of intervals, one for each dimension. Similarly for kmax and nk. xj, yj are the points coordinates and cj are the function values. Set iflag<0 for -i in exponent.

nufft3d1_anydomain!(output_storage, input_data, nk, kmax, iflag, eps; kwargs...)

In place version of nufft3d1_anydomain. Use nufft3d1_anydomain_precomp to precompute the input data and memory required.

nufft3d1_anydomain(cube, nk, kmax, xj, yj, zj, cj, iflag, eps; kwargs...)

Computes the approximate dft of a 3d function using the nufft algorithm, with points in any cube, at wavenumbers defined by kmax and nk. cube here should be a tuple of intervals, one for each dimension. Similarly for kmax and nk. xj, yj, zj are the points coordinates and cj are the function values. Set iflag<0 for -i in exponent.

nufft_anydomain(region, nk, kmax, points::PointSet, c, iflag, eps; kwargs...)

Computes the approximate dft of a function using the nufft algorithm, with points in any region, at wavenumbers defined by kmax and nk.

Currently only 1d, 2d and 3d are supported.

observations(sd::SpatialData, i)

Return the ith spatial observations from a SpatialData object.

observations(sd::SpatialData)

Return the spatial observations from a SpatialData object.

optimaltapers(region::Meshes.GeometryOrDomain, grid::CartesianGrid; wavenumber_region::Meshes.GeometryOrDomain, ntapers::Int, wavenumber_res, wavenumber_downsample=nothing, real_tapers=true, tol=0.0)

Function to compute the optimal tapers for a given region, grid and region in wavenumber.

Arguments:

• region: The spatial observation region. Should be of type Geometry.
• grid: A CartesianGrid containing the region on which we wish to have a taper.
• wavenumber_region: The region in wavenumber we want to concentrate the taper on. Typically a Ball centered at zero.
• ntapers: The number of desired tapers.
• wavenumber_res: The oversampling to be used in wavenumber.
• real_tapers: Optional argument to decided if real or complex tapers should be provided. Default is true which provides real tapers.
• tol: Optional argument passed to the eigs function of Arpack. You likely need to play with this.
• check_grid: Optional argument to decide if the grid should be checked for compatibility with the region. Default is true. Sometimes due to float error this can fail, so you may want to set this to false if you are sure the grid and region are compatible.

Note about errors:

If you have an error, it is likely a convergence problem. Try setting a larger value for tol.

padto(x::AbstractArray{T,D}, n::NTuple{D,Int}) where {T,D}
padto(x::AbstractArray{T,D}, n::Int) where {T,D}

Pad x with zeros to size n. Important: the resultant array is size n, not size(x) .+ n. Note that size(x) ≤ n must hold, and the new array is of size n not size(x).+n. If an integer is provided for n, it is interpreted as the same size in all dimensions.

partial_coherence(data, region; kwargs...)

Compute partial coherence estimate directly from data and region.

partial_coherence(x::Number)

Partial coherence for a single process (always returns 1).

partial_coherence(data::SpatialData; kwargs...)

Compute partial coherence estimate from spatial data.

partial_coherence(x::AbstractMatrix)

Compute partial coherence from a spectral matrix.

Partial coherence removes the linear effects of all other processes.

partial_coherence(spectrum::NormalOrRotationalSpectra)

Compute partial coherence from marginal spectral estimates.

Arguments

• spectrum: A marginal spectral estimate

Returns

A Coherence{PartialTrait} object containing partial coherence estimates.

Throws

• ArgumentError: If the spectrum does not have equal input and output process sets

partial_from_resampled(
    J_cross::NTuple{Q},
    mean_cross,
    J_marginal,
    mean_marginal,
    atoms,
    wavenumber,
    J_original,
    mean_original,
    f_inv_cross,
    f_inv_marginal,
    ntapers::Int) where {Q}
    J_cross::NTuple{Q, AbstractArray{T, N}},
    J_marginal::NTuple{P, AbstractArray{T, N}},
    atoms,
    wavenumber,
    J_original::NTuple{P, AbstractArray{T, N}},
    f_inv_cross::Matrix{AbstractArray{SMatrix{R, Q, T, S}, D}}},
    f_inv_marginal::NTuple{P, AbstractArray{SMatrix{Q, Q, T, L}, D}},
    λ,
    ntapers::Int) where {P, Q, R, L, S, T, N, D}

Takes the original DFTs, the inverses of the cross and marginal spectral estimates from the original data, and the DFTs from the resampled data, and returns the partial spectra.

The cross inverses will have been premultiplied so that they take the form f[idx, idx]^{-1}f[idx, jdx]

partial_k_function(data, region; kwargs...)

Compute partial Ripley's K function directly from data and region.

partial_k_function(data::SpatialData; kwargs...)

Compute partial Ripley's K function from spatial data.

Partial K functions remove the linear influence of other processes.

partial_l_function(data, region; kwargs...)

Compute partial L function directly from data and region.

partial_l_function(data::SpatialData; kwargs...)

Compute partial L function from spatial data.

partial_magnitude_coherence(spectrum::Spectra{MarginalTrait})

Compute partial magnitude coherence from marginal spectral estimates.

partial_shift_resample(rng::AbstractRNG, statistic, resampler; kwargs...)

Assumes that statistic is of the form statistic(partialspectra::PartialSpectra; kwargs...)

The resampler should be able to generate resamples of the partial spectra.

partial_spectra(x::AbstractMatrix, ::Nothing)

Compute uncorrected partial spectra from a general matrix.

Uses explicit indexing for maximum clarity and type stability.

partial_spectra(spectrum::RotationalSpectra{MarginalTrait})

Compute partial spectral estimates from rotational marginal spectral estimates.

For rotational estimates, the bias correction is handled differently and is not currently fully implemented.

partial_spectra(x::AbstractMatrix{T}, ntapers::Int) where {T}

Compute bias-corrected partial spectra from a general matrix.

Arguments

• x: Square spectral matrix
• ntapers: Number of tapers used in original estimation

partial_spectra(data, region; kwargs...)

Compute partial spectral estimates directly from data and region.

partial_spectra(x::Number, ntapers)

Partial spectra for a single process (identity operation).

For a single process, the partial spectrum is identical to the original spectrum since there are no other processes to partial out.

partial_spectra(data::SpatialData; kwargs...)

Compute partial spectral estimates from spatial data.

First computes marginal spectra, then converts to partial spectra.

partial_spectra(x::SMatrix, ::Nothing)

Compute uncorrected partial spectra from a static matrix.

This is the core mathematical transformation that converts a spectral matrix to partial spectral matrix through matrix inversion and element-wise operations.

Mathematical Formula

For spectral matrix S with inverse G = S⁻¹:

• Diagonal: Pᵢᵢ = 1/Gᵢᵢ
• Off-diagonal: Pᵢⱼ = -Gᵢⱼ/(GᵢᵢGⱼⱼ - |Gᵢⱼ|²)

partial_spectra(spectrum::Spectra{MarginalTrait})

Compute partial spectral estimates from marginal spectral estimates.

Partial spectra remove the linear influence of all other processes, providing a measure of the direct relationship between process pairs. The computation involves matrix inversion and bias correction based on the number of tapers.

Arguments

• spectrum::Spectra{MarginalTrait}: A marginal spectral estimate

Returns

A Spectra{PartialTrait} object containing the partial spectral estimates.

Throws

• ArgumentError: If the spectrum does not have equal input and output process sets

Mathematical Details

For a spectral matrix S, the partial spectrum P is computed as:

• Pᵢᵢ = 1/Gᵢᵢ (diagonal elements)
• Pᵢⱼ = -Gᵢⱼ/(GᵢᵢGⱼⱼ - |Gᵢⱼ|²) (off-diagonal elements)

where G = S⁻¹ is the inverse spectral matrix.

Examples

# Compute partial spectra from marginal estimates
partial_spec = partial_spectra(marginal_spec)

# Direct computation from data::SpatialData
partial_spec = partial_spectra(data; nk=nk, kmax=kmax, tapers=tapers)

partial_spectra(x::SMatrix{Q, Q, T, N}, ntapers::Int) where {Q, T, N}

Compute bias-corrected partial spectra from a static matrix.

Applies finite-sample bias correction based on the number of tapers used in the original spectral estimation.

Arguments

• x: Square spectral matrix
• ntapers: Number of tapers used in original estimation

Mathematical Details

The bias correction factor is (M)/(M - Q + δᵢⱼ) where:

• M = number of tapers
• Q = number of processes
• δᵢⱼ = 1 if i=j (diagonal), 2 if i≠j (off-diagonal)

partial_spectra(x::SMatrix{2, 2, T, 4}, ::Nothing) where {T}

Optimized partial spectra computation for 2×2 static matrices.

This specialized method provides optimized computation for the common case of two-process analysis with explicit element-wise calculations.

partial_spectra_uncorrected(spectrum::Spectra{MarginalTrait})

Compute partial spectra without bias correction from the number of tapers.

This function computes partial spectra using the raw inverse relationship without the finite-sample bias correction that accounts for the number of tapers. Use with caution as results may be biased for small numbers of tapers.

Arguments

• spectrum::Spectra{MarginalTrait}: A marginal spectral estimate

Returns

A Spectra{PartialTrait} object with uncorrected partial spectral estimates.

phase(spectrum::Union{Spectra, Coherence})

Compute phase of coherence or cross-spectral estimates.

Returns the phase angle in radians.

point2unitlesstype(p::Point)

Return the type of the unitless coordinates of a Point.

points2coords(points::PointSet)

Return a Tuple, with j'th entry a vector of thejth entry of eachpointinpoints`.

precompute_nufft1d1_anydomain_input!(interval, nk, kmax, xj, cj)

Precomputes the input data required for nufft1d1_anydomain!.

precompute_nufft1d1_anydomain_output(nk, xj, cj)

Precomputes the memory required for nufft1d1_anydomain!.

precompute_nufft2d1_anydomain_input!(box, nk, kmax, xj, yj, cj)

Precomputes the input data required for nufft2d1_anydomain!.

precompute_nufft2d1_anydomain_output(box, nk, kmax, xj, yj, cj)

Precomputes the memory required for nufft2d1_anydomain!.

precompute_nufft3d1_anydomain_input!(cube, nk, kmax, xj, yj, zj, cj)

Precomputes the input data required for nufft3d1_anydomain!.

precompute_nufft3d1_anydomain_output(cube, nk, kmax, xj, yj, zj, cj)

Precomputes the memory required for nufft3d1_anydomain!.

prediction_kernel_ft(spec)

Computes the Fourier transform of the prediction kernel given estimates of the spectral density function.

process_trait(estimate::AbstractEstimate)

Determine the trait for an estimate based on its process information.

process_trait(data::MultipleSpatialDataVec)

Determine the trait for input spatial data.

reprocess(h_large, grid, wavenumber_downsample)

Undoes the downsampling and padding.

rescale_points(x, nk, kmax, interval; maxoversample = 100)

Rescales points to be in the interval [-π, π] and returns the oversampling factor. Oversampling factor is the smallest integer so that the interval is rescaled to fit in [-π, π], and the corresponding kmax wavenumber is a multiple of kmax. Interval just needs to have minimum and maximum defined.

rotational_estimate(est::AnisotropicEstimate{E}; radii, kernel) where {E}

Compute a rotationally averaged estimate from an anisotropic estimate.

Rotational averaging integrates the estimate over circular annuli at different radii, producing an isotropic result that depends only on distance from the origin. This is useful for analyzing scale-dependent properties without directional bias.

Arguments

• est::AnisotropicEstimate: The anisotropic estimate to average
• radii: Radial distances for evaluation (default: default_rotational_radii(est))
• kernel: Smoothing kernel for averaging (default: default_rotational_kernel(est))

Returns

A RotationalEstimate containing the rotationally averaged values. The exception is if kernel is NoRotational, in which case the input estimate is returned. This is used to skip rotational averaging in some downstream cases.

Examples

# Basic rotational averaging
rot_est = rotational_estimate(spectrum)

# With custom parameters
radii = range(0.1, 2.0, length=50)
kernel = GaussKernel(0.1)  # Gaussian smoothing with bandwidth 0.1
rot_est = rotational_estimate(spectrum, radii=radii, kernel=kernel)

make_sin_tapers(lbw, region::Box, gridsize)

Construct the sin tapers for a given region and gridsize.

Arguments:

• ntapers: The number of tapers.
• region: The region to construct the tapers on. Must be a box.

Background:

The sin tapers on a unit interval are defined as:

\[h_m(x) = \sqrt{2} * \sin(πxm)\]

single_geotable(data::GeoTable, idx::Integer)

Extract the idx-th column of a GeoTable as a new single-column GeoTable.

single_prediction_kernel_ft(S_mat::AbstractMatrix, idx::Int, ::Nothing)

Computes the Fourier transform of the prediction kernel for a single variable given an estimate of the spectral density function.

single_prediction_kernel_ft(S_mat::SMatrix{P,P,T,L}, idx::Int, ::Nothing)

Computes the Fourier transform of the prediction kernel for a single variable given an estimate of the spectral density function.

spatial_data(data, region::R, [names])

Creates a SpatialData object from raw spatial data and a specified region. Automatically unpacks, validates, and masks the input data to the specified region.

spectra(data, region; nk, kmax, tapers, mean_method = DefaultMean())

Compute the multitaper spectral estimate from a tapered DFT.

Arguments

• data: The data to estimate the spectrum from
• region: The region to estimate the spectrum from
• nk: The number of wavenumbers in each dimension
• kmax: The maximum wavenumber in each dimension
• dk: The wavenumber spacing in each dimension
• tapers: A tuple of taper functions
• mean_method::MeanEstimationMethod: The method to estimate the mean (default: DefaultMean())

If one of nk and kmax is specified, dk will be set to a default based on the region. Otherwise, you only need to specify two of the three parameters nk, kmax, and dk. They can either be scalars (applied uniformly across all dimensions) or tuples specifying each dimension. You can mix the two styles, e.g. nk=100 and kmax=(0.5, 1.0) for 2D data. nk must be an Int or tuple of Ints and should be positive, kmax and dk must be Real or tuples of Reals and also positive.

Returns

A Spectra object with wavenumber and power fields:

• wavenumber: D-dimensional NTuple of wavenumber arrays for each dimension
• power: Power spectral density in one of the following forms:
  • Single process: n_1 × ... × n_D array
  • NTuple{P} data: n_1 × ... × n_D array of SMatrix{P, P}
  • Vector of P processes: P × P × n_1 × ... × n_D array

Notes

• Indexing into a Spectra object returns a subset with the same wavenumber
• Use KnownMean(x) to specify a known mean value

Examples

spec = spectra(data, region, kmax=(0.5, 0.5), tapers=tapers)

tapered_dft(data::SpatialData, tapers, nk, kmax, mean_method=DefaultMean())

Compute the tapered discrete Fourier transform of spatial data.

Returns different output formats depending on the input data structure:

Single SpatialData

For a single process, returns a single array of the transformed data.

MultipleSpatialDataVec

For P processes stored as a vector, returns an array of size P × M × n1 × ... × nD, where M is the number of tapers and n1, ..., nD are the wavenumber dimensions.

MultipleSpatialDataTuple

For P processes stored as a tuple, returns a tuple of length P, where each entry is an array of size n1 × ... × nD × M.

Arguments

• data::SpatialData: The spatial data to transform
• tapers: Collection of taper functions to apply
• nk: Number of wavenumbers in each dimension
• kmax: Maximum wavenumber in each dimension
• mean_method::MeanEstimationMethod = DefaultMean(): Method for mean estimation

Notes

The mean_method should have the same length as the number of processes in data.

tapers_on_grid(tapers::TaperFamily, grid)

Converts continuous tapers to a discrete taper family on the provided grid. If the grid is NoGrid, it returns the tapers unchanged. This is used for checking if the tapers are suitable for the combination of grids present in the data.

unitless_coords(p::Point)

Return the coordinates of a point without units.

unitless_maximum(grid::CartesianGrid)

Return the maximum coordinates of a CartesianGrid without units.

unitless_measure(region)

Return the measure of a region without units.

unitless_minimum(grid::CartesianGrid)

Return the minimum coordinates of a CartesianGrid without units.

unitless_origin(grid::CartesianGrid)

Return the origin of a CartesianGrid without units.

unitless_spacing(grid::CartesianGrid)

Return the spacing of a CartesianGrid without units.

unpack_observations(data)

Unpacks input data into a supported format for spatial analysis.

Behavior

• PointSet → returned as-is
• GeoTable with single column → returned as-is
• GeoTable with multiple columns → error for marked points, unpacked for grids
• Tuple of spatial data → recursively unpacked and flattened
• Vector of spatial data → recursively unpacked and concatenated

Examples

# Single process
unpack_observations(pointset) # → pointset

# Multiple grid fields
unpack_observations(grid_geotable) # → (field1_geotable, field2_geotable, ...)

# Mixed collection
unpack_observations((pointset, grid)) # → (pointset, field1, field2, ...)

Notes

Marked point processes with multiple marks are not supported - pass as separate GeoTables instead.

unwrap_fft_output(J::AbstractArray, kmaxrel::NTuple{D, Int}) where {D}

Copies the result of the fft to give something which goes up to some multiple of the nyquist wavenumber (whilst keeping the number of output wavenumbers).

This is needed because if we want wavenumbers higher than the nyquist wavenumber, we can use periodicity to recover wavenumbers up to kmax = kmaxrel * nyquist. This works because at this point we assume the fft input is centered at zero, so the output has this periodicity. The final output of the calling function is adjusted appropriately for shifts in the input.

unwrap_index(n::Int, a::Int)

The indices required to unwrap the fft output to go up to n with a given oversampling a.

upsample(x::AbstractArray, grid::CartesianGrid, wavenumber_downsample)

Upsamples an array using linear interpolation, checking the grid size against x. Note this is specific to downsample as used here, and shouldn't be used for general upsampling.

wavenumber_downsample_index(nk, oversample)

Returns the indices required to downsample the oversampled wavenumbers with a given oversampling. This assumes that the wavenumbers are fftshifted.

Arguments

• nk::Int: The number of wavenumbers for the desired output.
• oversample::Int: The oversampling factor used.

wavenumber_downsample_startindex(nk,oversample)

Returns the starting index for the oversampled wavenumbers to downsample.

LinearOperator

A linear operator. Assumes that inputspace and outputspace are defined. One should also define mul!(y, L::LinearOperator, x) and preallocate_output(L::LinearOperator, x) for efficiency.

compute_eigenfunctions(R::AbstractArray,K::AbstractArray,ntapers=1;imag_tol=1e-5,real_tapers=true, tol=0.0)

Computes the eigenfunctions of the concentration operator. Provide real_tapers=true for real valued tapers.

make_concentration_operator(R::AbstractArray,K::AbstractArray,real_tapers=Val{true}())

Constructs the concentration operator. Provide real_tapers=Val{true}() for real valued tapers. Choose real_tapers=Val{false}() for complex valued tapers.