Applying FIR filters to long timeseries

Created by: andrew-lundgren

Scipy seems to lack a method for applying FIR filters to very long time series (meaning longer than can be efficiently FFTed). One standard way to do this is the overlap-save method. It goes through the data by chunks, overlapping enough to avoid the corruption from filter wraparound.

I've implemented this using standard numpy and scipy calls. The sample code is below. I think this should be made a method of the TimeSeries class, maybe called .convolve, since it's equivalent to scipy.signal.convolve.

If this is OK, I can prepare a patch to TimeSeries in a more gwpy style. This could be used to close #358 (closed).

def overlap_save(in1, in2, fftlen=16384):
    """Convolve a very long data series with a short one
        using the overlap-save method. Result should be
        nearly identical to scipy.signal.convolve with
        mode='same'.
        
        Inputs:
        in1 (numpy array): long data series
        in2 (numpy array): short data series (e.g. FIR filter)
        fftlen (int, optional): number of points to use for FFT
        
        Output:
        Numpy array with same length as in1"""
    # If data is short enough, just do a single convolve
    if len(in1) <= 2*fftlen:
        return sig.convolve(in1, in2, mode='same')
    
    # Find how much is corrupted by filter wraparound
    pad = int(np.ceil(len(in2)/2.))
    step = int(fftlen) - 2*pad
    # Create room for the result
    result=np.zeros(len(in1), dtype=in1.dtype)
    # Do the first chunk on its own
    result[:fftlen-pad]=sig.fftconvolve(in1[:fftlen-pad], in2, mode='same')
    # Walk through chunk by chunk
    for idx in xrange(2*step, len(in1), step):
        tmp = sig.fftconvolve(in1[idx-fftlen:idx], in2, mode='same')
        result[idx-fftlen+pad:idx-pad] = tmp[pad:-pad]
    # Clean up the last piece which may be anything smaller than step
    # It's important that we made sure to have at least one step in
    # the loop above, otherwise idx is undefined
    len_final = len(in1) - idx
    tmp = sig.fftconvolve(in1[-fftlen:], in2, mode='same')
    result[-len_final-pad:] = tmp[-len_final-pad:]
    
    return result