API Reference

Notes

FIF_1D and FIF_ND no longer accept the old combined kernel_construction_method argument. Use kernel_construction_method_flux and kernel_construction_method_observable instead.

Naive FIF kernels remain available for comparison only and emit a warning because their outputs are not remotely accurate.

Simulation

Simulation module for multifractal fields

Contains methods for generating synthetic multifractal time series and fields.

scaleinvariance.fBm_1D_circulant(size, H, periodic=True)[source]

Generate fractional Brownian motion by spectral synthesis.

Creates fBm with power spectrum |k|^(-β) where β = 2H + 1. Uses random phases and inverse FFT.

Parameters:

  • size (int) – Number of points (must be power of 2)

  • H (float) – Hurst parameter. Typical range is (0, 1) for fBm, but negative values are also supported (H < 0 gives blue noise with β < 1).

  • periodic (bool, optional) – If True (default), returns full periodic simulation suitable for periodic boundary conditions. If False, doubles simulation size internally then returns only the first half to eliminate periodicity artifacts.

Returns:

Generated 1D fBm, zero mean, unit std

Return type:

numpy.ndarray

scaleinvariance.fBm_ND_circulant(size, H, periodic=True)[source]

Generate N-dimensional fractional Brownian motion by spectral synthesis.

Creates N-D fBm with rotationally symmetric power spectrum |k|^(-β) where β = 2H + D. Uses random phases and inverse FFT to generate N-D random field.

Parameters:

  • size (tuple of ints) – Size of N-D array. Each dimension must be power of 2.

  • H (float) – Hurst parameter. Typical range is (0, 1) for fBm, but negative values are also supported (H < 0 gives blue noise with β < D).

  • periodic (bool or tuple of bool, optional) – If True (default), returns full periodic simulation suitable for periodic boundary conditions. If False, doubles simulation size internally along each dimension then returns only the first octant to eliminate periodicity artifacts. Can also be a tuple of bool for per-axis control.

Returns:

Generated N-D fBm field, zero mean, unit std

Return type:

numpy.ndarray

scaleinvariance.fBm_1D(size, H, causal=True, outer_scale=None, outer_scale_width_factor=2.0, periodic=True, gaussian_noise=None, kernel_construction_method='LS2010')[source]

Generate 1D fractional Brownian motion by fractional integration of Gaussian noise.

This method generates fBm by convolving Gaussian white noise with a power-law kernel. Supports extended Hurst range H in (-0.5, 1.5).

Algorithm:

  • For H in (-0.5, 0.5): Convolve noise with kernel |x|^(H - 1/2)

  • For H = 0.5: Integrate noise (standard Brownian motion)

  • For H in (0.5, 1.5): Integrate noise, then convolve with kernel |x|^(H - 3/2)

Parameters:

  • size (int) – Length of simulation (must be even, power of 2 recommended)

  • H (float) – Hurst exponent in (-0.5, 1.5). Controls correlation structure. - H in (-0.5, 0.5): Fractional integration with anti-correlation to mild correlation - H = 0.5: Standard Brownian motion - H in (0.5, 1.5): Integration + fractional integration for smooth processes

  • causal (bool, optional) – Use causal kernels (future doesn’t affect past). Default True.

  • outer_scale (int, optional) – Large-scale cutoff. Defaults to size.

  • outer_scale_width_factor (float, optional) – Controls transition width for outer scale. Default is 2.0.

  • periodic (bool, optional) – If True (default), returns full periodic simulation suitable for periodic boundary conditions. If False, doubles simulation size internally then returns only the first half to eliminate periodicity artifacts.

  • gaussian_noise (torch.Tensor or numpy.ndarray, optional) – Pre-generated Gaussian noise for reproducibility. Must have same size as simulation.

  • kernel_construction_method (str, optional) – Method for constructing convolution kernels. Options: - ‘naive’: Simple power-law kernels without corrections - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections (default) - ‘spectral’: Exact spectral-response kernel for non-causal periodic runs

Returns:

1D array of simulated fBm values (zero mean)

Return type:

numpy.ndarray

Raises:

ValueError – If size is not even or H is outside valid range (-0.5, 1.5)

Examples

>>> # Standard fBm with H=0.7
>>> fbm = fBm_1D(2048, H=0.7)
>>>
>>> # Anti-correlated process with H=-0.3
>>> fbm_blue = fBm_1D(2048, H=-0.3)
>>>
>>> # Super-smooth process with H=1.2
>>> fbm_smooth = fBm_1D(2048, H=1.2)

Notes

  • Computational complexity is O(N log N) due to FFT-based convolutions

  • Uses LS2010 finite-size corrected kernels by default

  • For H in standard range (0, 1), fBm_1D_circulant may be more efficient

scaleinvariance.FIF_1D(size, alpha, C1, H, levy_noise=None, causal=True, outer_scale=None, outer_scale_width_factor=2.0, kernel_construction_method_flux='LS2010', kernel_construction_method_observable='LS2010', periodic=True)[source]

Generate a 1D Fractionally Integrated Flux (FIF) multifractal simulation.

FIF is a multiplicative cascade model that generates multifractal fields with specified Hurst exponent H, intermittency parameter C1, and Levy index alpha. The method follows Lovejoy & Schertzer 2010, including finite-size corrections.

Returns field normalized by mean.

Algorithm:

  1. Generate extremal Lévy noise with stability parameter alpha

  2. Convolve with corrected kernel (finite-size corrected |k|^(-1/alpha)) to create log-flux

  3. Scale by (C1)^(1/alpha) and exponentiate to get conserved flux

  4. For H ≠ 0: Convolve flux with kernel |k|^(-1+H) to get observable

  5. For H < 0: Apply differencing to handle negative Hurst exponents

Parameters:

  • size (int) – Length of simulation (must be even, power of 2 recommended)

  • alpha (float) – Lévy stability parameter in [0.5, 2] and != 1. Controls noise distribution. Values < 0.5 are currently not implemented for FIF simulation.

  • C1 (float) – Codimension of the mean, controls intermittency strength. Must be > 0 for multifractal behavior. Special case: C1 = 0 routes to fBm (no intermittency), but requires causal=False since fBm cannot be causal.

  • H (float) – Hurst exponent in (-1, 1). Controls correlation structure.

  • levy_noise (torch.Tensor, optional) – Pre-generated Lévy noise for reproducibility. Must have same size as simulation.

  • causal (bool, optional) – Use causal kernels (future doesn’t affect past). Default True. False gives symmetric (non-causal) kernels.

  • outer_scale (int, optional) – Large-scale cutoff. Defaults to size.

  • outer_scale_width_factor (float, optional) – Controls transition width for outer scale. Transition width = outer_scale * width_factor. Default is 2.0.

  • kernel_construction_method_flux (str, optional) –

    Method for constructing flux (cascade) kernel. Options: - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections (default) - ‘spectral’: spectral-response kernel for non-causal periodic runs,

    LS2010 fallback for causal runs

    • ’naive’: Simple power-law kernels without corrections

  • kernel_construction_method_observable (str, optional) –

    Method for constructing observable (H) kernel. Options: - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections (default) - ‘spectral’: spectral-response kernel for non-causal periodic runs,

    LS2010 fallback for causal runs

    • ’naive’: Simple power-law kernels without corrections

    • ’spectral_odd’: Odd (antisymmetric) spectral kernel — left half negated

  • periodic (bool, optional) – If True (default), returns full periodic simulation suitable for periodic boundary conditions. If False, doubles simulation size internally then returns only the first half to eliminate periodicity artifacts.

Returns:

1D array of simulated multifractal field values

Return type:

numpy.ndarray

Raises:

  • ValueError – If size is not even or H is outside valid range (-1, 1)

  • ValueError – If C1 < 0 or if C1 = 0 with causal=True (fBm cannot be causal)

  • ValueError – If provided levy_noise doesn’t match specified size

  • ValueError – If alpha < 0.5 or alpha == 1 (not implemented)

Examples

>>> # Basic multifractal with strong intermittency
>>> fif = FIF_1D(1024, alpha=1.8, C1=0.1, H=0.3)
>>>
>>> # Monofractal case (C1=0 routes to fBm, requires causal=False)
>>> fbm = FIF_1D(1024, alpha=1.8, C1=0, H=0.3, causal=False)

Notes

  • Computational complexity is O(N log N) due to FFT-based convolutions

  • Large C1 values (> 0.5) can produce extreme values requiring careful handling

  • alpha < 0.5 and alpha = 1 are not currently implemented

  • C1 = 0: Routes internally to fBm_1D_circulant() for monofractal case (requires causal=False since fBm cannot be causal)

scaleinvariance.FIF_ND(size, alpha, C1, H, levy_noise=None, outer_scale=None, outer_scale_width_factor=2.0, kernel_construction_method_flux='LS2010', kernel_construction_method_observable='LS2010', periodic=False, scale_metric=None, scale_metric_dim=None)[source]

Generate an N-D Fractionally Integrated Flux (FIF) multifractal simulation.

FIF is a multiplicative cascade model that generates multifractal fields with specified Hurst exponent H, intermittency parameter C1, and Levy index alpha. The method follows Lovejoy & Schertzer 2010, including finite-size corrections. Returns field normalized by mean.

Algorithm:

  1. Generate N-D extremal Lévy noise with stability parameter alpha

  2. Convolve with corrected kernel (finite-size corrected |r|^(-d/alpha)) to create log-flux

  3. Scale by (C1)^(1/alpha) and exponentiate to get conserved flux

  4. For H ≠ 0: Convolve flux with kernel |r|^(-d+H) to get observable

Parameters:

  • size (tuple of ints) – Size of simulation as tuple specifying dimensions (e.g., (height, width) for 2D, (depth, height, width) for 3D). All dimensions must be even numbers.

  • alpha (float) – Lévy stability parameter in [0.5, 2] and != 1. Controls noise distribution. Values < 0.5 are currently not implemented for FIF simulation.

  • C1 (float) – Codimension of the mean, controls intermittency strength. Must be > 0 for multifractal behavior. Special case: C1 = 0 routes to fBm_ND_circulant() (no intermittency).

  • H (float) – Hurst exponent in (0, 1). Controls correlation structure.

  • levy_noise (ndarray, optional) – Pre-generated N-D Lévy noise for reproducibility. Must have same shape as simulation.

  • outer_scale (int, optional) – Large-scale cutoff. Defaults to max(dimensions).

  • outer_scale_width_factor (float, optional) – Controls transition width for outer scale. Transition width = outer_scale * width_factor. Default is 2.0.

  • kernel_construction_method_flux (str, optional) – Method for constructing flux (cascade) kernel. Options: - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections (default) - ‘spectral’: exact isotropic spectral-response kernel for periodic N-D runs

  • kernel_construction_method_observable (str, optional) – Method for constructing observable (H) kernel. Options: - ‘LS2010’: Lovejoy & Schertzer 2010 finite-size corrections (default) - ‘spectral’: exact isotropic spectral-response kernel for periodic N-D runs

  • periodic (bool or tuple of bool, optional) – Controls periodicity behavior for each axis. - If bool: applies same periodicity to all axes (default False) - If tuple: specifies per-axis periodicity, must have length matching size tuple When an axis is non-periodic (False), the domain is doubled along that axis internally and one hyper-octant is returned to suppress periodic artifacts.

  • scale_metric (ndarray, optional) – Custom N-D array defining the scale metric (generalized distance) for GSI norms. If None (default), standard Euclidean distance is computed. Shape requirement: Must match simulation domain shape after accounting for periodicity (doubled along non-periodic axes). For example, if size=(256, 256) and periodic=False, scale_metric.shape must be (512, 512). Note: For LS2010 method, the metric should use dx=2 spacing to match the kernel construction grid.

  • scale_metric_dim (float, optional) – Dimension for scaling exponents in kernel calculations. This may differ from the spatial dimension when using GSI norms with non-integer dimensions. If None (default), uses the spatial dimension (inferred from size tuple length). Use case: In Generalized Scale Invariance (GSI), the scaling dimension may be non-integer and different from the embedding space dimension.

Returns:

N-D array of simulated multifractal field values

Return type:

numpy.ndarray

Raises:

  • ValueError – If size dimensions are not even or H is outside valid range (0, 1)

  • ValueError – If C1 < 0 (must be non-negative; C1 = 0 routes to fBm)

  • ValueError – If provided levy_noise doesn’t match specified size

  • ValueError – If alpha < 0.5 or alpha == 1 (not implemented)

Examples

>>> # 2D multifractal with strong intermittency
>>> fif = FIF_ND((512, 512), alpha=1.8, C1=0.1, H=0.3)
>>>
>>> # 3D multifractal
>>> fif = FIF_ND((128, 128, 128), alpha=1.5, C1=0.05, H=0.2)
>>>
>>> # 2D with per-axis periodicity (periodic in x, non-periodic in y)
>>> fif = FIF_ND((256, 256), alpha=1.6, C1=0.08, H=0.4, periodic=(True, False))
>>>
>>> # Monofractal case (C1=0 routes to fBm_ND_circulant)
>>> fbm = FIF_ND((512, 512), alpha=1.8, C1=0, H=0.3)
>>>
>>> # Using custom GSI norm (example with anisotropic scaling)
>>> import numpy as np
>>> size = (256, 256)
>>> sim_size = (512, 512)  # doubled for non-periodic
>>> # Create coordinate arrays with dx=2 spacing (to match LS2010)
>>> x = np.arange(-(sim_size[0]-1), sim_size[0], 2)
>>> y = np.arange(-(sim_size[1]-1), sim_size[1], 2)
>>> X, Y = np.meshgrid(x, y, indexing='ij')
>>> # Define anisotropic GSI metric with different scaling in x and y
>>> scale_metric = (X**2 + (Y/2)**2)**0.5  # y-direction scales faster
>>> # Use non-integer dimension for GSI
>>> fif = FIF_ND(size, alpha=1.7, C1=0.1, H=0.3, periodic=False,
...              scale_metric=scale_metric, scale_metric_dim=2.3)

Notes

  • Computational complexity is O(N log N) due to FFT-based convolutions

  • Large C1 values (> 0.5) can produce extreme values requiring careful handling

  • alpha < 0.5 and alpha = 1 are not currently implemented

  • N-D kernels are always non-causal

  • Spectral N-D kernels currently assume the default isotropic Fourier metric and do not support custom scale_metric

  • Does not support negative H values (use FIF_1D for H < 0)

  • C1 = 0: Routes internally to fBm_ND_circulant() for monofractal case

Analysis

Analysis module for multifractal fields

Contains methods for estimating scaling exponents and analyzing multifractal properties.

Hurst Estimation

scaleinvariance.haar_fluctuation_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False)[source]

Calculate the Hurst exponent using Haar fluctuation analysis.

For a self-similar process with Hurst exponent H, the Haar fluctuation scales as: F_H(r) ~ r^H

This function estimates H by performing a linear regression on log(F_H(r)) vs log(r).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • min_sep (int, optional) – Minimum separation/lag to include in the fit.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. Defaults to array_size - 1 (a warning is emitted for small arrays).

  • axis (int, optional) – Axis along which to compute the Haar fluctuation (default: 0).

  • return_fit (bool, optional) – If True, return the lags, fluctuation values, and fitted line.

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

    lagsnp.ndarray

    The lags used in the calculation.

    haar_valuesnp.ndarray

    The Haar fluctuation values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

scaleinvariance.structure_function_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False)[source]

Calculate the Hurst exponent using structure functions.

For a self-similar process with Hurst exponent H, the first order structure function scales as: S_1(r) ~ r^H

This function estimates H by performing a linear regression on log(S_1(r)) vs log(r).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • min_sep (int, optional) – Minimum separation/lag to include in the fit. If None, defaults to 1.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. If None, defaults to array_size/8 or array_size if array_size < 128 (with warning).

  • axis (int, optional) – Axis along which to compute the structure function (default: 0).

  • return_fit (bool, optional) – If True, return the lags, structure function values, and fitted line.

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

    lagsnp.ndarray

    The lags used in the calculation.

    sf_valuesnp.ndarray

    The structure function values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

scaleinvariance.spectral_hurst(data, max_wavelength=None, min_wavelength=None, nbins=50, axis=0, return_fit=False)[source]

Estimate the Hurst exponent from the power spectral density along specified axis.

For fractional Brownian motion, PSD ~ f^(-β) with β = 2H + 1. This function estimates H by performing linear regression on log(PSD) vs log(frequency).

Parameters:

  • data (array-like) – N-dimensional input data array

  • max_wavelength (float, optional) – Maximum wavelength to include in analysis. If None, defaults to domain_size/8 or domain_size if domain_size < 512 (with warning).

  • min_wavelength (float, optional) – Minimum wavelength to include in analysis. If None, defaults to 4*Nyquist or Nyquist if domain_size < 512.

  • nbins (int, optional) – Number of logarithmic bins for PSD (default: 50)

  • axis (int, optional) – Axis along which to compute FFT (default: 0)

  • return_fit (bool, optional) – If True, return frequencies, PSD values, and fitted line

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

    freqsnp.ndarray

    The frequencies used in the calculation

    psd_valuesnp.ndarray

    The power spectral density values

    fit_linenp.ndarray

    The fitted line used to estimate H

Analysis Functions

scaleinvariance.analysis.haar_fluctuation.haar_fluctuation_analysis(data, order=1, max_sep=None, axis=0, lags='powers of 1.2', nan_behavior='raise')[source]

Calculate the mean absolute Haar fluctuation for regularly spaced scalar data along a given axis.

Uses scipy.signal.convolve for efficient computation.

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • order (int, optional) – Order of the fluctuation (default: 1).

  • max_sep (int, optional) – Maximum separation/lag to compute. If None, defaults to (size along axis) - 1.

  • axis (int, optional) – Axis along which to compute the Haar fluctuation (default: 0).

  • lags (str, list, or np.ndarray, optional) –

    Option for selecting lag values. Can be:

    • A list or array of lags.

    • ’all’: all integer lags from 1 to max_sep.

    • ’powers of X’: lags as integers of powers of X, where X is any number (e.g., ‘powers of 2’, ‘powers of 1.2’).

  • nan_behavior (str, optional) –

    Behavior when encountering NaN values in the data. Options are:

    • ’raise’: Raise an error if NaN values are found (default).

    • ’ignore’: Ignore NaN values in the calculation.

Returns:

lags : 1-D array of lag values. haar : 1-D array of mean absolute haar fluctuation values corresponding to each lag.

Return type:

tuple (np.ndarray, np.ndarray)

scaleinvariance.analysis.haar_fluctuation.haar_fluctuation_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False)[source]

Calculate the Hurst exponent using Haar fluctuation analysis.

For a self-similar process with Hurst exponent H, the Haar fluctuation scales as: F_H(r) ~ r^H

This function estimates H by performing a linear regression on log(F_H(r)) vs log(r).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • min_sep (int, optional) – Minimum separation/lag to include in the fit.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. Defaults to array_size - 1 (a warning is emitted for small arrays).

  • axis (int, optional) – Axis along which to compute the Haar fluctuation (default: 0).

  • return_fit (bool, optional) – If True, return the lags, fluctuation values, and fitted line.

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

    lagsnp.ndarray

    The lags used in the calculation.

    haar_valuesnp.ndarray

    The Haar fluctuation values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

scaleinvariance.analysis.structure_function.structure_function_analysis(data, order=1, max_sep=None, axis=0, lags='powers of 1.2')[source]

Calculate the mean structure function for regularly spaced scalar data along a given axis.

The structure function S_n(r) of order n is defined as:

S_n(r) = <|u(x + r) - u(x)|^n>

where <> denotes the average over all valid pairs separated by lag r.

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • order (int, optional) – Order of the structure function (default: 1). Must be positive.

  • max_sep (int, optional) – Maximum lag to compute. If None, defaults to (size along axis) - 1.

  • axis (int, optional) – Axis along which to compute the structure function (default: 0).

  • lags (str, list, or np.ndarray, optional) –

    Option for selecting lag values. Can be:

    • A list or array of lags.

    • ’all’: all integer lags from 1 to max_sep.

    • ’powers of X’: lags as integers of powers of X, where X is any number (e.g., ‘powers of 2’, ‘powers of 1.2’).

Returns:

lags : 1-D array of lag values. sf : 1-D array of structure function values corresponding to each lag.

Return type:

tuple (np.ndarray, np.ndarray)

scaleinvariance.analysis.structure_function.structure_function_hurst(data, min_sep=None, max_sep=None, axis=0, return_fit=False)[source]

Calculate the Hurst exponent using structure functions.

For a self-similar process with Hurst exponent H, the first order structure function scales as: S_1(r) ~ r^H

This function estimates H by performing a linear regression on log(S_1(r)) vs log(r).

Parameters:

  • data (array-like) – N-dimensional scalar data array.

  • min_sep (int, optional) – Minimum separation/lag to include in the fit. If None, defaults to 1.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit. If None, defaults to array_size/8 or array_size if array_size < 128 (with warning).

  • axis (int, optional) – Axis along which to compute the structure function (default: 0).

  • return_fit (bool, optional) – If True, return the lags, structure function values, and fitted line.

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent.

    uncertaintyfloat

    Standard error of the Hurst exponent estimate.

    lagsnp.ndarray

    The lags used in the calculation.

    sf_valuesnp.ndarray

    The structure function values.

    fit_linenp.ndarray

    The fitted line used to estimate H.

scaleinvariance.analysis.spectral.spectral_analysis(data, max_wavelength=None, min_wavelength=None, nbins=50, axis=0)[source]

Compute binned power spectral density for data along specified axis.

Parameters:

  • data (array-like) – N-dimensional input data array

  • max_wavelength (float, optional) – Maximum wavelength to include in analysis

  • min_wavelength (float, optional) – Minimum wavelength to include in analysis

  • nbins (int, optional) – Number of logarithmic bins for PSD (default: 50)

  • axis (int, optional) – Axis along which to compute FFT (default: 0)

Returns:

binned_freq : 1-D array of binned frequencies binned_psd : 1-D array of binned power spectral density values

Return type:

tuple (np.ndarray, np.ndarray)

scaleinvariance.analysis.spectral.spectral_hurst(data, max_wavelength=None, min_wavelength=None, nbins=50, axis=0, return_fit=False)[source]

Estimate the Hurst exponent from the power spectral density along specified axis.

For fractional Brownian motion, PSD ~ f^(-β) with β = 2H + 1. This function estimates H by performing linear regression on log(PSD) vs log(frequency).

Parameters:

  • data (array-like) – N-dimensional input data array

  • max_wavelength (float, optional) – Maximum wavelength to include in analysis. If None, defaults to domain_size/8 or domain_size if domain_size < 512 (with warning).

  • min_wavelength (float, optional) – Minimum wavelength to include in analysis. If None, defaults to 4*Nyquist or Nyquist if domain_size < 512.

  • nbins (int, optional) – Number of logarithmic bins for PSD (default: 50)

  • axis (int, optional) – Axis along which to compute FFT (default: 0)

  • return_fit (bool, optional) – If True, return frequencies, PSD values, and fitted line

Returns:

  • If return_fit is False –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

  • If return_fit is True –

    Hfloat

    Estimated Hurst exponent

    uncertaintyfloat

    Standard error of the Hurst exponent estimate

    freqsnp.ndarray

    The frequencies used in the calculation

    psd_valuesnp.ndarray

    The power spectral density values

    fit_linenp.ndarray

    The fitted line used to estimate H

Utilities

Utility functions for scaleinvariance package

Common utilities used across analysis and simulation modules.

scaleinvariance.utils.estimate_hurst_from_scaling(lags, values, min_sep=None, max_sep=None, return_fit=False)[source]

Estimate Hurst exponent from scaling relationship.

Performs linear regression on log(values) vs log(lags) to estimate scaling exponent.

Parameters:

  • lags (array-like) – Lag values.

  • values (array-like) – Corresponding values (structure function or fluctuation values).

  • min_sep (int, optional) – Minimum separation/lag to include in the fit.

  • max_sep (int, optional) – Maximum separation/lag to include in the fit.

  • return_fit (bool, optional) – If True, return additional fit information.

Returns:

  • H (float) – Estimated Hurst exponent.

  • uncertainty (float) – Standard error of the Hurst exponent estimate.

  • lags_fit (np.ndarray (if return_fit=True)) – The lags used in the fit.

  • values_fit (np.ndarray (if return_fit=True)) – The values used in the fit.

  • fit_line (np.ndarray (if return_fit=True)) – The fitted line.

scaleinvariance.utils.process_lags(lags, max_sep, even_only=False)[source]

Process lag specification into array of lag values.

Parameters:

  • lags (str, list, or np.ndarray) –

    Option for selecting lag values. Can be:

    • A list or array of lags.

    • ’all’: all integer lags from 1 to max_sep.

    • ’powers of XXXX’: lags as closest intergers to powers of XXXX, e.g. 2 (~ 3 per decade) or 1.2 (~ 8 per decade).

  • max_sep (int) – Maximum separation/lag value.

  • even_only (bool, optional) – If True, filter to only even lags (default: False).

Returns:

Array of processed lag values.

Return type:

np.ndarray