Filter frequencies << sampling frequency are unstable

[QuLogic]

Hi,

I just found a possible question about butterworth-bandpass filtering. I download a seismogram from IRIS and try to filter it with a butterworth-bandpass filter from 0.02hz to 0.2hz, but found the result is incorrect with some extremely large number after filtering. I guess it might be caused by some kind of unstability of filtering algorithm. Here is my code:

from obspy.fdsn import Client
from obspy import UTCDateTime
from obspy.taup.taup import getTravelTimes
from obspy.core import AttribDict
client = Client('iris')
starttime = UTCDateTime('2014-08-24T10:20:44.0600')
p_arrival_time = 376.945875099  # p phase arrival time
st = client.get_waveforms(network='AK', station='BPAW', location='*', channel='BHZ', 
                          starttime=starttime+p_arrival_time-60, endtime=starttime+p_arrival_time+180,
                          attach_response=True)
st.remove_response(output='disp')  # remove the instrument response
st.detrend(type='demean')  # remove the mean
st.detrend(type='linear')  # remove the linear trend
st.taper(type='hann', max_percentage=0.05)  # taper
# st.resample(20)  # the decreasing samping rate can make filtering work well
st.filter('bandpass',freqmin=0.02, freqmax=0.2, zerophase=True)  # filter
st.write('tmp.sac',format='sac')

BTW, when I decrease the sampling rate from original 50hz to 20hz, the filtering will be good. I wonder whether it might also be related to sampling rate.

I also find a similar question about scipy:
http://stackoverflow.com/questions/21862777/bandpass-butterworth-filter-frequencies-in-scipy
Since the filtering in obspy is actually based on scipy, they possibly point to one same question....

Thx

[QuLogic]

Based on this answer, I think we should a) add a warning if the filtering frequencies are very small compared to the sampling rate, and b) switch to sosfilt once a new scipy is available.

[krischer]

I think we should backport the second-order filter sections from scipy. This has been an annoying issue for a long time and we cannot simply depend on scipy 0.16 once available.

That should be fairly quick to do and incredibly useful.

[QuLogic]

It looks like we'd need sosfilt, zpk2sos and all filters need to support output='zpk'. I'm not sure when zpk output was added? Also, zpk2sos is pretty long...

[QuLogic]

Well, there you go...

But you'd need to check whether zpk form is supported all the way back in scipy 0.7.2.

[krischer]

Nice :-) The results are quite spectacular. Its possible to choose much higher order filters than before. Eventually the new filters still become unstable but that is to be expected and it is now MUCH better.

Left is before, right is after. The title are the settings for the bandpass filter.

[megies]

✨

[megies]

Two questions..

• the idea is that this is considered a bug fix and goes into 0.10.x?
• should we make it a compatibility layer, trying to use scipy (for >= 0.16, for sake of upstream bug fixes/improvements) and only using our copied code for older scipy?

[krischer]

No to both in my opinion.

This is not a bugfix but rather a different algorithm.

I would also tend to avoid having a compatibility layer and rather just merge upstream fixes into our backport if we notice them. Otherwise the danger is that the same version of ObsPy with two different scipy yields very different results.

[ThomasLecocq]

👍 wooooow ! I can finally remove the low+high pass two lines in msnoise :)

[megies]

This is not a bugfix but rather a different algorithm.

Then the base branch would have to be changed..

rather just merge upstream fixes into our backport if we notice them.

means a maintenance overhead..

the same version of ObsPy with two different scipy yields very different results.

That's one way to look at it.. you could also say, having a newer scipy would immediately bring it's benefits without having to wait for a new obspy release. Also, that statement would apply to any of our dependencies.

[QuLogic]

the idea is that this is considered a bug fix and goes into 0.10.x?

Sorry, I snuck this one onto the releases target... 😛

should we make it a compatibility layer, trying to use scipy (for >= 0.16, for sake of upstream bug fixes/improvements) and only using our copied code for older scipy?

It already does so. zpk2sos and such are really long; I don't feel like maintaining them unless really necessary.

[krischer]

@megies @QuLogic Alright...point taken :-) Let's have a compatibility layer.

But this should still not got into releases I think as the change is fairly serious and will affect everyone.

[krischer]

This PR works fine with scipy 0.9.0: http://tests.obspy.org/24261/

Assuming #1031 is ok with everyone and we don't merge it in releases it should be good to go.

[XStargate]

After testing the backport sosfilt from scipy 0.16.0, the filtering waveform is the same with waveform filtered by sac even for some extremely small frequency bandpass (e.g. 0.02-0.2hz here). The first plot is filtered by sac and the second one is filtered by obspy sosfilt

The difference between two filtering is as small as 1e-9 which can be ignored.

[megies]

Interesting.. thanks for the feedback.