search

eqcorrscan / EQcorrscan

Earthquake detection and analysis in Python.
https://eqcorrscan.readthedocs.io/en/latest/
Other
166 stars 86 forks source link

Incorrect filter design for filter gain removal in mag_calc.py #376

Closed TomWinder closed 4 years ago

TomWinder commented 4 years ago

When correcting for gain of pre_filt in mag_calc, a digital filter design is used rather than analog, which is not appropriate to look up the filter gain for a given frequency.

In ObsPy the iirfilter designed to apply a butterworth bandpass filter is digital, so that it is suitable to be applied using scipy.signal.sosfilt().

To look up the filter gain for a given frequency (and correct the amplitude observation for this) an analog filter should be used (as expected by the ObsPy response functions).

Current code: 

z, p, k = iirfilter(corners, [lowcut / (0.5 * tr.stats.sampling_rate),
                           highcut / (0.5 * tr.stats.sampling_rate)],
                           btype='band', ftype='butter', output='zpk')
filt_paz = {'poles': list(p), 'zeros': list(z), 'gain': k,
                  'sensitivity': 1.0}
amplitude /= (paz_2_amplitude_value_of_freq_resp(
                    filt_paz, 1 / period) * filt_paz['sensitivity'])

Correct filter design:

z, p, k = iirfilter(corners, [lowcut * 2 * np.pi, highcut * 2 * np.pi],    # For analog filter, Wn in rad/s
                          btype='band', ftype='butter', output='zpk', analog=True)

See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.iirfilter.html

calum-chamberlain commented 4 years ago

Morning, thanks for looking at this. It looks like you are correct - apologies for this, I clearly didn't look into paz_2_amplitude_value_of_freq_resp properly. Is the analog equivalent actually equivalent, or should we be correcting using something like:

from scipy.signal import iirfilter, sosfreqz

sos = iirfilter(corners, [lowcut / (0.5 * tr.stats.sampling_rate),
                                   highcut / (0.5 * tr.stats.sampling_rate)],
                     btype='band', ftype='butter', output='sos')
_, gain = sosfreqz(sos, worN=[1 / period], fs=tr.stats.sampling_rate)
amplitude /= np.abs(gain)
calum-chamberlain commented 4 years ago

@twinder36 I had a go at fixing this using the method I suggested above using the digital filter in the PR here: #378

It would be great to get your feedback on that.

TomWinder commented 4 years ago

No worries @calum-chamberlain, I've been going through all of this with a fine tooth-comb the past few days to finalise the magnitudes module of QuakeMigrate - the magnitudes code in EQcorrscan has provided a useful reference!

Having played with the two options a bit it seems there is generally a negligible difference between using the equivalent analog response and reading it from sosfreqz, but yours is the more correct approach. Looks good to me.

calum-chamberlain commented 4 years ago

Great, I'm glad it hasn't actually screwed up a workflow - thanks for reporting the issue!

QuakeMigrate looks like its coming along well - what more are you guys planning with it?