def _apply_spectral_fusion(
    array_low, array_high, freq_array_low, freq_array_high, ds_factor
):
    """
    Apply spectral fusion to merge two arrays in the frequency domain.
    """

    # Validate inputs
    if array_low.shape != freq_array_low.shape:
        raise ValueError("Shape of array_low must match shape of freq_array_low.")
    if array_high.shape != freq_array_high.shape:
        raise ValueError("Shape of array_high must match shape of freq_array_high.")

    nax, _ = np.shape(array_low)
    nx, _ = np.shape(array_high)
    k0 = nax // 2

    # Calculate power spectral density at specific frequency
    def compute_psd(array, fft_size):
        return rapsd(array, fft_method=np.fft)[k0 - 1] * fft_size**2

    psd_low = compute_psd(array_low, nax)
    psd_high = compute_psd(array_high, nx)

    # Normalize high-resolution array
    normalization_factor = np.sqrt(psd_low / psd_high)
    array_high *= normalization_factor

    # Perform FFT on both arrays
    fft_low = np.fft.fft2(array_low)
    fft_high = np.fft.fft2(array_high)

    # Initialize the merged FFT array with low-resolution data
    fft_merged = np.zeros_like(fft_high, dtype=np.complex128)
    fft_merged[0:k0, 0:k0] = fft_low[0:k0, 0:k0]
    fft_merged[nx - k0 : nx, 0:k0] = fft_low[k0 : 2 * k0, 0:k0]
    fft_merged[0:k0, nx - k0 : nx] = fft_low[0:k0, k0 : 2 * k0]
    fft_merged[nx - k0 : nx, nx - k0 : nx] = fft_low[k0 : 2 * k0, k0 : 2 * k0]

    fft_merged[k0, 0] = np.conj(fft_merged[nx - k0, 0])
    fft_merged[0, k0] = np.conj(fft_merged[0, nx - k0])

    # Compute frequency arrays
    freq_i = np.fft.fftfreq(nx, d=1 / ds_factor)
    freq_i = np.tile(freq_i, (nx, 1))
    freq_j = freq_i.T

    # Compute frequency domain adjustment
    ddx = np.pi * (1 / nax - 1 / nx) / np.abs(freq_i[0, 1] - freq_i[0, 0])
    freq_squared_high = freq_array_high**2
    freq_squared_low_center = freq_array_low[k0, k0] ** 2

    # Fuse in the frequency domain
    mask_high = freq_squared_high > freq_squared_low_center
    mask_low = ~mask_high
    fft_merged = fft_high * mask_high + fft_merged * mask_low * np.exp(
        -1j * ddx * freq_i - 1j * ddx * freq_j
    )

    # Inverse FFT to obtain the merged array in the spatial domain
    merged = np.real(np.fft.ifftn(fft_merged)) / fft_merged.size

    return merged
Apply spectral fusion to merge two arrays in the frequency domain.
pysteps/downscaling/rainfarm.py