pysmo.tools.signal

@register_filter
def bandpass[T: Seismogram](
    seismogram: T,
    freqmin: float = 0.1,
    freqmax: float = 0.5,
    corners: int = 2,
    zerophase: bool = False,
    clone: bool = False,
) -> T | None:
    """
    Apply a bandpass filter to the input seismogram.

    Args:
        seismogram: The input seismogram to be filtered.
        freqmin: The minimum frequency of the bandpass filter (in Hz).
        freqmax: The maximum frequency of the bandpass filter (in Hz).
        corners: The number of corners (poles) for the Butterworth filter.
        zerophase: If True, apply the filter in both forward and reverse directions to achieve zero phase distortion.
        clone: If True, return a new Seismogram object with the filtered data. If False, modify the input seismogram in place.

    Returns:
        A new Seismogram object containing the filtered data when called with `clone=True`.
    """
    fe = 0.5 / seismogram.delta.total_seconds()
    low = freqmin / fe
    high = freqmax / fe

    if not (0 < low < 1):
        raise ValueError(
            f"freqmin ({freqmin}) is invalid for sampling rate {1 / seismogram.delta.total_seconds()} Hz."
        )
    if not (0 < high < 1):
        raise ValueError(
            f"freqmax ({freqmax}) is invalid for sampling rate {1 / seismogram.delta.total_seconds()} Hz."
        )
    if freqmin >= freqmax:
        raise ValueError("freqmin must be less than freqmax.")

    sos = iirfilter(corners, [low, high], btype="band", ftype="butter", output="sos")

    if clone:
        seismogram = deepcopy(seismogram)

    if zerophase:
        seismogram.data = sosfiltfilt(sos, seismogram.data)
    else:
        seismogram.data = sosfilt(sos, seismogram.data)

    return seismogram if clone else None

@register_filter
def bandstop[T: Seismogram](
    seismogram: T,
    freqmin: float = 0.1,
    freqmax: float = 0.5,
    corners: int = 2,
    zerophase: bool = False,
    clone: bool = False,
) -> T | None:
    """
    Apply a bandstop filter to the input seismogram.

    Args:
        seismogram: The input seismogram to be filtered.
        freqmin: The minimum frequency of the bandstop filter (in Hz).
        freqmax: The maximum frequency of the bandstop filter (in Hz).
        corners: The number of corners (poles) for the Butterworth filter.
        zerophase: If True, apply the filter in both forward and reverse directions to achieve zero phase distortion.
        clone: If True, return a new Seismogram object with the filtered data. If False, modify the input seismogram in place.

    Returns:
        A new Seismogram object containing the filtered data when called with `clone=True`.
    """
    fe = 0.5 / seismogram.delta.total_seconds()
    low = freqmin / fe
    high = freqmax / fe

    if not (0 < low < 1):
        raise ValueError(
            f"freqmin ({freqmin}) is invalid for sampling rate {1 / seismogram.delta.total_seconds()} Hz."
        )
    if not (0 < high < 1):
        raise ValueError(
            f"freqmax ({freqmax}) is invalid for sampling rate {1 / seismogram.delta.total_seconds()} Hz."
        )
    if freqmin >= freqmax:
        raise ValueError("freqmin must be less than freqmax.")

    sos = iirfilter(
        corners, [low, high], btype="bandstop", ftype="butter", output="sos"
    )

    if clone:
        seismogram = deepcopy(seismogram)

    if zerophase:
        seismogram.data = sosfiltfilt(sos, seismogram.data)
    else:
        seismogram.data = sosfilt(sos, seismogram.data)

    return seismogram if clone else None

def delay(
    seismogram1: Seismogram,
    seismogram2: Seismogram,
    total_delay: bool = False,
    max_shift: Timedelta | None = None,
    abs_max: bool = False,
) -> tuple[Timedelta, float]:
    """
    Cross correlates two seismograms to determine signal delay.

    This function is a wrapper around the
    [`scipy.signal.correlate`][scipy.signal.correlate] function. The default
    behaviour is to call the correlate function with `#!py mode="full"` using
    the full length data of the input seismograms. This is the most robust
    option, but also the slowest.

    If an approximate delay is known (e.g. because a particular phase is being
    targeted using a computed arrival time), the search space can be limited
    by setting the `max_shift` parameter to a value. The length of the
    seismogram data used for the cross-correlation is then set such that the
    calculated delay lies within +/- `max_shift`.

    Note:
        `max_shift` intentionally does *not* take the begin times of the
        seismograms into consideration. Thus, calling `#!py delay()` with
        `#!py total_delay=True` may return a delay that is larger than
        `max_shift`.

    Implications of setting the `max_shift` parameter are as follows:

    - This mode requires the seismograms to be of equal length.
    - If the true delay (i.e. the amount of time the seismograms _should_ be
      shifted by) lies within the `max_shift` range, and also produces the
      highest correlation, the delay time returned is identical for both
      methods.
    - If the true delay lies outside the `max_shift` range and produces the
      highest correlation, the delay time returned will be incorrect when
      `max_shift` is set.
    - In the event that the true delay lies within the `max_shift` range but
      the maximum signal correlation occurs outside, it will be correctly
      retrieved when the `max_shift` parameter is set, while not setting it
      yields an incorrect result.


    Args:
        seismogram1: First seismogram to use for cross correlation.
        seismogram2: Second seismogram to use for cross correlation.
        total_delay: Include the difference in `begin_time` in the delay.
        max_shift: Maximum (absolute) length of the delay.
        abs_max: Return the delay corresponding to absolute maximum.

    Returns:
        delay: Time delay of the second seismogram with respect to the first.
        ccnorm: Normalised cross correlation value of the overlapping
            seismograms *after* shifting (uses
            [`scipy.stats.mstats.pearsonr`][scipy.stats.mstats.pearsonr] for
            the calculation). This value ranges from -1 to 1, with 1 indicating
            a perfect correlation, 0 indicating no correlation, and -1
            indicating a perfect anti-correlation.

    Examples:
        To illustration the use of the `delay()` function, we read a seismogram
        from a SAC file and then generate a second seismogram from it by
        modifying it with a shift in the data and a shift in the begin time.
        This means the true delay is known, and we are able to compare the
        computed delays with this known delay:

        ```python
        >>> from pysmo import MiniSeismogram
        >>> from pysmo.classes import SAC
        >>> from pysmo.functions import detrend, clone_to_mini
        >>> from pysmo.tools.signal import delay
        >>> from datetime import timedelta
        >>> import numpy as np
        >>>
        >>> # Create a Seismogram from a SAC file and detrend it:
        >>> seis1 = SAC.from_file("example.sac").seismogram
        >>> detrend(seis1)
        >>>
        >>> # Create a second seismogram from the first with
        >>> # a different begin_time and a shift in the data:
        >>> seis2 = clone_to_mini(MiniSeismogram, seis1)
        >>> nroll = 1234
        >>> seis2.data = np.roll(seis2.data, nroll)
        >>> begin_time_delay = timedelta(seconds=100)
        >>> seis2.begin_time += begin_time_delay
        >>>
        >>> # The signal delay is the number of samples shifted * delta:
        >>> (signal_delay := nroll * seis1.delta).total_seconds()
        24.68
        >>>
        >>> # Call the delay function with the two seismograms and verify
        >>> # that the caclulated_delay is equal to the known signal delay:
        >>> calculated_delay, _ = delay(seis1, seis2)
        >>> calculated_delay == signal_delay
        True
        >>>
        ```

        As we know what the true delay is, we can mimic a scenario where an
        approximate delay is known prior to the cross-correlation, which can be
        used to limit the search space and speed up the calculation. Here we
        assign a value of the known signal delay plus 1 second to the
        `max_shift` parameter:

        ```python
        >>> max_shift = signal_delay + timedelta(seconds=1)
        >>> calculated_delay, _ = delay(seis1, seis2, max_shift=max_shift)
        >>> # As before, the calculated delay should be equal to the signal delay:
        >>> calculated_delay == signal_delay
        True
        >>>
        ```

        Setting `total_delay=True`, we also takes into account the difference
        in `begin_time` between the two seismograms:

        ```python
        >>> calculated_delay, _ = delay(seis1, seis2, total_delay=True, max_shift=signal_delay+timedelta(seconds=1))
        >>> # With `total_delay=True`, the calculated delay should be equal to
        >>> # the signal delay plus the begin time difference:
        >>> calculated_delay == signal_delay + (seis1.begin_time - seis2.begin_time)
        True
        >>>
        ```

        To demonstrate usage of the `abs_max` parameter, we flip the sign of
        the second seismogram data:

        ```python
        >>> seis2.data *= -1
        >>> calculated_delay, ccnorm = delay(seis1, seis2)
        >>> # Without `abs_max=True`, the calculated delay is no longer equal
        >>> # to the true signal delay (as expected):
        >>> calculated_delay == signal_delay
        False
        >>> # The normalised cross correlation value is also not very high
        >>> ccnorm
        np.float64(0.4267205)
        >>>
        >>> calculated_delay, ccnorm = delay(seis1, seis2, abs_max=True)
        >>> # with `abs_max=True`, the signal delay is again retrieved:
        >>> calculated_delay == signal_delay
        True
        >>> # And, as the signals are completely opposite, the normalised
        >>> # cross correlation value is -1:
        >>> np.testing.assert_approx_equal(ccnorm, -1)
        >>>
        ```
    """

    _check_same_delta(seismogram1, seismogram2)

    if max_shift is not None and len(seismogram1.data) != len(seismogram2.data):
        raise ValueError(
            "Input seismograms must be of equal length when using `max_shift`."
        )

    data1, data2 = seismogram1.data, seismogram2.data
    delta = seismogram1.delta

    if max_shift is not None:
        max_lag_in_samples = math.ceil(max_shift / delta)
        data1 = np.pad(data1, max_lag_in_samples)
        corr = _correlate(data1, data2, mode="valid")
    else:
        corr = _correlate(data1, data2, mode="full")

    corr_index = np.argmax(corr)

    if abs_max and np.max(corr) < -1 * np.min(corr):
        corr_index = np.argmin(corr)

    if max_shift is not None:
        shift = -int(corr_index - max_lag_in_samples)
    else:
        shift = int(len(data2) - 1 - corr_index)

    delay = shift * delta

    # find overlapping parts of seismograms after alignment
    if shift < 0:
        data1 = data1[-shift:]
    else:
        data2 = data2[shift:]
    if len(data1) > len(data2):
        data1 = data1[: len(data2)]
    else:
        data2 = data2[: len(data1)]

    covr, _ = _pearsonr(data1, data2)

    if total_delay:
        delay += seismogram1.begin_time - seismogram2.begin_time

    return delay, covr

@register_filter
def envelope[T: Seismogram](
    seismogram: T, Tn: float, alpha: float, clone: bool = False
) -> T | None:
    """
    Calculates the envelope of a gaussian filtered seismogram.

    Args:
        seismogram: Name of the seismogram object passed to this function.
        Tn: Center period of Gaussian filter [in seconds]
        alpha: Set alpha (which determines filterwidth)
        clone: If True, return a new Seismogram object with the filtered data. If False, modify the input seismogram in place.

    Returns:
        Seismogram containing the envelope

    Examples:
        ```python
        >>> from pysmo.classes import SAC
        >>> from pysmo.tools.signal import envelope
        >>> seis = SAC.from_file("example.sac").seismogram
        >>> Tn = 50 # Center Gaussian filter at 50s period
        >>> alpha = 50 # Set alpha (which determines filterwidth) to 50
        >>> envelope_seis = envelope(seis, Tn, alpha, clone=True)
        >>>
        ```
    """
    if clone:
        seismogram = deepcopy(seismogram)
    seismogram.data = _gauss(seismogram, Tn, alpha)[0]
    return seismogram if clone else None

def filter[T: Seismogram](
    seismogram: T,
    filter_name: FilterName,
    clone: bool = False,
    **filter_options: float | int | bool,
) -> T | None:
    """
    Apply a specified filter to the input seismogram.

    This function is a convenience wrapper that calls other filters in this module.

    Args:
        seismogram: The input seismogram to be filtered.
        filter_name: The type of filter to apply.
        clone: If True, return a new Seismogram object with the filtered data. If False, modify the input seismogram in place.
        **filter_options: Filter parameters passed to the specified filter function.

    Returns:
        A new Seismogram object containing the filtered data when called with `clone=True`.

    Examples:
        ```python
        >>> from pysmo.classes import SAC
        >>> from pysmo.tools.signal import filter
        >>> seis = SAC.from_file("example.sac").seismogram
        >>>
        >>> # create a new filtered seismogram with a lowpass filter
        >>> filtered_seis = filter(seis, "lowpass", freqmax=0.5, clone=True)
        >>>
        >>> # or update in place with a bandpass filter
        >>> filter(seis, "bandpass", freqmin=0.1, freqmax=0.5)
        >>>
        ```
    """

    try:
        filter_func = _FILTER_REGISTRY[filter_name]
    except KeyError:
        # This fallback handles cases where FilterName is updated but
        # the function isn't decorated yet.
        valid_filters = ", ".join(_FILTER_REGISTRY.keys())
        raise ValueError(
            f"Filter '{filter_name}' is not registered. Available: {valid_filters}"
        )

    if clone:
        return filter_func(seismogram, clone=True, **filter_options)
    filter_func(seismogram, clone=False, **filter_options)
    return None

@register_filter
def gauss[T: Seismogram](
    seismogram: T, Tn: float, alpha: float, clone: bool = False
) -> T | None:
    """
    Returns a gaussian filtered seismogram.

    Args:
        seismogram: Name of the SAC object passed to this function.
        Tn: Center period of Gaussian filter [in seconds]
        alpha: Set alpha (which determines filterwidth)

    Returns:
        Gaussian filtered seismogram.

    Examples:
        ```python
        >>> from pysmo.classes import SAC
        >>> from pysmo.tools.signal import gauss
        >>> seis = SAC.from_file("example.sac").seismogram
        >>> Tn = 50 # Center Gaussian filter at 50s period
        >>> alpha = 50 # Set alpha (which determines filterwidth) to 50
        >>> gauss_seis = gauss(seis, Tn, alpha, clone=True)
        >>>
        ```
    """
    if clone:
        seismogram = deepcopy(seismogram)
    seismogram.data = _gauss(seismogram, Tn, alpha)[1]
    return seismogram if clone else None

@register_filter
def highpass[T: Seismogram](
    seismogram: T,
    freqmin: float = 0.1,
    corners: int = 2,
    zerophase: bool = False,
    clone: bool = False,
) -> T | None:
    """
    Apply a highpass filter to the input seismogram.

    Args:
        seismogram: The input seismogram to be filtered.
        freqmin: The minimum frequency of the highpass filter (in Hz).
        corners: The number of corners (poles) for the Butterworth filter.
        zerophase: If True, apply the filter in both forward and reverse directions to achieve zero phase distortion.
        clone: If True, return a new Seismogram object with the filtered data. If False, modify the input seismogram in place.

    Returns:
        A new Seismogram object containing the filtered data when called with `clone=True`.
    """
    fe = 0.5 / seismogram.delta.total_seconds()
    low = freqmin / fe

    if not (0 < low < 1):
        raise ValueError(
            f"freqmin ({freqmin}) is invalid for sampling rate {1 / seismogram.delta.total_seconds()} Hz."
        )

    sos = iirfilter(corners, low, btype="high", ftype="butter", output="sos")

    if clone:
        seismogram = deepcopy(seismogram)

    if zerophase:
        seismogram.data = sosfiltfilt(sos, seismogram.data)
    else:
        seismogram.data = sosfilt(sos, seismogram.data)

    return seismogram if clone else None

@register_filter
def lowpass[T: Seismogram](
    seismogram: T,
    freqmax: float = 0.5,
    corners: int = 2,
    zerophase: bool = False,
    clone: bool = False,
) -> T | None:
    """
    Apply a lowpass filter to the input seismogram.

    Args:
        seismogram: The input seismogram to be filtered.
        freqmax: The maximum frequency of the lowpass filter (in Hz).
        corners: The number of corners (poles) for the Butterworth filter.
        zerophase: If True, apply the filter in both forward and reverse directions to achieve zero phase distortion.
        clone: If True, return a new Seismogram object with the filtered data. If False, modify the input seismogram in place.

    Returns:
        A new Seismogram object containing the filtered data when called with `clone=True`.
    """
    fe = 0.5 / seismogram.delta.total_seconds()
    high = freqmax / fe

    if not (0 < high < 1):
        raise ValueError(
            f"freqmax ({freqmax}) is invalid for sampling rate {1 / seismogram.delta.total_seconds()} Hz."
        )

    sos = iirfilter(corners, high, btype="low", ftype="butter", output="sos")

    if clone:
        seismogram = deepcopy(seismogram)

    if zerophase:
        seismogram.data = sosfiltfilt(sos, seismogram.data)
    else:
        seismogram.data = sosfilt(sos, seismogram.data)

    return seismogram if clone else None

def mccc(
    seismograms: Sequence[Seismogram], min_cc: float = 0.5, damping: float = 0.1
) -> tuple[list[Timedelta], list[Timedelta], Timedelta]:
    """
    Multi-Channel Cross-Correlation (MCCC) for relative arrival times.

    Computes all pairwise cross-correlation delays using
    [`multi_multi_delay`][pysmo.tools.signal.multi_multi_delay], then solves
    for self-consistent relative time shifts using a weighted least-squares
    inversion with a zero-mean constraint and Tikhonov regularization. Pairs
    whose correlation coefficient falls below `min_cc` are excluded from the
    inversion.

    The returned `times` list sums to zero by construction, so the values
    represent relative shifts around the group mean.

    Args:
        seismograms: Sequence of Seismogram objects. All must share the same
            sampling interval.
        min_cc: Minimum correlation coefficient required to include a pair
            in the inversion.
        damping: Tikhonov regularization strength. Set to 0 to disable.

    Returns:
        times: List of relative arrival times.
        errors: List of standard errors.
        rmse: Root-mean-square error of the fit.

    Raises:
        ValueError: If any seismogram has a different sampling rate than the
            others (raised by `multi_multi_delay`).

    Examples:
        Create seismograms with known shifts and recover them with `mccc`:

        ```python
        >>> from pysmo import MiniSeismogram
        >>> from pysmo.tools.signal import mccc
        >>> import numpy as np
        >>>
        >>> # Build three seismograms with known shifts (in samples):
        >>> data = np.sin(np.linspace(0, 8 * np.pi, 1000))
        >>> shifts = [0, 5, -10]
        >>> seismograms = [MiniSeismogram(data=np.roll(data, s)) for s in shifts]
        >>>
        >>> # Run MCCC inversion:
        >>> times, errors, rmse = mccc(seismograms)
        >>>
        >>> # The relative times sum to approximately zero:
        >>> abs(sum(t.total_seconds() for t in times)) < 1e-5
        True
        >>>
        >>> # Pairwise differences recover the known shifts
        >>> # (times[i] - times[j] ≈ (shifts[j] - shifts[i]) * delta):
        >>> round((times[1] - times[0]).total_seconds())
        -5
        >>> round((times[2] - times[0]).total_seconds())
        10
        >>>
        ```

    References:
        VanDecar, J. C., and R. S. Crosson. “Determination of Teleseismic
        Relative Phase Arrival Times Using Multi-Channel Cross-Correlation and
        Least Squares.” Bulletin of the Seismological Society of America,
        vol. 80, no. 1, Feb. 1990, pp. 150–69,
        <https://doi.org/10.1785/BSSA0800010150>.
    """
    n_traces = len(seismograms)
    zero_delta = Timedelta(0)

    if n_traces < 2:
        return [zero_delta] * n_traces, [zero_delta] * n_traces, zero_delta

    delay_matrix, coeff_matrix = multi_multi_delay(seismograms)

    # Convert numpy.timedelta64 array to float seconds for linear algebra
    delay_matrix_seconds = delay_matrix / np.timedelta64(1, "s")

    # Build linear system (g @ m = d)
    rows: list[np.ndarray] = []
    data_vec: list[float] = []
    weights: list[float] = []

    for i in range(n_traces):
        for j in range(i + 1, n_traces):
            cc = coeff_matrix[i, j]
            if cc < min_cc:
                continue

            lag_seconds = delay_matrix_seconds[i, j]

            row = np.zeros(n_traces)
            row[i] = 1.0
            row[j] = -1.0

            rows.append(row)
            data_vec.append(lag_seconds)
            weights.append(cc**2)

    if not rows:
        return [zero_delta] * n_traces, [zero_delta] * n_traces, zero_delta

    g = np.array(rows)
    d = np.array(data_vec)
    w = np.array(weights)

    # Apply weights
    g_weighted = g * w[:, np.newaxis]
    d_weighted = d * w

    # Zero-mean constraint (sum of times = 0)
    constraint_weight = np.sum(w)
    g_system = np.vstack([g_weighted, np.ones(n_traces) * constraint_weight])
    d_system = np.concatenate([d_weighted, [0.0]])

    # Tikhonov regularization
    if damping > 0:
        g_system = np.vstack([g_system, damping * np.eye(n_traces)])
        d_system = np.concatenate([d_system, np.zeros(n_traces)])

    # Solve least squares
    solution, _, _, _ = lstsq(g_system, d_system)

    # Compute statistics
    predicted = g @ solution
    residuals = d - predicted
    sse = np.sum((residuals**2) * w)
    dof = max(len(d) - n_traces, 1)
    sigma_squared = sse / dof

    try:
        cov_matrix = sigma_squared * inv(g_system.T @ g_system)
        std_errors = np.sqrt(np.abs(np.diag(cov_matrix)))
    except np.linalg.LinAlgError:
        std_errors = np.zeros(n_traces)

    times = [Timedelta(seconds=float(t)) for t in solution]
    errors = [Timedelta(seconds=float(e)) for e in std_errors]
    rmse = Timedelta(seconds=float(np.sqrt(sse / len(d))))

    return times, errors, rmse

def multi_delay(
    template: Seismogram, seismograms: Sequence[Seismogram], abs_max: bool = False
) -> tuple[list[Timedelta], list[float]]:
    """
    Calculates delays and correlation coefficients for a list of seismograms against a template.

    This function uses FFT-based cross-correlation to efficiently compute delays
    for multiple seismograms at once against a single template. This is faster
    than calling [`delay`][pysmo.tools.signal.delay] in a loop, as the template
    FFT is computed only once.

    Args:
        template: Template seismogram object.
        seismograms: Sequence of Seismogram objects.
        abs_max: If `True`, uses absolute max correlation (polarity insensitive).

    Returns:
        delays: Delays of each input seismogram relative to template.
        coeffs: Correlation coefficients at maximum correlation for each seismogram.

    Raises:
        ValueError: If any seismogram has a different sampling rate than the template.

    Examples:
        Create a template seismogram and several shifted copies, then use
        `multi_delay` to recover the known shifts:

        ```python
        >>> from pysmo import MiniSeismogram
        >>> from pysmo.tools.signal import multi_delay
        >>> import numpy as np
        >>>
        >>> # Create a template seismogram with sinusoidal data:
        >>> data = np.sin(np.linspace(0, 8 * np.pi, 1000))
        >>> template = MiniSeismogram(data=data)
        >>>
        >>> # Create shifted copies (shifts in samples):
        >>> shifts = [0, 10, -5]
        >>> seismograms = [MiniSeismogram(data=np.roll(data, s)) for s in shifts]
        >>>
        >>> # Calculate delays for all seismograms at once:
        >>> delays, coeffs = multi_delay(template, seismograms)
        >>> [d.total_seconds() for d in delays]
        [0.0, 10.0, -5.0]
        >>>
        ```

        Use `abs_max=True` for polarity-insensitive matching:

        ```python
        >>> flipped = MiniSeismogram(data=-np.roll(data, 10))
        >>> delays, coeffs = multi_delay(template, [flipped], abs_max=True)
        >>> delays[0].total_seconds()
        10.0
        >>> coeffs[0] < 0
        True
        >>>
        ```
    """
    if not seismograms:
        return [], []

    _check_same_delta(template, seismograms)

    # setup dimensions & FFT length
    n_traces = len(seismograms)
    len_t = len(template.data)
    len_s = max(len(s.data) for s in seismograms)

    # pad to avoid circular convolution artifacts
    # (length >= len_template + len_signal - 1)
    n_fft = next_fast_len(len_s + len_t - 1)

    # pad and normalise template
    t_data = template.data
    t_mean = np.mean(t_data)
    t_std = np.std(t_data)
    if t_std == 0:
        t_std = 1.0
    template_padded = np.zeros(n_fft, dtype=float)
    template_padded[:len_t] = (t_data - t_mean) / t_std

    # normalise *before* padding to keep stats valid
    seismogram_matrix = np.zeros((n_traces, n_fft), dtype=float)
    for i, s in enumerate(seismograms):
        data = s.data
        mean = np.mean(data)
        std = np.std(data)
        if std == 0:
            std = 1.0
        seismogram_matrix[i, : len(data)] = (data - mean) / std

    # forward FFT (rfft is faster for real data)
    t_freq = rfft(template_padded, n=n_fft)
    s_freq = rfft(seismogram_matrix, n=n_fft, axis=1)

    # cross-correlation in frequency domain
    cc_freq = s_freq * np.conj(t_freq)

    # inverse FFT & scale (div by len(template) for Pearson approx)
    cc_matrix = irfft(cc_freq, n=n_fft, axis=1) / len_t

    # find maxima
    if abs_max:
        max_indices = np.argmax(np.abs(cc_matrix), axis=1)
    else:
        max_indices = np.argmax(cc_matrix, axis=1)

    # convert circular indices to signed lags
    mid_point = n_fft // 2
    signed_lags = np.where(max_indices <= mid_point, max_indices, max_indices - n_fft)
    delta = template.delta
    delays = [int(lag) * delta for lag in signed_lags]
    coeffs: list[float] = cc_matrix[np.arange(n_traces), max_indices].tolist()

    return delays, coeffs