docs: def compute_broadening(wv_obj, sp_obj, wv_ref, sp_ref,
                         sigmalist, fpixminmax=None, naround_zero=None,
                         ax1=None, ax2=None,
                         debugplot=0):
      """Compute broadening to match 'sp_obj' with 'sp_ref'.
  
      The reference spectrum 'sp_ref' should have a better spectral
      resolution than the object spectrum 'sp_obj'.
  
      It is assumed that the wavelength arrays are provided in
      ascending order.
  
      The wavelength sampling 'wv_ref' and 'wv_ref' can be different.
  
      The comparison is performed in the common wavelength
      interval. Within this interval, both spectra are linearly
      resampled to the result of merging 'wv_obj' and 'wv_ref'.
  
      Parameters
      ----------
      wv_obj : numpy array
          Wavelength sampling of the object spectrum.
      sp_obj : numpy array
          Flux of the object spectrum.
      wv_ref : numpy array
          Wavelength sampling of the reference spectrum.
      sp_ref : numpy array
          Flux of the reference spectrum.
      sigmalist : numpy array or list
          Sigma broadening (in pixels).
      fpixminmax : tuple of floats or None
          Number of subintervals in which the pixel range (in the
          wavelength direction) is divided in order to determine
          de frequency filtering. If None, no frequency filtering is
          employed. For example, a tuple (8, 1000) means that the
          lower frequency cut is set to 8/npix, whereas the upper
          frequency cut is set to 1000/npix, being npix the number of
          pixels in the wavelength direction.
      naround_zero : int
          Half width of the window (around zero offset) to look for
          the correlation peak. If None, the whole correlation
          spectrum is employed. Otherwise, the peak will be sought
          in the interval [-naround_zero, +naround_zero].
      ax1 : matplotlib Axes object
          If not none, this plot represents the offset and fpeak
          variation as a function of sigma.
      ax2 : matplotlib Axes object
          If not none, this plot represents the two input spectra and
          the broadened reference spectrum.
      debugplot : int
          Determines whether intermediate computations and/or plots
          are displayed. The valid codes are defined in
          numina.array.display.pause_debugplot.
  
      Results
      -------
      offset_broad : float
          Offset (in pixels) between the two input spectra measured
          at the sigma value where fpeak is maximum.
      sigma_broad : float
          Sigma (in pixels) for which fpeak is maximum.
      sp_ref_broad : numpy array
          Reference spectrum broadened with sigma_broad (at the original
          wavelength sampling).
  
      """
  
      # determine wether there is an intersection between the two
      # wavelength ranges
      minmax_obj = [min(wv_obj), max(wv_obj)]
      minmax_ref = [min(wv_ref), max(wv_ref)]
      minimum_max = min(minmax_obj[1], minmax_ref[1])
      maximum_min = max(minmax_obj[0], minmax_ref[0])
      overlap = minimum_max - maximum_min
      if overlap <= 0:
          raise ValueError("Wavelength ranges of the two spectra don't overlap!")
  
      # merge wavelength ranges (removing duplicates)
      wv = np.unique(np.sort(np.concatenate((wv_obj, wv_ref))))
      # truncate to common interval
      wv = wv[(wv >= maximum_min) * (wv <= minimum_max)]
  
      # linear interpolation of input spectrum using the merged
      # wavelength sampling
      funinterp_obj = interp1d(wv_obj, sp_obj, kind='linear')
      funinterp_ref = interp1d(wv_ref, sp_ref, kind='linear')
      flux_obj = funinterp_obj(wv)
      flux_ref = funinterp_ref(wv)
  
      # normalize each spectrum dividing by its median
      flux_obj /= np.median(flux_obj)
      flux_ref /= np.median(flux_ref)
  
      # plot initial resampled spectra
      if abs(debugplot) in (21, 22):
          ax = ximplotxy(wv, flux_ref,
                         xlabel='wavelength (Angstrom)',
                         ylabel='flux (arbitrary units)',
                         label='flux_ref', show=False)
          ax.plot(wv, flux_obj, label='flux_obj')
          ax.legend()
          pause_debugplot(debugplot=debugplot, pltshow=True)
  
      # set parameters for periodic_corr1d
      if fpixminmax is None:
          fminmax = None
      else:
          fminmax = (fpixminmax[0]/wv.size, fpixminmax[1]/wv.size)
      if naround_zero is None:
          naround_zero = int(wv.size/20)
  
      nsigmas = len(list(sigmalist))
      offset = np.zeros(nsigmas)
      fpeak = np.zeros(nsigmas)
      for i, sigma in enumerate(sigmalist):
          # broaden reference spectrum
          flux_ref_broad = gaussian_filter(flux_ref, sigma)
          # plot the two spectra
          if abs(debugplot) in (21, 22):
              ax = ximplotxy(wv, flux_ref_broad,
                             xlabel='wavelength (Angstrom)',
                             ylabel='flux (arbitrary units)',
                             label='flux_ref', show=False)
              ax.plot(wv, flux_obj, label='flux_obj')
              ax.set_title('sigma: ' + str(sigma) + ' pixels')
              ax.legend()
              pause_debugplot(debugplot=debugplot, pltshow=True)
          # periodic correlation between the two spectra
          offset[i], fpeak[i] = periodic_corr1d(
              flux_ref_broad, flux_obj,
              fminmax=fminmax,
              naround_zero=naround_zero,
              norm_spectra=True,
              debugplot=debugplot
          )
  
      if ax1 is not None:
          ax1.plot(sigmalist, offset,
                   color='C0', marker='o', linestyle='', label='offset')
          ax1.set_xlabel('sigma (pixels)')
          ax1.set_ylabel('offset (pixels)', color='C0')
          ax1_ = ax1.twinx()
          ax1_.plot(sigmalist, fpeak,
                    color='C1', marker='o', linestyle='', label='fpeak')
          ax1_.set_ylabel('fpeak', color='C1')
  
      offset_broad = offset[np.argmax(fpeak)]
      sigma_broad = sigmalist[np.argmax(fpeak)]
      sp_ref_broad = gaussian_filter(sp_ref, sigma_broad)
  
      if ax2 is not None:
          ax2.plot(wv_obj, sp_obj, label='sp_obj')
          ax2.plot(wv_ref, sp_ref, color='#aaaaaa', label='sp_ref')
          ax2.plot(wv_ref, sp_ref_broad, label='sp_ref_broad')
          ax2.set_xlabel('wavelength (Angstrom)')
          ax2.set_ylabel('flux (arbitrary units)')
          ax2.legend()
  
      return offset_broad, sigma_broad, sp_ref_broad