def temporal_autocorrelation(
    x,
    d=0,
    domain="spatial",
    x_shape=None,
    mask=None,
    use_full_fft=False,
    window="gaussian",
    window_radius=np.inf,
):
    """
    Compute lag-l temporal autocorrelation coefficients
    :math:`\gamma_l=\mbox{corr}(x(t),x(t-l))`, :math:`l=1,2,\dots,n-1`,
    from a time series :math:`x_1,x_2,\dots,x_n`. If a multivariate time series
    is given, each element of :math:`x_i` is treated as one sample from the
    process generating the time series. Use
    :py:func:`temporal_autocorrelation_multivariate` if cross-correlations
    between different elements of the time series are desired.

    Parameters
    ----------
    x: array_like
        Array of shape (n, ...), where each row contains one sample from the
        time series :math:`x_i`. The inputs are assumed to be in increasing
        order with respect to time, and the time step is assumed to be regular.
        All inputs are required to have finite values. The remaining dimensions
        after the first one are flattened before computing the correlation
        coefficients.
    d: {0,1}
        The order of differencing. If d=1, the input time series is differenced
        before computing the correlation coefficients. In this case, a time
        series of length n+1 is needed for computing the n-1 coefficients.
    domain: {"spatial", "spectral"}
        The domain of the time series x. If domain is "spectral", the elements
        of x are assumed to represent the FFTs of the original elements.
    x_shape: tuple
        The shape of the original arrays in the spatial domain before applying
        the FFT. Required if domain is "spectral".
    mask: array_like
        Optional mask to use for computing the correlation coefficients. Input
        elements with mask==False are excluded from the computations. The shape
        of the mask is expected to be x.shape[1:]. Applicable if domain is
        "spatial".
    use_full_fft: bool
        If True, x represents the full FFTs of the original arrays. Otherwise,
        the elements of x are assumed to contain only the symmetric part, i.e.
        in the format returned by numpy.fft.rfft2. Applicable if domain is
        'spectral'. Defaults to False.
    window: {"gaussian", "uniform"}
        The weight function to use for the moving window. Applicable if
        window_radius < np.inf. Defaults to 'gaussian'.
    window_radius: float
        If window_radius < np.inf, the correlation coefficients are computed in
        a moving window. Defaults to np.inf (i.e. the coefficients are computed
        over the whole domain). If window is 'gaussian', window_radius is the
        standard deviation of the Gaussian filter. If window is 'uniform', the
        size of the window is 2*window_radius+1.

    Returns
    -------
    out: list
        List of length n-1 containing the temporal autocorrelation coefficients
        :math:`\gamma_i` for time lags :math:`l=1,2,...,n-1`. If
        window_radius<np.inf, the elements of the list are arrays of shape
        x.shape[1:]. In this case, nan values are assigned, when the sample size
        for computing the correlation coefficients is too small.

    Notes
    -----
    Computation of correlation coefficients in the spectral domain is currently
    implemented only for two-dimensional fields.

    """
    if len(x.shape) < 2:
        raise ValueError("the dimension of x must be >= 2")
    if len(x.shape) != 3 and domain == "spectral":
        raise NotImplementedError(
            "len(x.shape[1:]) = %d, but with domain == 'spectral', this function has only been implemented for two-dimensional fields"
            % len(x.shape[1:])
        )
    if mask is not None and mask.shape != x.shape[1:]:
        raise ValueError(
            "dimension mismatch between x and mask: x.shape[1:]=%s, mask.shape=%s"
            % (str(x.shape[1:]), str(mask.shape))
        )
    if np.any(~np.isfinite(x)):
        raise ValueError("x contains non-finite values")

    if d == 1:
        x = np.diff(x, axis=0)

    if domain == "spatial" and mask is None:
        mask = np.ones(x.shape[1:], dtype=bool)

    gamma = []
    for k in range(x.shape[0] - 1):
        if domain == "spatial":
            if window_radius == np.inf:
                cc = np.corrcoef(x[-1, :][mask], x[-(k + 2), :][mask])[0, 1]
            else:
                cc = _moving_window_corrcoef(
                    x[-1, :], x[-(k + 2), :], window_radius, mask=mask
                )
        else:
            cc = spectral.corrcoef(
                x[-1, :, :], x[-(k + 2), :], x_shape, use_full_fft=use_full_fft
            )
        gamma.append(cc)

    return gamma
Compute lag-l temporal autocorrelation coefficients
:math:`\gamma_l=\mbox{corr}(x(t),x(t-l))`, :math:`l=1,2,\dots,n-1`,
from a time series :math:`x_1,x_2,\dots,x_n`. If a multivariate time series
is given, each element of :math:`x_i` is treated as one sample from the
process generating the time series. Use
:py:func:`temporal_autocorrelation_multivariate` if cross-correlations
between different elements of the time series are desired.

Parameters
----------
x: array_like
Array of shape (n, ...), where each row contains one sample from the
time series :math:`x_i`. The inputs are assumed to be in increasing
order with respect to time, and the time step is assumed to be regular.
All inputs are required to have finite values. The remaining dimensions
after the first one are flattened before computing the correlation
coefficients.
d: {0,1}
The order of differencing. If d=1, the input time series is differenced
before computing the correlation coefficients. In this case, a time
series of length n+1 is needed for computing the n-1 coefficients.
domain: {"spatial", "spectral"}
The domain of the time series x. If domain is "spectral", the elements
of x are assumed to represent the FFTs of the original elements.
x_shape: tuple
The shape of the original arrays in the spatial domain before applying
the FFT. Required if domain is "spectral".
mask: array_like
Optional mask to use for computing the correlation coefficients. Input
elements with mask==False are excluded from the computations. The shape
of the mask is expected to be x.shape[1:]. Applicable if domain is
"spatial".
use_full_fft: bool
If True, x represents the full FFTs of the original arrays. Otherwise,
the elements of x are assumed to contain only the symmetric part, i.e.
in the format returned by numpy.fft.rfft2. Applicable if domain is
'spectral'. Defaults to False.
window: {"gaussian", "uniform"}
The weight function to use for the moving window. Applicable if
window_radius < np.inf. Defaults to 'gaussian'.
window_radius: float
If window_radius < np.inf, the correlation coefficients are computed in
a moving window. Defaults to np.inf (i.e. the coefficients are computed
over the whole domain). If window is 'gaussian', window_radius is the
standard deviation of the Gaussian filter. If window is 'uniform', the
size of the window is 2*window_radius+1.

Returns
-------
out: list
List of length n-1 containing the temporal autocorrelation coefficients
:math:`\gamma_i` for time lags :math:`l=1,2,...,n-1`. If
window_radius<np.inf, the elements of the list are arrays of shape
x.shape[1:]. In this case, nan values are assigned, when the sample size
for computing the correlation coefficients is too small.

Notes
-----
Computation of correlation coefficients in the spectral domain is currently
implemented only for two-dimensional fields.
pysteps/timeseries/correlation.py