def bias_term(ci, ref, meta_dict, n_range):
    """
    Compute the bias term to be subtracted off the cross spectrum to compute
    the covariance spectrum. Equation in Equation in footnote 4 (section 2.1.3,
    page 12) of Uttley et al. 2014.

    Assumes power spectra are raw (not at all normalized, and not Poisson-noise-
    subtracted).

    Parameters
    ----------
    ci : Band object
        The channel of interest or interest band. ci.power is raw (not
        normalized and not Poisson-noise-subtracted), of the frequencies to be
        averaged over, with size = 1 (if avg over freq) or n_bins/2+1 (if avg
        over energy). ci.mean_rate is in units of cts/s (size = 1 if avg over
        energy, size = detchans if avg over freq). Power and frequency
        have only the positive (tot en met Nyquist) frequencies.

    ref : Band object
        The reference band. ref.power is raw (not normalized and not Poisson-
        noise-subtracted), of the frequencies to be averaged over, with
        size = 1 (if avg over freq) or n_bins/2+1 (if avg over energy).
        ref.mean_rate is in units of cts/s. Power and frequency have only the
        positive (tot en met Nyquist) frequencies.

    meta_dict : dict
        Dictionary of meta-parameters needed for analysis.

    n_range : int
        Number of bins that will be averaged together for the lags. Energy bins
        for frequency lags, frequency bins for energy lags.

    Returns
    -------
    n_squared : float
        The bias term to be subtracted off the cross spectrum for computing the
        covariance spectrum. Equation in footnote 4 (section 2.1.3, page 12) of
        Uttley et al. 2014. If you get an undefined error for the bias term,
        just set it equal to zero.

    """
    ## Compute the Poisson noise level in absolute rms units
    Pnoise_ref = ref.mean_rate * 2.0
    Pnoise_ci = ci.mean_rate * 2.0

    ## Normalizing power spectra to absolute rms normalization
    ## Not subtracting the noise (yet)!
    abs_ci = ci.power * (2.0 * meta_dict['dt'] / float(n_range))
    abs_ref = ref.power * (2.0 * meta_dict['dt'] / float(n_range))

    temp_a = (abs_ref - Pnoise_ref) * Pnoise_ci
    temp_b = (abs_ci - Pnoise_ci) * Pnoise_ref
    temp_c = Pnoise_ref * Pnoise_ci

    n_squared = (temp_a + temp_b + temp_c) / (n_range * meta_dict['n_seg'])
    return n_squared