docs: def periodic_corr1d(sp_reference, sp_offset,
                      fminmax=None,
                      naround_zero=None,
                      norm_spectra=False,
                      plottitle=None,
                      pdf=None,
                      debugplot=0):
      """Periodic correlation between two spectra, implemented using FFT.
  
      Parameters
      ----------
      sp_reference : numpy array
          Reference spectrum.
      sp_offset : numpy array
          Spectrum which offset is going to be measured relative to the
          reference spectrum.
      fminmax : tuple of floats or None
          Minimum and maximum frequencies to be used. If None, no
          frequency filtering is employed.
      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].
      norm_spectra : bool
          If True, the filtered spectra are normalized before computing
          the correlation function. This can be important when comparing
          the peak value of this function using different spectra.
      plottitle : str
          Optional plot title.
      pdf : PdfFile object or None
          If not None, output is sent to PDF file.
      debugplot : int
          Determines whether intermediate computations and/or plots
          are displayed. The valid codes are defined in
          numina.array.display.pause_debugplot.
  
      Returns
      -------
      offset : float
          Offset between the two input spectra.
      fpeak : float
          Maximum of the cross-correlation function.
  
      """
  
      # protections
      if sp_reference.ndim != 1 or sp_offset.ndim != 1:
          raise ValueError("Invalid array dimensions")
      if sp_reference.shape != sp_offset.shape:
          raise ValueError("x and y shapes are different")
  
      if plottitle is None:
          plottitle = ' '
  
      naxis1 = len(sp_reference)
  
      xcorr = np.arange(naxis1)
      naxis1_half = int(naxis1 / 2)
      for i in range(naxis1_half):
          xcorr[i + naxis1_half] -= naxis1
      isort = xcorr.argsort()
      xcorr = xcorr[isort]
  
      if fminmax is not None:
          fmin, fmax = fminmax
          sp_reference_filtmask = filtmask(sp_reference, fmin=fmin, fmax=fmax,
                                           debugplot=debugplot)
          sp_offset_filtmask = filtmask(sp_offset, fmin=fmin, fmax=fmax,
                                        debugplot=debugplot)
          if abs(debugplot) in (21, 22):
              from numina.array.display.matplotlib_qt import plt
              xdum = np.arange(naxis1) + 1
              # reference spectrum
              ax = ximplotxy(xdum, sp_reference, show=False,
                             title='reference spectrum',
                             label='original spectrum')
              ax.plot(xdum, sp_reference_filtmask,
                      label='filtered and masked spectrum')
              ax.legend()
              plt.show()
              # offset spectrum
              ax = ximplotxy(xdum, sp_offset, show=False,
                             title='offset spectrum',
                             label='original spectrum')
              ax.plot(xdum, sp_offset_filtmask,
                      label='filtered and masked spectrum')
              ax.legend()
              plt.show()
      else:
          sp_reference_filtmask = sp_reference
          sp_offset_filtmask = sp_offset
  
      if (abs(debugplot) in (21, 22)) or (pdf is not None):
          xdum = np.arange(naxis1) + 1
          ax = ximplotxy(xdum, sp_reference_filtmask, show=False,
                         title=plottitle,
                         label='reference spectrum')
          ax.plot(xdum, sp_offset_filtmask, label='offset spectrum')
          ax.legend()
          if pdf is not None:
              pdf.savefig()
          else:
              pause_debugplot(debugplot=debugplot, pltshow=True)
  
      # normalize spectra if required
      if norm_spectra:
          sp_reference_norm = np.copy(sp_reference_filtmask)
          sp_offset_norm = np.copy(sp_offset_filtmask)
          sp_dum = np.concatenate((sp_reference_norm, sp_offset_norm))
          spmin = min(sp_dum)
          spmax = max(sp_dum)
          idum = np.where(sp_reference_norm > 0)
          sp_reference_norm[idum] /= spmax
          idum = np.where(sp_reference_norm < 0)
          sp_reference_norm[idum] /= -spmin
          idum = np.where(sp_offset_norm > 0)
          sp_offset_norm[idum] /= spmax
          idum = np.where(sp_offset_norm < 0)
          sp_offset_norm[idum] /= -spmin
          if (abs(debugplot) in (21, 22)) or (pdf is not None):
              xdum = np.arange(naxis1) + 1
              ax = ximplotxy(xdum, sp_reference_norm, show=False,
                             title=plottitle + ' [normalized]',
                             label='reference spectrum')
              ax.plot(xdum, sp_offset_norm, label='offset spectrum')
              ax.legend()
              if pdf is not None:
                  pdf.savefig()
              else:
                  pause_debugplot(debugplot=debugplot, pltshow=True)
      else:
          sp_reference_norm = sp_reference_filtmask
          sp_offset_norm = sp_offset_filtmask
  
      corr = np.fft.ifft(np.fft.fft(sp_offset_norm) *
                         np.fft.fft(sp_reference_norm).conj()).real
      corr = corr[isort]
  
      # determine correlation peak
      if naround_zero is None:
          iminpeak = 0
          imaxpeak = naxis1
      else:
          iminpeak = max(int(naxis1 / 2 - naround_zero), 0)
          imaxpeak = min(int(naxis1 / 2 + naround_zero), naxis1)
      ixpeak = corr[iminpeak:imaxpeak].argmax() + iminpeak
  
      # fit correlation peak with 2nd order polynomial
      nfit = 7
      nmed = nfit // 2
      imin = ixpeak - nmed
      imax = ixpeak + nmed
      lpeak_ok = True
      if imin < 0 or imax > len(corr):
          x_refined_peak = 0
          y_refined_peak = 0
          lpeak_ok = False
          poly_peak = Polynomial([0.0])
      else:
          x_fit = np.arange(-nmed, nmed + 1, dtype=np.float)
          y_fit = corr[imin:(imax+1)]
          poly_peak = Polynomial.fit(x_fit, y_fit, 2)
          poly_peak = Polynomial.cast(poly_peak)
          coef = poly_peak.coef
          if coef[2] != 0:
              x_refined_peak = -coef[1] / (2.0 * coef[2])
          else:
              x_refined_peak = 0.0
          y_refined_peak = poly_peak(x_refined_peak)
          x_refined_peak += ixpeak
  
      offset = x_refined_peak - naxis1_half
      fpeak = y_refined_peak
  
      if (abs(debugplot) % 10 != 0) or (pdf is not None):
          ax = ximplotxy(xcorr, corr,
                         xlabel='offset (pixels)',
                         ylabel='cross-correlation function',
                         title=plottitle,
                         xlim=(-naxis1/2, naxis1/2), show=False)
          ax.axvline(offset, color='grey', linestyle='dashed')
          coffset = "(offset:{0:6.2f} pixels)".format(offset)
          ax.text(0.01, 0.99, coffset,
                  horizontalalignment='left',
                  verticalalignment='top',
                  transform=ax.transAxes)
          if naround_zero is not None:
              cwindow = "(peak region: [{},{}] pixels)".format(-naround_zero,
                                                               naround_zero)
              ax.text(0.01, 0.93, cwindow,
                      horizontalalignment='left',
                      verticalalignment='top',
                      transform=ax.transAxes)
          # inset plot
          inset_ax = inset_axes(
              ax,
              width="40%",
              height="40%",
              loc=1
          )
          inset_ax.plot(xcorr, corr)
          if naround_zero is not None:
              inset_ax.set_xlim([-naround_zero, naround_zero])
          else:
              inset_ax.set_xlim([-50, 50])
          if lpeak_ok:
              xplot = np.arange(-nmed, nmed, 0.5)
              yplot = poly_peak(xplot)
              xplot += ixpeak - naxis1_half
              inset_ax.plot(xplot, yplot, '-')
              inset_ax.plot([x_refined_peak - naxis1_half],
                            [y_refined_peak], 'o')
          inset_ax.axvline(offset, color='grey', linestyle='dashed')
          if pdf is not None:
              pdf.savefig()
          else:
              pause_debugplot(debugplot=debugplot,
                              tight_layout=False, pltshow=True)
  
      return offset, fpeak