def filter_gaussian(
    shape,
    n,
    gauss_scale=0.5,
    d=1.0,
    normalize=True,
    return_weight_funcs=False,
    include_mean=True,
):
    """
    Implements a set of Gaussian bandpass filters in logarithmic frequency
    scale.

    Parameters
    ----------
    shape: int or tuple
        The dimensions (height, width) of the input field. If shape is an int,
        the domain is assumed to have square shape.
    n: int
        The number of frequency bands to use. Must be greater than 2.
    gauss_scale: float
        Optional scaling parameter. Proportional to the standard deviation of
        the Gaussian weight functions.
    d: scalar, optional
        Sample spacing (inverse of the sampling rate). Defaults to 1.
    normalize: bool
        If True, normalize the weights so that for any given wavenumber
        they sum to one.
    return_weight_funcs: bool
        If True, add callable weight functions to the output dictionary with
        the key 'weight_funcs'.
    include_mean: bool
        If True, include the first Fourier wavenumber (corresponding to the
        field mean) to the first filter.

    Returns
    -------
    out: dict
        A dictionary containing the bandpass filters corresponding to the
        specified frequency bands.

    References
    ----------
    :cite:`PCH2018`

    """
    if n < 3:
        raise ValueError("n must be greater than 2")

    try:
        height, width = shape
    except TypeError:
        height, width = (shape, shape)

    max_length = max(width, height)

    rx = np.s_[: int(width / 2) + 1]

    if (height % 2) == 1:
        ry = np.s_[-int(height / 2) : int(height / 2) + 1]
    else:
        ry = np.s_[-int(height / 2) : int(height / 2)]

    y_grid, x_grid = np.ogrid[ry, rx]
    dy = int(height / 2) if height % 2 == 0 else int(height / 2) + 1

    r_2d = np.roll(np.sqrt(x_grid * x_grid + y_grid * y_grid), dy, axis=0)

    r_max = int(max_length / 2) + 1
    r_1d = np.arange(r_max)

    wfs, central_wavenumbers = _gaussweights_1d(
        max_length,
        n,
        gauss_scale=gauss_scale,
    )

    weights_1d = np.empty((n, r_max))
    weights_2d = np.empty((n, height, int(width / 2) + 1))

    for i, wf in enumerate(wfs):
        weights_1d[i, :] = wf(r_1d)
        weights_2d[i, :, :] = wf(r_2d)

    if normalize:
        weights_1d_sum = np.sum(weights_1d, axis=0)
        weights_2d_sum = np.sum(weights_2d, axis=0)
        for k in range(weights_2d.shape[0]):
            weights_1d[k, :] /= weights_1d_sum
            weights_2d[k, :, :] /= weights_2d_sum

    for i in range(len(wfs)):
        if i == 0 and include_mean:
            weights_1d[i, 0] = 1.0
            weights_2d[i, 0, 0] = 1.0
        else:
            weights_1d[i, 0] = 0.0
            weights_2d[i, 0, 0] = 0.0

    out = {"weights_1d": weights_1d, "weights_2d": weights_2d}
    out["shape"] = shape

    central_wavenumbers = np.array(central_wavenumbers)
    out["central_wavenumbers"] = central_wavenumbers

    # Compute frequencies
    central_freqs = 1.0 * central_wavenumbers / max_length
    central_freqs[0] = 1.0 / max_length
    central_freqs[-1] = 0.5  # Nyquist freq
    central_freqs = 1.0 * d * central_freqs
    out["central_freqs"] = central_freqs

    if return_weight_funcs:
        out["weight_funcs"] = wfs

    return out
Implements a set of Gaussian bandpass filters in logarithmic frequency
scale.

Parameters
----------
shape: int or tuple
The dimensions (height, width) of the input field. If shape is an int,
the domain is assumed to have square shape.
n: int
The number of frequency bands to use. Must be greater than 2.
gauss_scale: float
Optional scaling parameter. Proportional to the standard deviation of
the Gaussian weight functions.
d: scalar, optional
Sample spacing (inverse of the sampling rate). Defaults to 1.
normalize: bool
If True, normalize the weights so that for any given wavenumber
they sum to one.
return_weight_funcs: bool
If True, add callable weight functions to the output dictionary with
the key 'weight_funcs'.
include_mean: bool
If True, include the first Fourier wavenumber (corresponding to the
field mean) to the first filter.

Returns
-------
out: dict
A dictionary containing the bandpass filters corresponding to the
specified frequency bands.

References
----------
:cite:`PCH2018`
pysteps/cascade/bandpass_filters.py