anomalously expensive to resample 300s 500Hz trace

Hi

(tested on obspy 1.0.3, linux)

Trying to resample a 500 Hz 300s trace in ObsPy I'm amazed that the opeation takes 16+ seconds. Funny thing is that if I first slice the trace to 299s or less, operation goes in a split second so something seems fishy here.

Attached is the data file and the script I've used to test it, on my machine running the script produces the output:

:~> test_resample.py 300
reading trace took 0.010s
1 Trace(s) in Stream:
HU.HUSM.00.HHE | 2012-02-24T00:00:00.000000Z - 2012-02-24T00:05:00.000000Z | 500.0 Hz, 150001 samples
slicing trace took 0.001s
HU.HUSM.00.HHE | 2012-02-24T00:00:00.000000Z - 2012-02-24T00:05:00.000000Z | 500.0 Hz, 150001 samples
resampling trace took 16.196s
HU.HUSM.00.HHE | 2012-02-24T00:00:00.000000Z - 2012-02-24T00:04:59.990000Z | 100.0 Hz, 30000 samples
:~> test_resample.py 299
reading trace took 0.010s
1 Trace(s) in Stream:
HU.HUSM.00.HHE | 2012-02-24T00:00:00.000000Z - 2012-02-24T00:05:00.000000Z | 500.0 Hz, 150001 samples
slicing trace took 0.001s
HU.HUSM.00.HHE | 2012-02-24T00:00:00.000000Z - 2012-02-24T00:04:59.000000Z | 500.0 Hz, 149501 samples
resampling trace took 0.505s
HU.HUSM.00.HHE | 2012-02-24T00:00:00.000000Z - 2012-02-24T00:04:58.990000Z | 100.0 Hz, 29900 samples
:~>

Best \p

test_resample.py (761 Bytes)

trace.mseed (189 KB)

Hi Peter,

This might be one of these cases where the used numpy FFT is super slow (there are a few number of samples where this is the case - we already deal with this in parts of our code base but probably not for the resample() method).

In general I would use the sinc-based Lanczos interpolation for any operation that changes the sampling rate. It is also linear in time with the number of samples.

tr.interpolate(sampling_rate=100, method=“lanczos”, a=12, window=“blackman”)

The a value determines the half-width of the sinc kernel in terms of original samples. This also defines the extent of any boundary effects. The larger the a, the better and slower the interpolation but there is almost no error with a value of 12.

Cheers!

Lion

Hi Lion

Will this properly handle aliasing when downsampling (I assume that the resample() method of the trace object do but perhaps this is not the case?)

regP

Hi Peter,

no - but the anti-aliasing filter of the normal .resample() method is also not the greatest to be honest. We’ve been meaning to improve this for a while but there is no one-size fits all solution for this type of thing.

If a zerophase filter is acceptable to you: I personally like to do the following:

    def zerophase_chebychev_lowpass_filter(trace, freqmax):
        """
        Custom Chebychev type two zerophase lowpass filter useful for
        decimation filtering.

        This filter is stable up to a reduction in frequency with a factor of
        10. If more reduction is desired, simply decimate in steps.

        Partly based on a filter in ObsPy.

        :param trace: The trace to be filtered.
        :param freqmax: The desired lowpass frequency.
        """
        # rp - maximum ripple of passband, rs - attenuation of stopband
        rp, rs, order = 1, 96, 1e99
        ws = freqmax / (trace.stats.sampling_rate * 0.5)  # stop band frequency
        wp = ws  # pass band frequency

        while True:
            if order <= 12:
                break
            wp *= 0.99
            order, wn = signal.cheb2ord(wp, ws, rp, rs, analog=0)

        b, a = signal.cheby2(order, rs, wn, btype="low", analog=0, output="ba")

        # Apply twice to get rid of the phase distortion.
        trace.data = signal.filtfilt(b, a, trace.data)

signal is scipy.signal. This attenuates the stop-band by like 200 dB which is plenty. The downside is that it does not have the sharpest transition between pass and stop-band.

Hope it helps,

Lion

Well, probably needs to think a bit more about this myself, for now thanks for the insight and pointer

regP