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