def normal(start, end, seed=0):
    """ Generate data with a white Gaussian (normal) distribution

    Parameters
    ----------
    start_time : int
        Start time in GPS seconds to generate noise
    end_time : int
        End time in GPS seconds to generate nosie
    seed : {None, int}
        The seed to generate the noise.

    Returns
    --------
    noise : TimeSeries
        A TimeSeries containing gaussian noise
    """
    # This is reproduceable because we used fixed seeds from known values
    s = int(start / BLOCK_SIZE)
    e = int(end / BLOCK_SIZE)

    # The data evenly divides so the last block would be superfluous
    if end % BLOCK_SIZE == 0:
        e -= 1

    sv = RandomState(seed).randint(-2**50, 2**50)
    data = numpy.concatenate([block(i + sv) for i in numpy.arange(s, e + 1, 1)])
    ts = TimeSeries(data, delta_t=1.0 / SAMPLE_RATE, epoch=start)
    return ts.time_slice(start, end)

def line_model(freq, data, tref, amp=1, phi=0):
    """ Simple time-domain model for a frequency line.

    Parameters
    ----------
    freq: float
        Frequency of the line.
    data: pycbc.types.TimeSeries
        Reference data, to get delta_t, start_time, duration and sample_times.
    tref: float
        Reference time for the line model.
    amp: {1., float}, optional
        Amplitude of the frequency line.
    phi: {0. float}, optional
        Phase of the frequency line (radians).

    Returns
    -------
    freq_line: pycbc.types.TimeSeries
        A timeseries of the line model with frequency 'freq'. The returned
        data are complex to allow measuring the amplitude and phase of the
        corresponding frequency line in the strain data. For extraction, use
        only the real part of the data.
    """
    freq_line = TimeSeries(zeros(len(data)), delta_t=data.delta_t,
                           epoch=data.start_time)

    times = data.sample_times - float(tref)
    alpha = 2 * numpy.pi * freq * times + phi
    freq_line.data = amp * numpy.exp(1.j * alpha)

    return freq_line

def noise_from_psd(length, delta_t, psd, seed=None):
    """ Create noise with a given psd.

    Return noise with a given psd. Note that if unique noise is desired
    a unique seed should be provided.

    Parameters
    ----------
    length : int
        The length of noise to generate in samples.
    delta_t : float
        The time step of the noise.
    psd : FrequencySeries
        The noise weighting to color the noise.
    seed : {0, int}
        The seed to generate the noise.

    Returns
    --------
    noise : TimeSeries
        A TimeSeries containing gaussian noise colored by the given psd.
    """
    noise_ts = TimeSeries(zeros(length), delta_t=delta_t)

    if seed is None:
        seed = numpy.random.randint(2**32)

    randomness = lal.gsl_rng("ranlux", seed)

    N = int (1.0 / delta_t / psd.delta_f)
    n = N/2+1
    stride = N/2

    if n > len(psd):
        raise ValueError("PSD not compatible with requested delta_t")

    psd = (psd[0:n]).lal()
    psd.data.data[n-1] = 0

    segment = TimeSeries(zeros(N), delta_t=delta_t).lal()
    length_generated = 0

    SimNoise(segment, 0, psd, randomness)
    while (length_generated < length):
        if (length_generated + stride) < length:
            noise_ts.data[length_generated:length_generated+stride] = segment.data.data[0:stride]
        else:
            noise_ts.data[length_generated:length] = segment.data.data[0:length-length_generated]

        length_generated += stride
        SimNoise(segment, stride, psd, randomness)

    return noise_ts

 def test_injection_presence(self):
     """Verify presence of signals at expected times"""
     injections = InjectionSet(self.inj_file.name)
     for det in self.detectors:
         for inj in self.injections:
             ts = TimeSeries(numpy.zeros(10 * self.sample_rate),
                             delta_t=1/self.sample_rate,
                             epoch=lal.LIGOTimeGPS(inj.end_time - 5),
                             dtype=numpy.float64)
             injections.apply(ts, det.name)
             max_amp, max_loc = ts.abs_max_loc()
             # FIXME could test amplitude and time more precisely
             self.assertTrue(max_amp > 0 and max_amp < 1e-10)
             time_error = ts.sample_times.numpy()[max_loc] - inj.end_time
             self.assertTrue(abs(time_error) < 1.5 * self.earth_time)

 def test_injection_absence(self):
     """Verify absence of signals outside known injection times"""
     clear_times = [
         self.injections[0].end_time - 86400,
         self.injections[-1].end_time + 86400
     ]
     injections = InjectionSet(self.inj_file.name)
     for det in self.detectors:
         for epoch in clear_times:
             ts = TimeSeries(numpy.zeros(10 * self.sample_rate),
                             delta_t=1/self.sample_rate,
                             epoch=lal.LIGOTimeGPS(epoch),
                             dtype=numpy.float64)
             injections.apply(ts, det.name)
             max_amp, max_loc = ts.abs_max_loc()
             self.assertEqual(max_amp, 0)

    def __init__(self, frame_src, 
                       channel_name,
                       start_time,
                       max_buffer=2048, 
                       force_update_cache=True,
                       increment_update_cache=None):
        """ Create a rolling buffer of frame data

        Parameters
        ---------
        frame_src: str of list of strings
            Strings that indicate where to read from files from. This can be a
        list of frame files, a glob, etc.
        channel_name: str
            Name of the channel to read from the frame files
        start_time: 
            Time to start reading from.
        max_buffer: {int, 2048}, Optional
            Length of the buffer in seconds
        """
        self.frame_src = frame_src
        self.channel_name = channel_name
        self.read_pos = start_time
        self.force_update_cache = force_update_cache
        self.increment_update_cache = increment_update_cache

        self.update_cache()
        self.channel_type, self.raw_sample_rate = self._retrieve_metadata(self.stream, self.channel_name)

        raw_size = self.raw_sample_rate * max_buffer
        self.raw_buffer = TimeSeries(zeros(raw_size, dtype=numpy.float64),
                                     copy=False,
                                     epoch=start_time - max_buffer,
                                     delta_t=1.0/self.raw_sample_rate)

def interpolate_complex_frequency(series, delta_f, zeros_offset=0, side='right'):
    """Interpolate complex frequency series to desired delta_f.

    Return a new complex frequency series that has been interpolated to the
    desired delta_f.

    Parameters
    ----------
    series : FrequencySeries
        Frequency series to be interpolated.
    delta_f : float
        The desired delta_f of the output
    zeros_offset : optional, {0, int}
        Number of sample to delay the start of the zero padding
    side : optional, {'right', str}
        The side of the vector to zero pad
        
    Returns
    -------
    interpolated series : FrequencySeries
        A new FrequencySeries that has been interpolated.
    """
    new_n = int( (len(series)-1) * series.delta_f / delta_f + 1)
    samples = numpy.arange(0, new_n) * delta_f
    old_N = int( (len(series)-1) * 2 )
    new_N = int( (new_n - 1) * 2 )
    time_series = TimeSeries(zeros(old_N), delta_t =1.0/(series.delta_f*old_N),
                             dtype=real_same_precision_as(series))
                             
    ifft(series, time_series)

    time_series.roll(-zeros_offset)
    time_series.resize(new_N)
    
    if side == 'left':
        time_series.roll(zeros_offset + new_N - old_N)
    elif side == 'right':
        time_series.roll(zeros_offset)

    out_series = FrequencySeries(zeros(new_n), epoch=series.epoch,
                           delta_f=delta_f, dtype=series.dtype)
    fft(time_series, out_series)

    return out_series

class DataBuffer(object):
    """ A linear buffer that acts as a FILO for reading in frame data
    """
    def __init__(self, frame_src, 
                       channel_name,
                       start_time,
                       max_buffer=2048):
        """ Create a rolling buffer of frame data

        Parameters
        ---------
        frame_src: str of list of strings
            Strings that indicate where to read from files from. This can be a
        list of frame files, a glob, etc.
        channel_name: str
            Name of the channel to read from the frame files
        start_time: 
            Time to start reading from.
        max_buffer: {int, 2048}, Optional
            Length of the buffer in seconds
        """
        self.frame_src = frame_src
        self.channel_name = channel_name
        self.read_pos = start_time

        self.update_cache()
        self.channel_type, self.sample_rate = self._retrieve_metadata(self.stream, self.channel_name)

        raw_size = self.sample_rate * max_buffer
        self.raw_buffer = TimeSeries(zeros(raw_size, dtype=numpy.float64),
                                     copy=False,
                                     epoch=start_time - max_buffer,
                                     delta_t=1.0/self.sample_rate)

    def update_cache(self):
        """ Reset the lal cache. This can be used to update the cache if the 
        result may change due to more files being added to the filesystem, 
        for example.
        """
        cache = locations_to_cache(self.frame_src)
        stream = lalframe.FrStreamCacheOpen(cache)
        self.stream = stream

    def _retrieve_metadata(self, stream, channel_name):
        """ Retrieve basic metadata by reading the first file in the cache
    
        Parameters
        ----------
        stream: lal stream object
            Stream containing a channel we want to learn about
        channel_name: str
            The name of the channel we want to know the dtype and sample rate of

        Returns
        -------
        channel_type: lal type enum
            Enum value which indicates the dtype of the channel
        sample_rate: int
            The sample rate of the data within this channel
        """
        data_length = lalframe.FrStreamGetVectorLength(channel_name, stream)
        channel_type = lalframe.FrStreamGetTimeSeriesType(channel_name, stream)
        create_series_func = _fr_type_map[channel_type][2]
        get_series_metadata_func = _fr_type_map[channel_type][3]
        series = create_series_func(channel_name, stream.epoch, 0, 0,
                            lal.ADCCountUnit, 0)
        get_series_metadata_func(series, stream)
        return channel_type, int(1.0/series.deltaT)        

    def _read_frame(self, blocksize):
        """ Try to read the block of data blocksize seconds long

        Parameters
        ----------
        blocksize: int
            The number of seconds to attempt to read from the channel

        Returns
        -------
        data: TimeSeries
            TimeSeries containg 'blocksize' seconds of frame data

        Raises
        ------
        RuntimeError:
            If data cannot be read for any reason
        """
        try:
            read_func = _fr_type_map[self.channel_type][0]
            dtype = _fr_type_map[self.channel_type][1]
            data = read_func(self.stream, self.channel_name, self.read_pos, blocksize, 0)
            return TimeSeries(data.data.data, delta_t=data.deltaT,
                              epoch=self.read_pos, 
                              dtype=dtype)     
        except:
            raise RuntimeError('Cannot read requested frame data') 

    def null_advance(self, blocksize):
        """ Advance and insert zeros

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

class DataBuffer(object):

    """A linear buffer that acts as a FILO for reading in frame data
    """

    def __init__(self, frame_src, 
                       channel_name,
                       start_time,
                       max_buffer=2048, 
                       force_update_cache=True,
                       increment_update_cache=None):
        """ Create a rolling buffer of frame data

        Parameters
        ---------
        frame_src: str of list of strings
            Strings that indicate where to read from files from. This can be a
        list of frame files, a glob, etc.
        channel_name: str
            Name of the channel to read from the frame files
        start_time: 
            Time to start reading from.
        max_buffer: {int, 2048}, Optional
            Length of the buffer in seconds
        """
        self.frame_src = frame_src
        self.channel_name = channel_name
        self.read_pos = start_time
        self.force_update_cache = force_update_cache
        self.increment_update_cache = increment_update_cache

        self.update_cache()
        self.channel_type, self.raw_sample_rate = self._retrieve_metadata(self.stream, self.channel_name)

        raw_size = self.raw_sample_rate * max_buffer
        self.raw_buffer = TimeSeries(zeros(raw_size, dtype=numpy.float64),
                                     copy=False,
                                     epoch=start_time - max_buffer,
                                     delta_t=1.0/self.raw_sample_rate)

    def update_cache(self):
        """Reset the lal cache. This can be used to update the cache if the 
        result may change due to more files being added to the filesystem, 
        for example.
        """
        cache = locations_to_cache(self.frame_src)
        stream = lalframe.FrStreamCacheOpen(cache)
        self.stream = stream

    def _retrieve_metadata(self, stream, channel_name):
        """Retrieve basic metadata by reading the first file in the cache
    
        Parameters
        ----------
        stream: lal stream object
            Stream containing a channel we want to learn about
        channel_name: str
            The name of the channel we want to know the dtype and sample rate of

        Returns
        -------
        channel_type: lal type enum
            Enum value which indicates the dtype of the channel
        sample_rate: int
            The sample rate of the data within this channel
        """
        data_length = lalframe.FrStreamGetVectorLength(channel_name, stream)
        channel_type = lalframe.FrStreamGetTimeSeriesType(channel_name, stream)
        create_series_func = _fr_type_map[channel_type][2]
        get_series_metadata_func = _fr_type_map[channel_type][3]
        series = create_series_func(channel_name, stream.epoch, 0, 0,
                            lal.ADCCountUnit, 0)
        get_series_metadata_func(series, stream)
        return channel_type, int(1.0/series.deltaT)

    def _read_frame(self, blocksize):
        """Try to read the block of data blocksize seconds long

        Parameters
        ----------
        blocksize: int
            The number of seconds to attempt to read from the channel

        Returns
        -------
        data: TimeSeries
            TimeSeries containg 'blocksize' seconds of frame data

        Raises
        ------
        RuntimeError:
            If data cannot be read for any reason
        """
        try:
            read_func = _fr_type_map[self.channel_type][0]
            dtype = _fr_type_map[self.channel_type][1]
            data = read_func(self.stream, self.channel_name,
                             self.read_pos, int(blocksize), 0)
            return TimeSeries(data.data.data, delta_t=data.deltaT,
                              epoch=self.read_pos, 
#.........这里部分代码省略.........

def get_td_from_freqtau(template=None, taper=None, **kwargs):
    """Return time domain ringdown with all the modes specified.

    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.
    taper: {None, float}, optional
        Tapering at the beginning of the waveform with duration taper * tau.
        This option is recommended with timescales taper=1./2 or 1. for
        time-domain ringdown-only injections.
        The abrupt turn on of the ringdown can cause issues on the waveform
        when doing the fourier transform to the frequency domain. Setting
        taper will add a rapid ringup with timescale tau/10.
        Each mode and overtone will have a different taper depending on its tau,
        the final taper being the superposition of all the tapers.
    lmns : list
        Desired lmn modes as strings (lm modes available: 22, 21, 33, 44, 55).
        The n specifies the number of overtones desired for the corresponding
        lm pair (maximum n=8).
        Example: lmns = ['223','331'] are the modes 220, 221, 222, and 330
    f_lmn: float
        Central frequency of the lmn overtone, as many as number of modes.
    tau_lmn: float
        Damping time of the lmn overtone, as many as number of modes.
    amp220 : float
        Amplitude of the fundamental 220 mode.
    amplmn : float
        Fraction of the amplitude of the lmn overtone relative to the
        fundamental mode, as many as the number of subdominant modes.
    philmn : float
        Phase of the lmn overtone, as many as the number of modes. Should also
        include the information from the azimuthal angle (phi + m*Phi).
    inclination : {None, float}, optional
        Inclination of the system in radians. If None, the spherical harmonics
        will be set to 1.
    delta_t : {None, float}, optional
        The time step used to generate the ringdown.
        If None, it will be set to the inverse of the frequency at which the
        amplitude is 1/1000 of the peak amplitude (the minimum of all modes).
    t_final : {None, float}, optional
        The ending time of the output frequency series.
        If None, it will be set to the time at which the amplitude
        is 1/1000 of the peak amplitude (the maximum of all modes).

    Returns
    -------
    hplustilde: FrequencySeries
        The plus phase of a ringdown with the lm modes specified and
        n overtones in frequency domain.
    hcrosstilde: FrequencySeries
        The cross phase of a ringdown with the lm modes specified and
        n overtones in frequency domain.
    """

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

    # Get required args
    f_0, tau = lm_freqs_taus(**input_params)
    lmns = input_params['lmns']
    for lmn in lmns:
        if int(lmn[2]) == 0:
            raise ValueError('Number of overtones (nmodes) must be greater '
                             'than zero.')
    # following may not be in input_params
    inc = input_params.pop('inclination', None)
    delta_t = input_params.pop('delta_t', None)
    t_final = input_params.pop('t_final', None)

    if not delta_t:
        delta_t = lm_deltat(f_0, tau, lmns)
    if not t_final:
        t_final = lm_tfinal(tau, lmns)

    kmax = int(t_final / delta_t) + 1
    # Different overtones will have different tapering window-size
    # Find maximum window size to create long enough output vector
    if taper:
        taper_window = int(taper*max(tau.values())/delta_t)
        kmax += taper_window

    outplus = TimeSeries(zeros(kmax, dtype=float64), delta_t=delta_t)
    outcross = TimeSeries(zeros(kmax, dtype=float64), delta_t=delta_t)
    if taper:
        start = - taper * max(tau.values())
        outplus._epoch, outcross._epoch = start, start

    for lmn in lmns:
        l, m, nmodes = int(lmn[0]), int(lmn[1]), int(lmn[2])
        hplus, hcross = get_td_lm(freqs=f_0, taus=tau, l=l, m=m, nmodes=nmodes,
                             taper=taper, inclination=inc, delta_t=delta_t,
                             t_final=t_final, **input_params)
        if not taper:
            outplus.data += hplus.data
            outcross.data += hcross.data
        else:
            outplus = taper_shift(hplus, outplus)
            outcross = taper_shift(hcross, outcross)

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

def get_td_lm(template=None, taper=None, **kwargs):
    """Return time domain lm mode with the given number of overtones.

    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.
    taper: {None, float}, optional
        Tapering at the beginning of the waveform with duration taper * tau.
        This option is recommended with timescales taper=1./2 or 1. for
        time-domain ringdown-only injections.
        The abrupt turn on of the ringdown can cause issues on the waveform
        when doing the fourier transform to the frequency domain. Setting
        taper will add a rapid ringup with timescale tau/10.
        Each overtone will have a different taper depending on its tau, the
        final taper being the superposition of all the tapers.
    freqs : dict
        {lmn:f_lmn} Dictionary of the central frequencies for each overtone,
        as many as number of modes.
    taus : dict
        {lmn:tau_lmn} Dictionary of the damping times for each overtone,
        as many as number of modes.
    l : int
        l mode (lm modes available: 22, 21, 33, 44, 55).
    m : int
        m mode (lm modes available: 22, 21, 33, 44, 55).
    nmodes: int
        Number of overtones desired (maximum n=8)
    amp220 : float
        Amplitude of the fundamental 220 mode, needed for any lm.
    amplmn : float
        Fraction of the amplitude of the lmn overtone relative to the
        fundamental mode, as many as the number of subdominant modes.
    philmn : float
        Phase of the lmn overtone, as many as the number of modes. Should also
        include the information from the azimuthal angle (phi + m*Phi).
    inclination : {None, float}, optional
        Inclination of the system in radians for the spherical harmonics.
    delta_t : {None, float}, optional
        The time step used to generate the ringdown.
        If None, it will be set to the inverse of the frequency at which the
        amplitude is 1/1000 of the peak amplitude (the minimum of all modes).
    t_final : {None, float}, optional
        The ending time of the output time series.
        If None, it will be set to the time at which the amplitude is
        1/1000 of the peak amplitude (the maximum of all modes).

    Returns
    -------
    hplus: TimeSeries
        The plus phase of a lm mode with overtones (n) in time domain.
    hcross: TimeSeries
        The cross phase of a lm mode with overtones (n) in time domain.
    """

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

    # Get required args
    amps, phis = lm_amps_phases(**input_params)
    f_0 = input_params.pop('freqs')
    tau = input_params.pop('taus')
    inc = input_params.pop('inclination', None)
    l, m = input_params.pop('l'), input_params.pop('m')
    nmodes = input_params.pop('nmodes')
    if int(nmodes) == 0:
        raise ValueError('Number of overtones (nmodes) must be greater '
                         'than zero.')
    # The following may not be in input_params
    delta_t = input_params.pop('delta_t', None)
    t_final = input_params.pop('t_final', None)

    if not delta_t:
        delta_t = lm_deltat(f_0, tau, ['%d%d%d' %(l,m,nmodes)])
    if not t_final:
        t_final = lm_tfinal(tau, ['%d%d%d' %(l, m, nmodes)])

    kmax = int(t_final / delta_t) + 1
    # Different overtones will have different tapering window-size
    # Find maximum window size to create long enough output vector
    if taper:
        taper_window = int(taper*max(tau.values())/delta_t)
        kmax += taper_window

    outplus = TimeSeries(zeros(kmax, dtype=float64), delta_t=delta_t)
    outcross = TimeSeries(zeros(kmax, dtype=float64), delta_t=delta_t)
    if taper:
        start = - taper * max(tau.values())
        outplus._epoch, outcross._epoch = start, start

    for n in range(nmodes):
        hplus, hcross = get_td_qnm(template=None, taper=taper,
                            f_0=f_0['%d%d%d' %(l,m,n)],
                            tau=tau['%d%d%d' %(l,m,n)],
                            phi=phis['%d%d%d' %(l,m,n)],
                            amp=amps['%d%d%d' %(l,m,n)],
                            inclination=inc, l=l, m=m,
                            delta_t=delta_t, t_final=t_final)
        if not taper:
            outplus.data += hplus.data
#.........这里部分代码省略.........

def get_td_qnm(template=None, taper=None, **kwargs):
    """Return a time 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.
    taper: {None, float}, optional
        Tapering at the beginning of the waveform with duration taper * tau.
        This option is recommended with timescales taper=1./2 or 1. for
        time-domain ringdown-only injections.
        The abrupt turn on of the ringdown can cause issues on the waveform
        when doing the fourier transform to the frequency domain. Setting
        taper will add a rapid ringup with timescale tau/10.
    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.
    delta_t : {None, float}, optional
        The time step used to generate the ringdown.
        If None, it will be set to the inverse of the frequency at which the
        amplitude is 1/1000 of the peak amplitude.
    t_final : {None, float}, optional
        The ending time of the output time series.
        If None, it will be set to the time at which the amplitude is
        1/1000 of the peak amplitude.

    Returns
    -------
    hplus: TimeSeries
        The plus phase of the ringdown in time domain.
    hcross: TimeSeries
        The cross phase of the ringdown in time 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 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_t = input_params.pop('delta_t', None)
    t_final = input_params.pop('t_final', None)

    if not delta_t:
        delta_t = 1. / qnm_freq_decay(f_0, tau, 1./1000)
        if delta_t < min_dt:
            delta_t = min_dt
    if not t_final:
        t_final = qnm_time_decay(tau, 1./1000)

    kmax = int(t_final / delta_t) + 1
    times = numpy.arange(kmax) * delta_t
    if inc is not None:
        Y_plus, Y_cross = spher_harms(l, m, inc)
    else:
        Y_plus, Y_cross = 1, 1

    hplus = amp * Y_plus * numpy.exp(-times/tau) * \
                                numpy.cos(two_pi*f_0*times + phi)
    hcross = amp * Y_cross * numpy.exp(-times/tau) * \
                                numpy.sin(two_pi*f_0*times + phi)

    if taper and delta_t < taper*tau:
        taper_window = int(taper*tau/delta_t)
        kmax += taper_window

    outplus = TimeSeries(zeros(kmax), delta_t=delta_t)
    outcross = TimeSeries(zeros(kmax), delta_t=delta_t)

    # If size of tapering window is less than delta_t, do not apply taper.
    if not taper or delta_t > taper*tau:
        outplus.data[:kmax] = hplus
        outcross.data[:kmax] = hcross

        return outplus, outcross

    else:
        taper_hp, taper_hc = apply_taper(delta_t, taper, f_0, tau, amp, phi,
                                                                    l, m, inc)
        start = - taper * tau
        outplus.data[:taper_window] = taper_hp
        outplus.data[taper_window:] = hplus
        outcross.data[:taper_window] = taper_hc
#.........这里部分代码省略.........

    def apply(self, strain, detector_name, f_lower=None, distance_scale=1):
        """Add injections (as seen by a particular detector) to a time series.

        Parameters
        ----------
        strain : TimeSeries
            Time series to inject signals into, of type float32 or float64.
        detector_name : string
            Name of the detector used for projecting injections.
        f_lower : {None, float}, optional
            Low-frequency cutoff for injected signals. If None, use value
            provided by each injection.
        distance_scale: {1, foat}, optional
            Factor to scale the distance of an injection with. The default is
            no scaling.

        Returns
        -------
        None

        Raises
        ------
        TypeError
            For invalid types of `strain`.
        """

        if not strain.dtype in (float32, float64):
            raise TypeError("Strain dtype must be float32 or float64, not " \
                    + str(strain.dtype))

        lalstrain = strain.lal()
        detector = Detector(detector_name)
        earth_travel_time = lal.REARTH_SI / lal.C_SI
        t0 = float(strain.start_time) - earth_travel_time
        t1 = float(strain.end_time) + earth_travel_time

        # pick lalsimulation tapering function
        taper = taper_func_map[strain.dtype]

        # pick lalsimulation injection function
        add_injection = injection_func_map[strain.dtype]

        for inj in self.table:
            # roughly estimate if the injection may overlap with the segment
            end_time = inj.get_time_geocent()
            #CHECK: This is a hack (10.0s); replace with an accurate estimate
            inj_length = 10.0
            eccentricity = 0.0
            polarization = 0.0
            start_time = end_time - 2 * inj_length
            if end_time < t0 or start_time > t1:
               continue

            # compute the waveform time series
            hp, hc = sim.SimBurstSineGaussian(float(inj.q),
                float(inj.frequency),float(inj.hrss),float(eccentricity),
                float(polarization),float(strain.delta_t))
            hp = TimeSeries(hp.data.data[:], delta_t=hp.deltaT, epoch=hp.epoch)
            hc = TimeSeries(hc.data.data[:], delta_t=hc.deltaT, epoch=hc.epoch)
            hp._epoch += float(end_time)
            hc._epoch += float(end_time)
            if float(hp.start_time) > t1:
               continue

            # compute the detector response, taper it if requested
            # and add it to the strain
            #signal = detector.project_wave(
            #        hp, hc, inj.longitude, inj.latitude, inj.polarization)
            signal_lal = hp.astype(strain.dtype).lal()
            if taper_map['TAPER_NONE'] is not None:
                taper(signal_lal.data, taper_map['TAPER_NONE'])
            add_injection(lalstrain, signal_lal, None)

        strain.data[:] = lalstrain.data.data[:]

def get_td_qnm(template=None, delta_t=None, t_lower=None, t_final=None, **kwargs):
    """Return a time 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_t : {None, float}, optional
        The time step used to generate the ringdown.
        If None, it will be set to the inverse of the frequency at which the
        amplitude is 1/100 of the peak amplitude.
    t_lower: {None, float}, optional
        The starting time of the output time series.
        If None, it will be set to delta_t.
    t_final : {None, float}, optional
        The ending time of the output time series.
        If None, it will be set to the time at which the amplitude is 
        1/1000 of the peak amplitude.

    Returns
    -------
    hplus: TimeSeries
        The plus phase of the ringdown in time domain.
    hcross: TimeSeries
        The cross phase of the ringdown in time domain.
    """

    input_params = props_ringdown(template,**kwargs)

    f_0 = input_params['f_0']
    tau = input_params['tau']
    t_0 = input_params['t_0']
    phi_0 = input_params['phi_0']
    Amp = input_params['Amp']
    if delta_t is None:
        delta_t = 1. / qnm_freq_decay(f_0, tau, 1./100)
    if t_lower is None:
        t_lower = delta_t
        kmin = 0
    else:
        kmin=int(t_lower / delta_t)
    if t_final is None:
        t_final = qnm_time_decay(tau, 1./1000)
    kmax = int(t_final / delta_t)
    n = int(t_final / delta_t) + 1

    two_pi = 2 * numpy.pi

    times = numpy.arange(t_lower, t_final, delta_t)

    hp = Amp * numpy.exp(-times/tau) * numpy.cos(two_pi*f_0*times + phi_0)
    hc = Amp * numpy.exp(-times/tau) * numpy.sin(two_pi*f_0*times + phi_0)

    hplus = TimeSeries(zeros(n), delta_t=delta_t)
    hcross = TimeSeries(zeros(n), delta_t=delta_t)
    hplus.data[kmin:kmax] = hp
    hcross.data[kmin:kmax] = hc

    return hplus, hcross

    def apply(self, strain, detector_name, f_lower=None, distance_scale=1,
              simulation_ids=None):
        """Add injections (as seen by a particular detector) to a time series.

        Parameters
        ----------
        strain : TimeSeries
            Time series to inject signals into, of type float32 or float64.
        detector_name : string
            Name of the detector used for projecting injections.
        f_lower : {None, float}, optional
            Low-frequency cutoff for injected signals. If None, use value
            provided by each injection.
        distance_scale: {1, float}, optional
            Factor to scale the distance of an injection with. The default is 
            no scaling. 
        simulation_ids: iterable, optional
            If given, only inject signals with the given simulation IDs.

        Returns
        -------
        None

        Raises
        ------
        TypeError
            For invalid types of `strain`.
        """

        if not strain.dtype in (float32, float64):
            raise TypeError("Strain dtype must be float32 or float64, not " \
                    + str(strain.dtype))

        lalstrain = strain.lal()    
        detector = Detector(detector_name)
        earth_travel_time = lal.REARTH_SI / lal.C_SI
        t0 = float(strain.start_time) - earth_travel_time
        t1 = float(strain.end_time) + earth_travel_time

        # pick lalsimulation injection function
        add_injection = injection_func_map[strain.dtype]

        injections = self.table
        if simulation_ids:
            injections = [inj for inj in injections \
                          if inj.simulation_id in simulation_ids]

        for inj in injections:
            if f_lower is None:
                f_l = inj.f_lower
            else:
                f_l = f_lower

            if inj.numrel_data != None and inj.numrel_data != "":
                # performing NR waveform injection
                # reading Hp and Hc from the frame files
                swigrow = self.getswigrow(inj)
                import lalinspiral
                Hp, Hc = lalinspiral.NRInjectionFromSimInspiral(swigrow,
                                                                strain.delta_t)
                # converting to pycbc timeseries
                hp = TimeSeries(Hp.data.data[:], delta_t=Hp.deltaT,
                                epoch=Hp.epoch)
                hc = TimeSeries(Hc.data.data[:], delta_t=Hc.deltaT,
                                epoch=Hc.epoch)
                hp /= distance_scale
                hc /= distance_scale
                end_time = float(hp.get_end_time())
                start_time = float(hp.get_start_time())
                if end_time < t0 or start_time > t1:
                    continue
            else:
                # roughly estimate if the injection may overlap with the segment
                end_time = inj.get_time_geocent()
                inj_length = sim.SimInspiralTaylorLength(
                    strain.delta_t, inj.mass1 * lal.MSUN_SI,
                    inj.mass2 * lal.MSUN_SI, f_l, 0)
                start_time = end_time - 2 * inj_length
                if end_time < t0 or start_time > t1:
                   continue
                   
                name, phase_order = legacy_approximant_name(inj.waveform)

                # compute the waveform time series
                hp, hc = get_td_waveform(
                    inj, approximant=name, delta_t=strain.delta_t,
                    phase_order=phase_order,
                    f_lower=f_l, distance=inj.distance * distance_scale,
                    **self.extra_args)

                hp._epoch += float(end_time)
                hc._epoch += float(end_time)
                if float(hp.start_time) > t1:
                   continue

            # compute the detector response, taper it if requested
            # and add it to the strain
            signal = detector.project_wave(
                    hp, hc, inj.longitude, inj.latitude, inj.polarization)
            # the taper_timeseries function converts to a LAL TimeSeries
#.........这里部分代码省略.........

def get_td_qnm(template=None, **kwargs):
    """Return a time 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.
    phi_0 : {0, float}, optional
        The initial phase of the ringdown.
    amp : {1, float}, optional
        The amplitude of the ringdown (constant for now).
    delta_t : {None, float}, optional
        The time step used to generate the ringdown.
        If None, it will be set to the inverse of the frequency at which the
        amplitude is 1/1000 of the peak amplitude.
    t_final : {None, float}, optional
        The ending time of the output time series.
        If None, it will be set to the time at which the amplitude is 
        1/1000 of the peak amplitude.

    Returns
    -------
    hplus: TimeSeries
        The plus phase of the ringdown in time domain.
    hcross: TimeSeries
        The cross phase of the ringdown in time 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
    phi_0 = input_params.pop('phi_0')
    amp = input_params.pop('amp')
    # the following may not be in input_params
    delta_t = input_params.pop('delta_t', None)
    t_final = input_params.pop('t_final', None)

    if delta_t is None:
        delta_t = 1. / qnm_freq_decay(f_0, tau, 1./1000)
    if t_final is None:
        t_final = qnm_time_decay(tau, 1./1000)
    kmax = int(t_final / delta_t) + 1

    two_pi = 2 * numpy.pi

    times = numpy.arange(kmax)*delta_t

    hp = amp * numpy.exp(-times/tau) * numpy.cos(two_pi*f_0*times + phi_0)
    hc = amp * numpy.exp(-times/tau) * numpy.sin(two_pi*f_0*times + phi_0)

    hplus = TimeSeries(zeros(kmax), delta_t=delta_t)
    hcross = TimeSeries(zeros(kmax), delta_t=delta_t)
    hplus.data[:kmax] = hp
    hcross.data[:kmax] = hc

    return hplus, hcross