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