def from_string(psd_name, length, delta_f, low_freq_cutoff):
    """Generate a frequency series containing a LALSimulation PSD specified by name.

    Parameters
    ----------
    psd_name : string
        PSD name as found in LALSimulation, minus the SimNoisePSD prefix.
    length : int
        Length of the frequency series in samples.
    delta_f : float
        Frequency resolution of the frequency series.
    low_freq_cutoff : float
        Frequencies below this value are set to zero.

    Returns
    -------
    psd : FrequencySeries
        The generated frequency series.
    """
    if psd_name not in _psd_list:
        raise ValueError(psd_name + ' not found among LALSimulation PSD functions.')
    kmin = int(low_freq_cutoff / delta_f)
    lalseries = lal.CreateREAL8FrequencySeries(
        '', lal.LIGOTimeGPS(0), 0, delta_f, lal.DimensionlessUnit, length)
    try:
        func = lalsimulation.__dict__[_name_prefix + psd_name + _name_suffix]
    except KeyError:
        func = lalsimulation.__dict__[_name_prefix + psd_name]
        func(lalseries, low_freq_cutoff)
    else:
        lalsimulation.SimNoisePSD(lalseries, 0, func)
    psd = FrequencySeries(lalseries.data.data, delta_f=delta_f)
    psd.data[:kmin] = 0
    return psd

def td_waveform_to_fd_waveform(waveform, out=None, length=None,
                               buffer_length=100):
    """ Convert a time domain into a frequency domain waveform by FFT.
        As a waveform is assumed to "wrap" in the time domain one must be
        careful to ensure the waveform goes to 0 at both "boundaries". To
        ensure this is done correctly the waveform must have the epoch set such
        the merger time is at t=0 and the length of the waveform should be
        shorter than the desired length of the FrequencySeries (times 2 - 1)
        so that zeroes can be suitably pre- and post-pended before FFTing.
        If given, out is a memory array to be used as the output of the FFT.
        If not given memory is allocated internally.
        If present the length of the returned FrequencySeries is determined
        from the length out. If out is not given the length can be provided
        expicitly, or it will be chosen as the nearest power of 2. If choosing
        length explicitly the waveform length + buffer_length is used when
        choosing the nearest binary number so that some zero padding is always
        added.
    """
    # Figure out lengths and set out if needed
    if out is None:
        if length is None:
            N = pnutils.nearest_larger_binary_number(len(waveform) + \
                                                     buffer_length)
            n = int(N//2) + 1
        else:
            n = length
            N = (n-1)*2
        out = zeros(n, dtype=complex_same_precision_as(waveform))
    else:
        n = len(out)
        N = (n-1)*2
    delta_f =  1. / (N * waveform.delta_t)

    # total duration of the waveform
    tmplt_length = len(waveform) * waveform.delta_t
    if len(waveform) > N:
        err_msg = "The time domain template is longer than the intended "
        err_msg += "duration in the frequency domain. This situation is "
        err_msg += "not supported in this function. Please shorten the "
        err_msg += "waveform appropriately before calling this function or "
        err_msg += "increase the allowed waveform length. "
        err_msg += "Waveform length (in samples): {}".format(len(waveform))
        err_msg += ". Intended length: {}.".format(N)
        raise ValueError(err_msg)
    # for IMR templates the zero of time is at max amplitude (merger)
    # thus the start time is minus the duration of the template from
    # lower frequency cutoff to merger, i.e. minus the 'chirp time'
    tChirp = - float( waveform.start_time )  # conversion from LIGOTimeGPS
    waveform.resize(N)
    k_zero = int(waveform.start_time / waveform.delta_t)
    waveform.roll(k_zero)
    htilde = FrequencySeries(out, delta_f=delta_f, copy=False)
    fft(waveform.astype(real_same_precision_as(htilde)), htilde)
    htilde.length_in_time = tmplt_length
    htilde.chirp_length = tChirp
    return htilde

 def _shift_and_ifft(self, fdsinx, tshift, fseries=None):
     """Calls apply_fd_time_shift, and iFFTs to the time domain.
     """
     start_time = self.time_series.start_time
     tdshift = apply_fd_time_shift(fdsinx, start_time+tshift,
                                   fseries=fseries)
     if not isinstance(tdshift, FrequencySeries):
         # cast to FrequencySeries so time series will work
         tdshift = FrequencySeries(tdshift, delta_f=fdsinx.delta_f,
                                   epoch=fdsinx.epoch)
     return tdshift.to_timeseries()

def get_waveform_filter(out, template=None, **kwargs):
    """Return a frequency domain waveform filter for the specified approximant
    """
    n = len(out)

    input_params = props(template, **kwargs)

    if input_params['approximant'] in filter_approximants(_scheme.mgr.state):
        wav_gen = filter_wav[type(_scheme.mgr.state)]
        htilde = wav_gen[input_params['approximant']](out=out, **input_params)
        htilde.resize(n)
        htilde.chirp_length = get_waveform_filter_length_in_time(**input_params)
        htilde.length_in_time = htilde.chirp_length
        return htilde

    if input_params['approximant'] in fd_approximants(_scheme.mgr.state):
        wav_gen = fd_wav[type(_scheme.mgr.state)]
        hp, hc = wav_gen[input_params['approximant']](**input_params)
        hp.resize(n)
        hp.chirp_length = get_waveform_filter_length_in_time(**input_params)
        hp.length_in_time = hp.chirp_length
        return hp

    elif input_params['approximant'] in td_approximants(_scheme.mgr.state):
        # N: number of time samples required
        N = (n-1)*2
        delta_f = 1.0 / (N * input_params['delta_t'])
        wav_gen = td_wav[type(_scheme.mgr.state)]
        hp, hc = wav_gen[input_params['approximant']](**input_params)
        # taper the time series hp if required
        if ('taper' in input_params.keys() and \
            input_params['taper'] is not None):
            hp = wfutils.taper_timeseries(hp, input_params['taper'],
                                          return_lal=False)
        # total duration of the waveform
        tmplt_length = len(hp) * hp.delta_t
        # for IMR templates the zero of time is at max amplitude (merger)
        # thus the start time is minus the duration of the template from
        # lower frequency cutoff to merger, i.e. minus the 'chirp time'
        tChirp = - float( hp.start_time )  # conversion from LIGOTimeGPS
        hp.resize(N)
        k_zero = int(hp.start_time / hp.delta_t)
        hp.roll(k_zero)
        htilde = FrequencySeries(out, delta_f=delta_f, copy=False)
        fft(hp.astype(real_same_precision_as(htilde)), htilde)
        htilde.length_in_time = tmplt_length
        htilde.chirp_length = tChirp
        return htilde

    else:
        raise ValueError("Approximant %s not available" %
                            (input_params['approximant']))

def from_string(psd_name, length, delta_f, low_freq_cutoff):
    """Generate a frequency series containing a LALSimulation PSD specified
    by name.

    Parameters
    ----------
    psd_name : string
        PSD name as found in LALSimulation, minus the SimNoisePSD prefix.
    length : int
        Length of the frequency series in samples.
    delta_f : float
        Frequency resolution of the frequency series.
    low_freq_cutoff : float
        Frequencies below this value are set to zero.

    Returns
    -------
    psd : FrequencySeries
        The generated frequency series.
    """

    # check if valid PSD model
    if psd_name not in get_psd_model_list():
        raise ValueError(psd_name + ' not found among analytical '
                         'PSD functions.')

    # if PSD model is in LALSimulation
    if psd_name in get_lalsim_psd_list():
        lalseries = lal.CreateREAL8FrequencySeries(
            '', lal.LIGOTimeGPS(0), 0, delta_f, lal.DimensionlessUnit, length)
        try:
            func = lalsimulation.__dict__[
                                        _name_prefix + psd_name + _name_suffix]
        except KeyError:
            func = lalsimulation.__dict__[_name_prefix + psd_name]
            func(lalseries, low_freq_cutoff)
        else:
            lalsimulation.SimNoisePSD(lalseries, 0, func)
        psd = FrequencySeries(lalseries.data.data, delta_f=delta_f)

    # if PSD model is coded in PyCBC
    else:
        func = pycbc_analytical_psds[psd_name]
        psd = func(length, delta_f, low_freq_cutoff)

    # zero-out content below low-frequency cutoff
    kmin = int(low_freq_cutoff / delta_f)
    psd.data[:kmin] = 0

    return psd

def get_waveform_filter(out, template=None, **kwargs):
    """Return a frequency domain waveform filter for the specified approximant
    """
    n = len(out)

    input_params = props(template, **kwargs)

    if input_params['approximant'] in filter_approximants(_scheme.mgr.state):
        wav_gen = filter_wav[type(_scheme.mgr.state)]
        htilde = wav_gen[input_params['approximant']](out=out, **input_params)
        htilde.resize(n)
        htilde.chirp_length = get_waveform_filter_length_in_time(**input_params)
        htilde.length_in_time = htilde.chirp_length
        return htilde

    if input_params['approximant'] in fd_approximants(_scheme.mgr.state):
        wav_gen = fd_wav[type(_scheme.mgr.state)]
        hp, hc = wav_gen[input_params['approximant']](**input_params)
        hp.resize(n)
        out[0:len(hp)] = hp[:]
        hp = FrequencySeries(out, delta_f=hp.delta_f, copy=False)
        hp.chirp_length = get_waveform_filter_length_in_time(**input_params)
        hp.length_in_time = hp.chirp_length
        return hp

    elif input_params['approximant'] in td_approximants(_scheme.mgr.state):
        # N: number of time samples required
        N = (n-1)*2
        delta_f = 1.0 / (N * input_params['delta_t'])
        wav_gen = td_wav[type(_scheme.mgr.state)]
        hp, hc = wav_gen[input_params['approximant']](**input_params)
        # taper the time series hp if required
        if ('taper' in input_params.keys() and \
            input_params['taper'] is not None):
            hp = wfutils.taper_timeseries(hp, input_params['taper'],
                                          return_lal=False)
        return td_waveform_to_fd_waveform(hp, out=out)

    else:
        raise ValueError("Approximant %s not available" %
                            (input_params['approximant']))

def flat_unity(length, delta_f, low_freq_cutoff):
    """ Returns a FrequencySeries of ones above the low_frequency_cutoff.

    Parameters
    ----------
    length : int
        Length of output Frequencyseries.
    delta_f : float
        Frequency step for output FrequencySeries.
    low_freq_cutoff : int
        Low-frequency cutoff for output FrequencySeries.

    Returns
    -------
    FrequencySeries
        Returns a FrequencySeries containing the unity PSD model.
    """
    fseries = FrequencySeries(numpy.ones(length), delta_f=delta_f)
    kmin = int(low_freq_cutoff / fseries.delta_f)
    fseries.data[:kmin] = 0
    return fseries

def _imrphenombfreq(**p):
    import lalinspiral
    params = lalinspiral.InspiralTemplate()
    m1 = p['mass1']
    m2 = p['mass2']

    mc, et = pnutils.mass1_mass2_to_mchirp_eta(m1, m2)
    params.approximant = lalsimulation.IMRPhenomB
    params.fLower = p['f_lower']
    params.eta = et
    params.distance = p['distance'] * lal.PC_SI * 1e6
    params.mass1 = m1
    params.mass2 = m2
    params.spin1[2] = p['spin1z']
    params.spin2[2] = p['spin2z']
    params.startPhase = p['coa_phase']*2 - lal.PI
    params.startTime = 0

    params.tSampling = 8192
    N = int(params.tSampling / p['delta_f'])
    n = N / 2

    # Create temporary memory to hold the results and call the generator
    hpt = zeros(N, dtype=float32)
    hct = zeros(N, dtype=float32)
    hpt=hpt.lal()
    hct=hct.lal()
    lalinspiral.BBHPhenWaveBFreqDomTemplates(hpt, hct, params)

    # Copy the results to a complex frequencyseries format
    hctc = FrequencySeries(zeros(n, dtype=complex64), delta_f=p['delta_f'])
    hptc = FrequencySeries(zeros(n, dtype=complex64), delta_f=p['delta_f'])

    hptc.data += hpt.data[0:n]
    hptc.data[1:n] += hpt.data[N:N-n:-1] * 1j

    hctc.data += hct.data[0:n]
    hctc.data[1:n] += hct.data[N:N-n:-1] * 1j

    return hptc.astype(complex128),  hctc.astype(complex128)

def from_lalsimulation(func, length, delta_f, low_freq_cutoff):
    """Generate a frequency series containing the specified LALSimulation PSD.

    Parameters
    ----------
    func : function
        LALSimulation PSD function.
    length : int
        Length of the frequency series in samples.
    delta_f : float
        Frequency resolution of the frequency series.
    low_freq_cutoff : float
        Frequencies below this value are set to zero.

    Returns
    -------
    psd : FrequencySeries
        The generated frequency series.
    """
    psd = FrequencySeries(zeros(length), delta_f=delta_f)
    kmin = int(low_freq_cutoff / delta_f)
    psd.data[kmin:] = map(func, numpy.arange(length)[kmin:] * delta_f)
    return psd

def welch(timeseries, seg_len=4096, seg_stride=2048, window='hann',
          avg_method='median', num_segments=None, require_exact_data_fit=False):
    """PSD estimator based on Welch's method.

    Parameters
    ----------
    timeseries : TimeSeries
        Time series for which the PSD is to be estimated.
    seg_len : int
        Segment length in samples.
    seg_stride : int
        Separation between consecutive segments, in samples.
    window : {'hann'}
        Function used to window segments before Fourier transforming.
    avg_method : {'median', 'mean', 'median-mean'}
        Method used for averaging individual segment PSDs.

    Returns
    -------
    psd : FrequencySeries
        Frequency series containing the estimated PSD.

    Raises
    ------
    ValueError
        For invalid choices of `seg_len`, `seg_stride` `window` and
        `avg_method` and for inconsistent combinations of len(`timeseries`),
        `seg_len` and `seg_stride`.

    Notes
    -----
    See arXiv:gr-qc/0509116 for details.
    """
    window_map = {
        'hann': numpy.hanning
    }

    # sanity checks
    if not window in window_map:
        raise ValueError('Invalid window')
    if not avg_method in ('mean', 'median', 'median-mean'):
        raise ValueError('Invalid averaging method')
    if type(seg_len) is not int or type(seg_stride) is not int \
        or seg_len <= 0 or seg_stride <= 0:
        raise ValueError('Segment length and stride must be positive integers')

    if timeseries.precision == 'single':
        fs_dtype = numpy.complex64
    elif timeseries.precision == 'double':
        fs_dtype = numpy.complex128
        
    num_samples = len(timeseries)
    if num_segments is None:
        num_segments = int(num_samples // seg_stride)
        # NOTE: Is this not always true?
        if (num_segments - 1) * seg_stride + seg_len > num_samples:
            num_segments -= 1

    if not require_exact_data_fit:
        data_len = (num_segments - 1) * seg_stride + seg_len

        # Get the correct amount of data
        if data_len < num_samples:
            diff = num_samples - data_len
            start = diff // 2
            end = num_samples - diff // 2
            # Want this to be integers so if diff is odd, catch it here.
            if diff % 2:
                start = start + 1

            timeseries = timeseries[start:end]
            num_samples = len(timeseries)
        if data_len > num_samples:
            err_msg = "I was asked to estimate a PSD on %d " %(data_len)
            err_msg += "data samples. However the data provided only contains "
            err_msg += "%d data samples." %(num_samples)

    if num_samples != (num_segments - 1) * seg_stride + seg_len:
        raise ValueError('Incorrect choice of segmentation parameters')
        
    w = Array(window_map[window](seg_len).astype(timeseries.dtype))

    # calculate psd of each segment
    delta_f = 1. / timeseries.delta_t / seg_len
    segment_tilde = FrequencySeries(numpy.zeros(seg_len / 2 + 1), \
        delta_f=delta_f, dtype=fs_dtype)

    segment_psds = []
    for i in xrange(num_segments):
        segment_start = i * seg_stride
        segment_end = segment_start + seg_len
        segment = timeseries[segment_start:segment_end]
        assert len(segment) == seg_len
        fft(segment * w, segment_tilde)
        seg_psd = abs(segment_tilde * segment_tilde.conj()).numpy()

        #halve the DC and Nyquist components to be consistent with TO10095
        seg_psd[0] /= 2
        seg_psd[-1] /= 2

#.........这里部分代码省略.........

def inverse_spectrum_truncation(psd, max_filter_len, low_frequency_cutoff=None, trunc_method=None):
    """Modify a PSD such that the impulse response associated with its inverse
    square root is no longer than `max_filter_len` time samples. In practice
    this corresponds to a coarse graining or smoothing of the PSD.

    Parameters
    ----------
    psd : FrequencySeries
        PSD whose inverse spectrum is to be truncated.
    max_filter_len : int
        Maximum length of the time-domain filter in samples.
    low_frequency_cutoff : {None, int}
        Frequencies below `low_frequency_cutoff` are zeroed in the output.
    trunc_method : {None, 'hann'}
        Function used for truncating the time-domain filter.
        None produces a hard truncation at `max_filter_len`.

    Returns
    -------
    psd : FrequencySeries
        PSD whose inverse spectrum has been truncated.

    Raises
    ------
    ValueError
        For invalid types or values of `max_filter_len` and `low_frequency_cutoff`.

    Notes
    -----
    See arXiv:gr-qc/0509116 for details.
    """
    # sanity checks
    if type(max_filter_len) is not int or max_filter_len <= 0:
        raise ValueError('max_filter_len must be a positive integer')
    if low_frequency_cutoff is not None and low_frequency_cutoff < 0 \
        or low_frequency_cutoff > psd.sample_frequencies[-1]:
        raise ValueError('low_frequency_cutoff must be within the bandwidth of the PSD')

    N = (len(psd)-1)*2

    inv_asd = FrequencySeries((1. / psd)**0.5, delta_f=psd.delta_f, \
        dtype=complex_same_precision_as(psd))
        
    inv_asd[0] = 0
    inv_asd[N/2] = 0
    q = TimeSeries(numpy.zeros(N), delta_t=(N / psd.delta_f), \
        dtype=real_same_precision_as(psd))

    if low_frequency_cutoff:
        kmin = int(low_frequency_cutoff / psd.delta_f)
        inv_asd[0:kmin] = 0

    ifft(inv_asd, q)
    
    trunc_start = max_filter_len / 2
    trunc_end = N - max_filter_len / 2

    if trunc_method == 'hann':
        trunc_window = Array(numpy.hanning(max_filter_len), dtype=q.dtype)
        q[0:trunc_start] *= trunc_window[max_filter_len/2:max_filter_len]
        q[trunc_end:N] *= trunc_window[0:max_filter_len/2]

    q[trunc_start:trunc_end] = 0
    psd_trunc = FrequencySeries(numpy.zeros(len(psd)), delta_f=psd.delta_f, \
                                dtype=complex_same_precision_as(psd))
    fft(q, psd_trunc)
    psd_trunc *= psd_trunc.conj()
    psd_out = 1. / abs(psd_trunc)

    return psd_out

def get_fd_qnm(template=None, **kwargs):
    """Return a frequency domain damped sinusoid.

    Parameters
    ----------
    template: object
        An object that has attached properties. This can be used to substitute
        for keyword arguments. A common example would be a row in an xml table.
    f_0 : float
        The ringdown-frequency.
    tau : float
        The damping time of the sinusoid.
    amp : float
        The amplitude of the ringdown (constant for now).
    phi : float
        The initial phase of the ringdown. Should also include the information
        from the azimuthal angle (phi_0 + m*Phi).
    inclination : {None, float}, optional
        Inclination of the system in radians for the spherical harmonics.
    l : {2, int}, optional
        l mode for the spherical harmonics. Default is l=2.
    m : {2, int}, optional
        m mode for the spherical harmonics. Default is m=2.
    t_0 :  {0, float}, optional
        The starting time of the ringdown.
    delta_f : {None, float}, optional
        The frequency step used to generate the ringdown.
        If None, it will be set to the inverse of the time at which the
        amplitude is 1/1000 of the peak amplitude.
    f_lower: {None, float}, optional
        The starting frequency of the output frequency series.
        If None, it will be set to delta_f.
    f_final : {None, float}, optional
        The ending frequency of the output frequency series.
        If None, it will be set to the frequency at which the amplitude is
        1/1000 of the peak amplitude.

    Returns
    -------
    hplustilde: FrequencySeries
        The plus phase of the ringdown in frequency domain.
    hcrosstilde: FrequencySeries
        The cross phase of the ringdown in frequency domain.
    """

    input_params = props(template, qnm_required_args, **kwargs)

    f_0 = input_params.pop('f_0')
    tau = input_params.pop('tau')
    amp = input_params.pop('amp')
    phi = input_params.pop('phi')
    # the following have defaults, and so will be populated
    t_0 = input_params.pop('t_0')
    # the following may not be in input_params
    inc = input_params.pop('inclination', None)
    l = input_params.pop('l', 2)
    m = input_params.pop('m', 2)
    delta_f = input_params.pop('delta_f', None)
    f_lower = input_params.pop('f_lower', None)
    f_final = input_params.pop('f_final', None)

    if not delta_f:
        delta_f = 1. / qnm_time_decay(tau, 1./1000)
    if not f_lower:
        f_lower = delta_f
        kmin = 0
    else:
        kmin = int(f_lower / delta_f)
    if not f_final:
        f_final = qnm_freq_decay(f_0, tau, 1./1000)
    if f_final > max_freq:
        f_final = max_freq
    kmax = int(f_final / delta_f) + 1

    freqs = numpy.arange(kmin, kmax)*delta_f
    if inc is not None:
        Y_plus, Y_cross = spher_harms(l, m, inc)
    else:
        Y_plus, Y_cross = 1, 1

    denominator = 1 + (4j * pi * freqs * tau) - (4 * pi_sq * ( freqs*freqs - f_0*f_0) * tau*tau)
    norm = amp * tau / denominator
    if t_0 != 0:
        time_shift = numpy.exp(-1j * two_pi * freqs * t_0)
        norm *= time_shift

    # Analytical expression for the Fourier transform of the ringdown (damped sinusoid)
    hp_tilde = norm * Y_plus * ( (1 + 2j * pi * freqs * tau) * numpy.cos(phi)
                                        - two_pi * f_0 * tau * numpy.sin(phi) )
    hc_tilde = norm * Y_cross * ( (1 + 2j * pi * freqs * tau) * numpy.sin(phi)
                                        + two_pi * f_0 * tau * numpy.cos(phi) )

    outplus = FrequencySeries(zeros(kmax, dtype=complex128), delta_f=delta_f)
    outcross = FrequencySeries(zeros(kmax, dtype=complex128), delta_f=delta_f)
    outplus.data[kmin:kmax] = hp_tilde
    outcross.data[kmin:kmax] = hc_tilde

    return outplus, outcross

def get_two_pol_waveform_filter(outplus, outcross, template, **kwargs):
    """Return a frequency domain waveform filter for the specified approximant.
    Unlike get_waveform_filter this function returns both h_plus and h_cross
    components of the waveform, which are needed for searches where h_plus
    and h_cross are not related by a simple phase shift.
    """
    n = len(outplus)

    # If we don't have an inclination column alpha3 might be used
    if not hasattr(template, 'inclination') and 'inclination' not in kwargs:
        if hasattr(template, 'alpha3'):
            kwargs['inclination'] = template.alpha3

    input_params = props(template, **kwargs)

    if input_params['approximant'] in fd_approximants(_scheme.mgr.state):
        wav_gen = fd_wav[type(_scheme.mgr.state)]
        hp, hc = wav_gen[input_params['approximant']](**input_params)
        hp.resize(n)
        hc.resize(n)
        outplus[0:len(hp)] = hp[:]
        hp = FrequencySeries(outplus, delta_f=hp.delta_f, copy=False)
        outcross[0:len(hc)] = hc[:]
        hc = FrequencySeries(outcross, delta_f=hc.delta_f, copy=False)
        hp.chirp_length = get_waveform_filter_length_in_time(**input_params)
        hp.length_in_time = hp.chirp_length
        hc.chirp_length = hp.chirp_length
        hc.length_in_time = hp.length_in_time
        return hp, hc
    elif input_params['approximant'] in td_approximants(_scheme.mgr.state):
        # N: number of time samples required
        N = (n-1)*2
        delta_f = 1.0 / (N * input_params['delta_t'])
        wav_gen = td_wav[type(_scheme.mgr.state)]
        hp, hc = wav_gen[input_params['approximant']](**input_params)
        # taper the time series hp if required
        if 'taper' in input_params.keys() and \
                input_params['taper'] is not None:
            hp = wfutils.taper_timeseries(hp, input_params['taper'],
                                          return_lal=False)
            hc = wfutils.taper_timeseries(hc, input_params['taper'],
                                          return_lal=False)
        # total duration of the waveform
        tmplt_length = len(hp) * hp.delta_t
        # for IMR templates the zero of time is at max amplitude (merger)
        # thus the start time is minus the duration of the template from
        # lower frequency cutoff to merger, i.e. minus the 'chirp time'
        tChirp = - float( hp.start_time )  # conversion from LIGOTimeGPS
        hp.resize(N)
        hc.resize(N)
        k_zero = int(hp.start_time / hp.delta_t)
        hp.roll(k_zero)
        hc.roll(k_zero)
        hp_tilde = FrequencySeries(outplus, delta_f=delta_f, copy=False)
        hc_tilde = FrequencySeries(outcross, delta_f=delta_f, copy=False)
        fft(hp.astype(real_same_precision_as(hp_tilde)), hp_tilde)
        fft(hc.astype(real_same_precision_as(hc_tilde)), hc_tilde)
        hp_tilde.length_in_time = tmplt_length
        hp_tilde.chirp_length = tChirp
        hc_tilde.length_in_time = tmplt_length
        hc_tilde.chirp_length = tChirp
        return hp_tilde, hc_tilde
    else:
        raise ValueError("Approximant %s not available" %
                            (input_params['approximant']))

def fd_decompress(amp, phase, sample_frequencies, out=None, df=None,
        f_lower=None, interpolation='linear'):
    """Decompresses an FD waveform using the given amplitude, phase, and the
    frequencies at which they are sampled at.

    Parameters
    ----------
    amp : array
        The amplitude of the waveform at the sample frequencies.
    phase : array
        The phase of the waveform at the sample frequencies.
    sample_frequencies : array
        The frequency (in Hz) of the waveform at the sample frequencies.
    out : {None, FrequencySeries}
        The output array to save the decompressed waveform to. If this contains
        slots for frequencies > the maximum frequency in sample_frequencies,
        the rest of the values are zeroed. If not provided, must provide a df.
    df : {None, float}
        The frequency step to use for the decompressed waveform. Must be
        provided if out is None.
    f_lower : {None, float}
        The frequency to start the decompression at. If None, will use whatever
        the lowest frequency is in sample_frequencies. All values at
        frequencies less than this will be 0 in the decompressed waveform.
    interpolation : {'linear', str}
        The interpolation to use for the amplitude and phase. Default is
        'linear'. If 'linear' a custom interpolater is used. Otherwise,
        ``scipy.interpolate.interp1d`` is used; for other options, see
        possible values for that function's ``kind`` argument.

    Returns
    -------
    out : FrqeuencySeries
        If out was provided, writes to that array. Otherwise, a new
        FrequencySeries with the decompressed waveform.
    """
        
    if out is None:
        if df is None:
            raise ValueError("Either provide output memory or a df")
        flen = int(numpy.ceil(sample_frequencies.max()/df+1))
        out = FrequencySeries(numpy.zeros(flen,
            dtype=numpy.complex128), copy=False, delta_f=df)
    else:
        df = out.delta_f
        flen = len(out)
    if f_lower is None:
        jmin = 0
        f_lower = sample_frequencies[0]
    else:
        if f_lower >= sample_frequencies.max():
            raise ValueError("f_lower is > than the maximum sample frequency")
        jmin = int(numpy.searchsorted(sample_frequencies, f_lower))
    imin = int(numpy.floor(f_lower/df))
    # interpolate the amplitude and the phase
    if interpolation == "linear":
        # use custom interpolation
        sflen = len(sample_frequencies)
        h = numpy.array(out.data, copy=False)
        # make sure df is a float
        df = float(df)
        code = r"""
        # include <math.h>
        # include <stdio.h>
        int j = jmin-1;
        double sf = 0.;
        double A = 0.;
        double nextA = 0.;
        double phi = 0.;
        double nextPhi = 0.;
        double next_sf = sample_frequencies[jmin];
        double f = 0.;
        double invsdf = 0.;
        double mAmp = 0.;
        double bAmp = 0.;
        double mPhi = 0.;
        double bPhi = 0.;
        double interpAmp = 0.;
        double interpPhi = 0.;
        // zero-out beginning of array
        std::fill(h, h+imin, std::complex<double>(0., 0.));
        // cycle over desired samples
        for (int i=imin; i<flen; i++){
            f = i*df;
            if (f >= next_sf){
                // update linear interpolations
                j += 1;
                // if we have gone beyond the sampled frequencies, just break
                if ((j+1) == sflen) {
                    // zero-out rest the rest of the array & exit
                    std::fill(h+i, h+flen, std::complex<double>(0., 0.));
                    break;
                }
                sf = (double) sample_frequencies[j];
                next_sf = (double) sample_frequencies[j+1];
                A = (double) amp[j];
                nextA = (double) amp[j+1];
                phi = (double) phase[j];
                nextPhi = (double) phase[j+1];
                invsdf = 1./(next_sf - sf);
#.........这里部分代码省略.........

def fd_decompress(amp, phase, sample_frequencies, out=None, df=None,
        f_lower=None, interpolation='linear'):
    """Decompresses an FD waveform using the given amplitude, phase, and the
    frequencies at which they are sampled at.

    Parameters
    ----------
    amp : array
        The amplitude of the waveform at the sample frequencies.
    phase : array
        The phase of the waveform at the sample frequencies.
    sample_frequencies : array
        The frequency (in Hz) of the waveform at the sample frequencies.
    out : {None, FrequencySeries}
        The output array to save the decompressed waveform to. If this contains
        slots for frequencies > the maximum frequency in sample_frequencies,
        the rest of the values are zeroed. If not provided, must provide a df.
    df : {None, float}
        The frequency step to use for the decompressed waveform. Must be
        provided if out is None.
    f_lower : {None, float}
        The frequency to start the decompression at. If None, will use whatever
        the lowest frequency is in sample_frequencies. All values at
        frequencies less than this will be 0 in the decompressed waveform.
    interpolation : {'linear', str}
        The interpolation to use for the amplitude and phase. Default is
        'linear'. If 'linear' a custom interpolater is used. Otherwise,
        ``scipy.interpolate.interp1d`` is used; for other options, see
        possible values for that function's ``kind`` argument.

    Returns
    -------
    out : FrqeuencySeries
        If out was provided, writes to that array. Otherwise, a new
        FrequencySeries with the decompressed waveform.
    """
    precision = _precision_map[sample_frequencies.dtype.name]
    if _precision_map[amp.dtype.name] != precision or \
            _precision_map[phase.dtype.name] != precision:
        raise ValueError("amp, phase, and sample_points must all have the "
            "same precision")
    if out is None:
        if df is None:
            raise ValueError("Either provide output memory or a df")
        hlen = int(numpy.ceil(sample_frequencies.max()/df+1))
        out = FrequencySeries(numpy.zeros(hlen,
            dtype=_complex_dtypes[precision]), copy=False,
            delta_f=df)
    else:
        # check for precision compatibility
        if out.precision == 'double' and precision == 'single':
            raise ValueError("cannot cast single precision to double")
        df = out.delta_f
        hlen = len(out)
    if f_lower is None:
        imin = 0
        f_lower = sample_frequencies[0]
    else:
        if f_lower >= sample_frequencies.max():
            raise ValueError("f_lower is > than the maximum sample frequency")
        imin = int(numpy.searchsorted(sample_frequencies, f_lower))
    start_index = int(numpy.floor(f_lower/df))
    # interpolate the amplitude and the phase
    if interpolation == "linear":
        if precision == 'single':
            code = _linear_decompress_code32
        else:
            code = _linear_decompress_code
        # use custom interpolation
        sflen = len(sample_frequencies)
        h = numpy.array(out.data, copy=False)
        delta_f = float(df)
        inline(code, ['h', 'hlen', 'sflen', 'delta_f', 'sample_frequencies',
                      'amp', 'phase', 'start_index', 'imin'],
               extra_compile_args=[WEAVE_FLAGS + '-march=native -O3 -w'] +\
                                  omp_flags,
               libraries=omp_libs)
    else:
        # use scipy for fancier interpolation
        outfreq = out.sample_frequencies.numpy()
        amp_interp = interpolate.interp1d(sample_frequencies.numpy(),
            amp.numpy(), kind=interpolation, bounds_error=False, fill_value=0.,
            assume_sorted=True)
        phase_interp = interpolate.interp1d(sample_frequencies.numpy(),
            phase.numpy(), kind=interpolation, bounds_error=False,
            fill_value=0., assume_sorted=True)
        A = amp_interp(outfreq)
        phi = phase_interp(outfreq)
        out.data[:] = A*numpy.cos(phi) + (1j)*A*numpy.sin(phi)
    return out

def get_fd_qnm(template=None, **kwargs):
    """Return a frequency domain damped sinusoid.

    Parameters
    ----------
    template: object
        An object that has attached properties. This can be used to substitute
        for keyword arguments. A common example would be a row in an xml table.
    f_0 : float
        The ringdown-frequency.
    tau : float
        The damping time of the sinusoid.
    t_0 :  {0, float}, optional
        The starting time of the ringdown.
    phi_0 : {0, float}, optional
        The initial phase of the ringdown.
    amp : {1, float}, optional
        The amplitude of the ringdown (constant for now).
    delta_f : {None, float}, optional
        The frequency step used to generate the ringdown.
        If None, it will be set to the inverse of the time at which the
        amplitude is 1/1000 of the peak amplitude.
    f_lower: {None, float}, optional
        The starting frequency of the output frequency series.
        If None, it will be set to delta_f.
    f_final : {None, float}, optional
        The ending frequency of the output frequency series.
        If None, it will be set to the frequency at which the amplitude is 
        1/1000 of the peak amplitude.

    Returns
    -------
    hplustilde: FrequencySeries
        The plus phase of the ringdown in frequency domain.
    hcrosstilde: FrequencySeries
        The cross phase of the ringdown in frequency domain.
    """

    input_params = props_ringdown(template,**kwargs)

    # get required args
    try:
        f_0 = input_params['f_0']
    except KeyError:
        raise ValueError('f_0 is required')
    try:
        tau = input_params['tau']
    except KeyError:
        raise ValueError('tau is required')
    # get optional args
    # the following have defaults, and so will be populated
    t_0 = input_params.pop('t_0')
    phi_0 = input_params.pop('phi_0')
    amp = input_params.pop('amp')
    # the following may not be in input_params
    delta_f = input_params.pop('delta_f', None)
    f_lower = input_params.pop('f_lower', None)
    f_final = input_params.pop('f_final', None)

    if delta_f is None:
        delta_f = 1. / qnm_time_decay(tau, 1./1000)
    if f_lower is None:
        f_lower = delta_f
        kmin = 0
    else:
        kmin = int(f_lower / delta_f)
    if f_final is None:
        f_final = qnm_freq_decay(f_0, tau, 1./1000)
    kmax = int(f_final / delta_f) + 1

    pi = numpy.pi
    two_pi = 2 * numpy.pi
    pi_sq = numpy.pi * numpy.pi

    freqs = numpy.arange(kmin, kmax)*delta_f

    denominator = 1 + (4j * pi * freqs * tau) - (4 * pi_sq * ( freqs*freqs - f_0*f_0) * tau*tau)
    norm = amp * tau / denominator
    if t_0 != 0:
        time_shift = numpy.exp(-1j * two_pi * freqs * t_0) 
        norm *= time_shift

    hp_tilde = norm * ( (1 + 2j * pi * freqs * tau) * numpy.cos(phi_0) - two_pi * f_0 * tau * numpy.sin(phi_0) )
    hc_tilde = norm * ( (1 + 2j * pi * freqs * tau) * numpy.sin(phi_0) + two_pi * f_0 * tau * numpy.cos(phi_0) )

    hplustilde = FrequencySeries(zeros(kmax, dtype=complex128), delta_f=delta_f)
    hcrosstilde = FrequencySeries(zeros(kmax, dtype=complex128), delta_f=delta_f)
    hplustilde.data[kmin:kmax] = hp_tilde
    hcrosstilde.data[kmin:kmax] = hc_tilde

    return hplustilde, hcrosstilde

def calculate_acf(data, delta_t=1.0, unbiased=False):
    """ Calculates the autocorrelation function (ACF) and returns the one-sided
    ACF.

    The ACF is defined as the autocovariance divided by the variance. The ACF
    can be estimated using

        \hat{R}(k) = \frac{1}{\left( n \sigma^{2}} \right) \sum_{t=1}^{n-k} \left( X_{t} - \mu \right) \left( X_{t+k} - \mu \right) 

    Where \hat{R}(k) is the ACF, X_{t} is the data series at time t, \mu is the
    mean of X_{t}, and \sigma^{2} is the variance of X_{t}.

    Parameters
    -----------
    data : {TimeSeries, numpy.array}
        A TimeSeries or numpy.array of data.
    delta_t : float
        The time step of the data series if it is not a TimeSeries instance.
    unbiased : bool
        If True the normalization of the autocovariance function is n-k
        instead of n. This is called the unbiased estimation of the
        autocovariance. Note that this does not mean the ACF is unbiased.

    Returns
    -------
    acf : numpy.array
        If data is a TimeSeries then acf will be a TimeSeries of the
        one-sided ACF. Else acf is a numpy.array.
    """

    # if given a TimeSeries instance then get numpy.array
    if isinstance(data, TimeSeries):
        y = data.numpy()
        delta_t = data.delta_t
    else:
        y = data

    # FFT data minus the mean
    fdata = TimeSeries(y-y.mean(), delta_t=delta_t).to_frequencyseries()

    # correlate
    # do not need to give the congjugate since correlate function does it
    cdata = FrequencySeries(zeros(len(fdata), dtype=numpy.complex64),
                           delta_f=fdata.delta_f, copy=False)
    correlate(fdata, fdata, cdata)

    # IFFT correlated data to get unnormalized autocovariance time series
    acf = cdata.to_timeseries()

    # normalize the autocovariance
    # note that dividing by acf[0] is the same as ( y.var() * len(acf) )
    if unbiased:
        acf /= ( y.var() * numpy.arange(len(acf), 0, -1) )
    else:
        acf /= acf[0]

    # return