def initialize_nonparam_2d_nested_filter(field, gridres=1.0, **kwargs):
    """
    Function to compute the local Fourier filters using a nested approach.

    Parameters
    ----------
    field: array-like
        Two- or three-dimensional array containing one or more input fields.
        All values are required to be finite.
        If more than one field are passed, the average fourier filter is returned.
        It assumes that fields are stacked by the first axis: [nr_fields, y, x].
    gridres: float
        Grid resolution in km.

    Other Parameters
    ----------------
    max_level: int
        Localization parameter. 0: global noise, >0: increasing degree of
        localization (default 3).
    win_fun: {'hann', 'tukey', None}
        Optional tapering function to be applied to the input field, generated with
        :py:func:`pysteps.utils.tapering.compute_window_function`
        (default 'tukey').
    war_thr: float [0;1]
        Threshold for the minimum fraction of rain needed for computing the FFT
        (default 0.1).
    rm_rdisc: bool
        Whether or not to remove the rain/no-rain disconituity. It assumes no-rain
        pixels are assigned with lowest value.
    fft_method: str or tuple
        A string or a (function,kwargs) tuple defining the FFT method to use
        (see "FFT methods" in :py:func:`pysteps.utils.interface.get_method`).
        Defaults to "numpy".

    Returns
    -------
    field: array-like
        Four-dimensional array containing the 2d fourier filters distributed over
        a 2d spatial grid.
        It can be passed to
        :py:func:`pysteps.noise.fftgenerators.generate_noise_2d_ssft_filter`.
    """

    if len(field.shape) < 2 or len(field.shape) > 3:
        raise ValueError("the input is not two- or three-dimensional array")
    if np.any(np.isnan(field)):
        raise ValueError("field must not contain NaNs")

    # defaults
    max_level = kwargs.get("max_level", 3)
    win_fun = kwargs.get("win_fun", "tukey")
    war_thr = kwargs.get("war_thr", 0.1)
    rm_rdisc = kwargs.get("rm_disc", True)
    fft = kwargs.get("fft_method", "numpy")
    if type(fft) == str:
        fft_shape = field.shape if len(field.shape) == 2 else field.shape[1:]
        fft = utils.get_method(fft, shape=fft_shape)

    field = field.copy()

    # remove rain/no-rain discontinuity
    if rm_rdisc:
        field[field > field.min()] -= field[field > field.min()].min() - field.min()

    # dims
    if len(field.shape) == 2:
        field = field[None, :, :]
    nr_fields = field.shape[0]
    dim = field.shape[1:]
    dim_x = dim[1]
    dim_y = dim[0]

    # make sure non-rainy pixels are set to zero
    field -= field.min(axis=(1, 2))[:, None, None]

    # Nested algorithm

    # prepare indices
    Idxi = np.array([[0, dim_y]])
    Idxj = np.array([[0, dim_x]])
    Idxipsd = np.array([[0, 2**max_level]])
    Idxjpsd = np.array([[0, 2**max_level]])

    # generate the FFT sample frequencies
    freqx = fft.fftfreq(dim_x, gridres)
    freqy = fft.fftfreq(dim_y, gridres)
    fx, fy = np.meshgrid(freqx, freqy)
    freq_grid = np.sqrt(fx**2 + fy**2)

    # domain fourier filter
    F0 = initialize_nonparam_2d_fft_filter(
        field, win_fun=win_fun, donorm=True, use_full_fft=True, fft_method=fft
    )["field"]
    # and allocate it to the final grid
    F = np.zeros((2**max_level, 2**max_level, F0.shape[0], F0.shape[1]))
    F += F0[np.newaxis, np.newaxis, :, :]

    # now loop levels and build composite spectra
    level = 0
    while level < max_level:
        for m in range(len(Idxi)):
            # the indices of rainfall field
            Idxinext, Idxjnext = _split_field(Idxi[m, :], Idxj[m, :], 2)
            # the indices of the field of fourier filters
            Idxipsdnext, Idxjpsdnext = _split_field(Idxipsd[m, :], Idxjpsd[m, :], 2)

            for n in range(len(Idxinext)):
                mask = _get_mask(dim, Idxinext[n, :], Idxjnext[n, :], win_fun)
                war = np.sum((field * mask[None, :, :]) > 0.01) / float(
                    (Idxinext[n, 1] - Idxinext[n, 0])
                    * (Idxjnext[n, 1] - Idxjnext[n, 0])
                    * nr_fields
                )

                if war > war_thr:
                    # the new filter
                    newfilter = initialize_nonparam_2d_fft_filter(
                        field * mask[None, :, :],
                        win_fun=None,
                        donorm=True,
                        use_full_fft=True,
                        fft_method=fft,
                    )["field"]

                    # compute logistic function to define weights as function of frequency
                    # k controls the shape of the weighting function
                    # TODO: optimize parameters
                    k = 0.05
                    x0 = (
                        Idxinext[n, 1] - Idxinext[n, 0]
                    ) / 2.0  # TODO: consider y dimension, too
                    merge_weights = 1 / (
                        1 + np.exp(-k * (1 / freq_grid - x0 * gridres))
                    )
                    newfilter *= 1 - merge_weights

                    # perform the weighted average of previous and new fourier filters
                    F[
                        Idxipsdnext[n, 0] : Idxipsdnext[n, 1],
                        Idxjpsdnext[n, 0] : Idxjpsdnext[n, 1],
                        :,
                        :,
                    ] *= merge_weights[np.newaxis, np.newaxis, :, :]
                    F[
                        Idxipsdnext[n, 0] : Idxipsdnext[n, 1],
                        Idxjpsdnext[n, 0] : Idxjpsdnext[n, 1],
                        :,
                        :,
                    ] += newfilter[np.newaxis, np.newaxis, :, :]

        # update indices
        level += 1
        Idxi, Idxj = _split_field((0, dim[0]), (0, dim[1]), 2**level)
        Idxipsd, Idxjpsd = _split_field((0, 2**max_level), (0, 2**max_level), 2**level)

    return {"field": F, "input_shape": field.shape[1:], "use_full_fft": True}
Function to compute the local Fourier filters using a nested approach.

Parameters
----------
field: array-like
Two- or three-dimensional array containing one or more input fields.
All values are required to be finite.
If more than one field are passed, the average fourier filter is returned.
It assumes that fields are stacked by the first axis: [nr_fields, y, x].
gridres: float
Grid resolution in km.

Other Parameters
----------------
max_level: int
Localization parameter. 0: global noise, >0: increasing degree of
localization (default 3).
win_fun: {'hann', 'tukey', None}
Optional tapering function to be applied to the input field, generated with
:py:func:`pysteps.utils.tapering.compute_window_function`
(default 'tukey').
war_thr: float [0;1]
Threshold for the minimum fraction of rain needed for computing the FFT
(default 0.1).
rm_rdisc: bool
Whether or not to remove the rain/no-rain disconituity. It assumes no-rain
pixels are assigned with lowest value.
fft_method: str or tuple
A string or a (function,kwargs) tuple defining the FFT method to use
(see "FFT methods" in :py:func:`pysteps.utils.interface.get_method`).
Defaults to "numpy".

Returns
-------
field: array-like
Four-dimensional array containing the 2d fourier filters distributed over
a 2d spatial grid.
It can be passed to
:py:func:`pysteps.noise.fftgenerators.generate_noise_2d_ssft_filter`.
pysteps/noise/fftgenerators.py