def detect_loud_glitches(strain, psd_duration=4., psd_stride=2.,
                         psd_avg_method='median', low_freq_cutoff=30.,
                         threshold=50., cluster_window=5., corrupt_time=4.,
                         high_freq_cutoff=None, output_intermediates=False):
    """Automatic identification of loud transients for gating purposes.

    This function first estimates the PSD of the input time series using the
    FindChirp Welch method. Then it whitens the time series using that
    estimate. Finally, it computes the magnitude of the whitened series,
    thresholds it and applies the FindChirp clustering over time to the
    surviving samples.

    Parameters
    ----------
    strain : TimeSeries
        Input strain time series to detect glitches over.
    psd_duration : {float, 4}
        Duration of the segments for PSD estimation in seconds.
    psd_stride : {float, 2}
        Separation between PSD estimation segments in seconds.
    psd_avg_method : {string, 'median'}
        Method for averaging PSD estimation segments.
    low_freq_cutoff : {float, 30}
        Minimum frequency to include in the whitened strain.
    threshold : {float, 50}
        Minimum magnitude of whitened strain for considering a transient to
        be present.
    cluster_window : {float, 5}
        Length of time window to cluster surviving samples over, in seconds.
    corrupt_time : {float, 4}
        Amount of time to be discarded at the beginning and end of the input
        time series.
    high_frequency_cutoff : {float, None}
        Maximum frequency to include in the whitened strain. If given, the
        input series is downsampled accordingly. If omitted, the Nyquist
        frequency is used.
    output_intermediates : {bool, False}
        Save intermediate time series for debugging.
    """
    logging.info('Autogating: tapering strain')
    taper_length = int(corrupt_time * strain.sample_rate)
    w = numpy.arange(taper_length) / float(taper_length)
    strain[0:taper_length] *= pycbc.types.Array(w, dtype=strain.dtype)
    strain[(len(strain)-taper_length):] *= pycbc.types.Array(w[::-1],
                                                             dtype=strain.dtype)

    # don't waste time trying to optimize a single FFT
    pycbc.fft.fftw.set_measure_level(0)

    if high_freq_cutoff:
        logging.info('Autogating: downsampling strain')
        strain = resample_to_delta_t(strain, 0.5 / high_freq_cutoff,
                                     method='ldas')
    if output_intermediates:
        strain.save_to_wav('strain_conditioned.wav')

    corrupt_length = int(corrupt_time * strain.sample_rate)

    # zero-pad strain to a power-of-2 length
    strain_pad_length = next_power_of_2(len(strain))
    pad_start = strain_pad_length/2 - len(strain)/2
    pad_end = pad_start + len(strain)
    strain_pad = pycbc.types.TimeSeries(
            pycbc.types.zeros(strain_pad_length, dtype=strain.dtype),
            delta_t=strain.delta_t, copy=False,
            epoch=strain.start_time-pad_start/strain.sample_rate)
    strain_pad[pad_start:pad_end] = strain[:]

    logging.info('Autogating: estimating PSD')
    psd = pycbc.psd.welch(strain[corrupt_length:(len(strain)-corrupt_length)],
                          seg_len=int(psd_duration * strain.sample_rate),
                          seg_stride=int(psd_stride * strain.sample_rate),
                          avg_method=psd_avg_method,
                          require_exact_data_fit=False)
    psd = pycbc.psd.interpolate(psd, 1./strain_pad.duration)
    psd = pycbc.psd.inverse_spectrum_truncation(
            psd, int(psd_duration * strain.sample_rate),
            low_frequency_cutoff=low_freq_cutoff,
            trunc_method='hann')
    kmin = int(low_freq_cutoff / psd.delta_f)
    psd[0:kmin] = numpy.inf
    if high_freq_cutoff:
        kmax = int(high_freq_cutoff / psd.delta_f)
        psd[kmax:] = numpy.inf

    logging.info('Autogating: time -> frequency')
    strain_tilde = pycbc.types.FrequencySeries(
            pycbc.types.zeros(len(strain_pad) / 2 + 1,
                              dtype=pycbc.types.complex_same_precision_as(strain)),
            delta_f=psd.delta_f, copy=False)
    pycbc.fft.fft(strain_pad, strain_tilde)

    logging.info('Autogating: whitening')
    if high_freq_cutoff:
        norm = high_freq_cutoff - low_freq_cutoff
    else:
        norm = strain.sample_rate/2. - low_freq_cutoff
    strain_tilde *= (psd * norm) ** (-0.5)

    logging.info('Autogating: frequency -> time')
    pycbc.fft.ifft(strain_tilde, strain_pad)

    pycbc.fft.fftw.set_measure_level(pycbc.fft.fftw._default_measurelvl)

    if output_intermediates:
        strain_pad[pad_start:pad_end].save_to_wav('strain_whitened.wav')

    logging.info('Autogating: computing magnitude')
    mag = abs(strain_pad[pad_start:pad_end])
    if output_intermediates:
        mag.save('strain_whitened_mag.npy')

    mag = mag.numpy()

    # remove strain corrupted by filters at the ends
    mag[0:corrupt_length] = 0
    mag[-1:-corrupt_length-1:-1] = 0

    logging.info('Autogating: finding loud peaks')
    indices = numpy.where(mag > threshold)[0]
    cluster_idx = pycbc.events.findchirp_cluster_over_window(
            indices, numpy.array(mag[indices]),
            int(cluster_window*strain.sample_rate))
    times = [idx * strain.delta_t + strain.start_time \
             for idx in indices[cluster_idx]]
    return times