adjtomo / pysep

Seismogram Extraction and Processing: Seismic data retrieval and record sections
https://pysep.readthedocs.io
MIT License
29 stars 15 forks source link

instrument response removal for very low frequencies #96

Open carltape opened 1 year ago

carltape commented 1 year ago

I'm not sure if this is a pysep issue or obspy issue, since we are using obspy functions within pysep. The testing notebook is response_removal_test.ipynb. This contains more than is needed to demonstrate the issue, which is this: When testing a range of options for removing the instrument response, we cannot seem to achieve the option to directly divide FFT(c(t)) by Ia(f) to get a(f), where c(t) is the time series in counts, Ia(f) is the instrument response to acceleration, and a(f) is the acceleration in the frequency domain.

Here are four subplots (from the notebook above) showing the 0.2-1.0 mHz spectra for II.NNA.00.LHZ derived from a 10-day time series of the 2004 Sumatra earthquake. From the top, they are: 1) instrument response to acceleration (simple flat line) 2) raw 3) acceleration from manual removal of response 4) acceleration via obspy remove_response (via pysep) using the best input options we found (water_level = None, pre_filt = default).

Screen Shot 2023-03-04 at 3 05 29 PM

Here is II.SUR.00.BHZ:

Screen Shot 2023-03-04 at 3 06 12 PM

For this class problem, we would like to achieve the cleanest results for instrument response removal, which seem to be the direct division. We hope that this is possible by using the remove_response function in obspy (and, in turn, pysep), but we can't see how. Perhaps there is some taper or operation that is happening in obspy that is causing the issue. Or perhaps pysep is applying some additional operation. Thanks.

carltape commented 1 year ago

Offering up some improved results (from Aakash) by using the taper_percentage option in pysep.

Here is NNA for using a 20 percent taper (taper_percentage = 0.2): 20_percent_taper

Here is NNA for using a 50 percent taper (taper_percentage = 0.5): 50_percent_taper

Still, the pysep version is slightly worse than the manual removal, so there must be something else happening inside pysep (or remove_response).

carltape commented 1 year ago

Attaching a test script here, since the notebook at the top of this thread is no longer available. test_response_removal.txt