def downscale(
    precip,
    ds_factor,
    alpha=None,
    threshold=None,
    return_alpha=False,
    kernel_type=None,
    spectral_fusion=False,
):
    """
    Downscale a rainfall field by increasing its spatial resolution by a positive
    integer factor.

    Parameters
    ----------
    precip: array_like
        Array of shape (m, n) containing the input field.
        The input is expected to contain rain rate values.
        All values are required to be finite.
    alpha: float, optional
        Spectral slope. If None, the slope is estimated from
        the input array.
    ds_factor: positive int
        Downscaling factor, it specifies by how many times
        to increase the initial grid resolution.
    threshold: float, optional
        Set all values lower than the threshold to zero.
    return_alpha: bool, optional
        Whether to return the estimated spectral slope ``alpha``.
    kernel_type: {None, "gaussian", "uniform", "tophat"}
        The name of the smoothing operator. If None no smoothing is applied.
    spectral_fusion: bool, optional
        Whether to apply spectral merging as in :cite:`DOnofrio2014`.

    Returns
    -------
    precip_highres: ndarray
        Array of shape (m * ds_factor, n * ds_factor) containing
        the downscaled field.
    alpha: float
        Returned only when ``return_alpha=True``.

    Notes
    -----
    Currently, the pysteps implementation of RainFARM only covers spatial downscaling.
    That is, it can improve the spatial resolution of a rainfall field. However, unlike
    the original algorithm from Rebora et al. (2006), it cannot downscale the temporal
    dimension. It implements spectral merging from D'Onofrio et al. (2014).

    References
    ----------
    :cite:`Rebora2006`
    :cite:`DOnofrio2014`

    """

    # Validate inputs
    if not np.isfinite(precip).all():
        raise ValueError("All values in 'precip' must be finite.")
    if not isinstance(ds_factor, int) or ds_factor <= 0:
        raise ValueError("'ds_factor' must be a positive integer.")

    # Preprocess the input field if spectral fusion is enabled
    precip_transformed = _gaussianize(precip) if spectral_fusion else precip

    # Compute frequency arrays for the original and high-resolution fields
    freq_array = _compute_freq_array(precip_transformed)
    freq_array_highres = _compute_freq_array(precip_transformed, ds_factor)

    # Estimate spectral slope alpha if not provided
    if alpha is None:
        alpha = _estimate_alpha(precip_transformed, freq_array)

    # Generate noise field
    noise_field = _compute_noise_field(freq_array_highres, alpha)

    # Apply spectral fusion if enabled
    if spectral_fusion:
        noise_field /= noise_field.shape[0] ** 2
        noise_field = np.exp(noise_field)
        noise_field = _apply_spectral_fusion(
            precip_transformed, noise_field, freq_array, freq_array_highres, ds_factor
        )

    # Normalize and exponentiate the noise field
    noise_field /= noise_field.std()
    noise_field = np.exp(noise_field)

    # Aggregate the noise field to low resolution
    noise_lowres = aggregate_fields(noise_field, ds_factor, axis=(0, 1))

    # Expand input and noise fields to high resolution
    precip_expanded = np.kron(precip, np.ones((ds_factor, ds_factor)))
    noise_lowres_expanded = np.kron(noise_lowres, np.ones((ds_factor, ds_factor)))

    # Apply smoothing if a kernel type is provided
    if kernel_type:
        if kernel_type not in _make_kernel:
            raise ValueError(
                f"kernel type '{kernel_type}' is invalid, available kernels: {list(_make_kernel)}"
            )
        kernel = _make_kernel[kernel_type](ds_factor)
        precip_expanded = _balanced_spatial_average(precip_expanded, kernel)
        noise_lowres_expanded = _balanced_spatial_average(noise_lowres_expanded, kernel)

    # Normalize the high-res precipitation field by the low-res noise field
    norm_k0 = precip_expanded / noise_lowres_expanded
    precip_highres = noise_field * norm_k0

    # Apply thresholding if specified
    if threshold is not None:
        precip_highres[precip_highres < threshold] = 0

    # Return the downscaled field and optionally the spectral slope alpha
    if return_alpha:
        return precip_highres, alpha

    return precip_highres
Downscale a rainfall field by increasing its spatial resolution by a positive
integer factor.

Parameters
----------
precip: array_like
Array of shape (m, n) containing the input field.
The input is expected to contain rain rate values.
All values are required to be finite.
alpha: float, optional
Spectral slope. If None, the slope is estimated from
the input array.
ds_factor: positive int
Downscaling factor, it specifies by how many times
to increase the initial grid resolution.
threshold: float, optional
Set all values lower than the threshold to zero.
return_alpha: bool, optional
Whether to return the estimated spectral slope ``alpha``.
kernel_type: {None, "gaussian", "uniform", "tophat"}
The name of the smoothing operator. If None no smoothing is applied.
spectral_fusion: bool, optional
Whether to apply spectral merging as in :cite:`DOnofrio2014`.

Returns
-------
precip_highres: ndarray
Array of shape (m * ds_factor, n * ds_factor) containing
the downscaled field.
alpha: float
Returned only when ``return_alpha=True``.

Notes
-----
Currently, the pysteps implementation of RainFARM only covers spatial downscaling.
That is, it can improve the spatial resolution of a rainfall field. However, unlike
the original algorithm from Rebora et al. (2006), it cannot downscale the temporal
dimension. It implements spectral merging from D'Onofrio et al. (2014).

References
----------
:cite:`Rebora2006`
:cite:`DOnofrio2014`
pysteps/downscaling/rainfarm.py