def csd(timeseries, other, segmentlength, noverlap=None, **kwargs):
    """Calculate the CSD of two `TimeSeries` using Welch's method

    Parameters
    ----------
    timeseries : `~gwpy.timeseries.TimeSeries`
        time-series of data

    other : `~gwpy.timeseries.TimeSeries`
        time-series of data

    segmentlength : `int`
        number of samples in single average.

    noverlap : `int`
        number of samples to overlap between segments, defaults to 50%.

    **kwargs
        other keyword arguments are passed to :meth:`scipy.signal.csd`

    Returns
    -------
    spectrum : `~gwpy.frequencyseries.FrequencySeries`
        average power `FrequencySeries`

    See also
    --------
    scipy.signal.csd
    """
    # calculate CSD
    try:
        freqs, csd_ = scipy.signal.csd(
            timeseries.value, other.value, noverlap=noverlap,
            fs=timeseries.sample_rate.decompose().value,
            nperseg=segmentlength, **kwargs)
    except AttributeError as exc:
        exc.args = ('{}, scipy>=0.16 is required'.format(str(exc)),)
        raise

    # generate FrequencySeries and return
    unit = scale_timeseries_unit(timeseries.unit,
                                 kwargs.get('scaling', 'density'))
    return FrequencySeries(
        csd_, unit=unit, frequencies=freqs,
        name=str(timeseries.name)+'---'+str(other.name),
        epoch=timeseries.epoch, channel=timeseries.channel)