Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Causal instrument response correction #46

Open
dylanmikesell opened this issue May 31, 2020 · 11 comments
Open

Causal instrument response correction #46

dylanmikesell opened this issue May 31, 2020 · 11 comments

Comments

@dylanmikesell
Copy link

I would like to implement the response deconvolution algorithm in this paper from Matt Haney et al. (2012) in BSSA. I was reading through the developer's guide on processing algorithms. This is a published algorithm and I will be translating directly from Matt's MATLAB code. This is an algorithm that ensures causality and has worked great for my purposes of signal correlations across a heterogeneous network of seismic stations. I am curious if this is of use for the SeisIO.jl package and if so, how should I go about adding this (keyword for "type" in remove_resp!). Should I add it to SeisIO.jl/src/Processing/Resp/0_resp_aux.jl?

@jpjones76
Copy link
Owner

It's definitely useful for SeisIO and I encourage it.

Recommendations:

  • Probably best in a new file in src/Processing/Resp/.
  • Definitely best to use a new keyword in translate_resp! or remove_resp!. Does the algorithm flatten to DC?
  • We shouldn't use "type" for a variable name in Julia. Too much potential confusion with Type, typeof, etc.

If you have a Julia implementation of Matt's algorithm, we can work out insertion points for the new function.

@dylanmikesell
Copy link
Author

dylanmikesell commented May 31, 2020

Okay. I am working on the Julia version now. Will probably be late tomorrow or Tuesday before I have something worth sharing.

Working through his Cleveland Volcano example has brought up another issue with SeisData.resp.f0, but I will post that in another issue....if it turns out I can't figure it out.

@dylanmikesell
Copy link
Author

I am back to translating this MATLAB code. Any opinions on interpolation packages?

Interpolations.jl
Dierckx.jl
GridInterpolations.jl
ApproXD.jl

I need to do a 1D nearest neighbor interpolation (to follow Matt's code). I don't see an interpolation dep in SeisIO yet. It will be an interpolation on evenly gridded data. Happy to use whatever people suggest. I am playing with Interpolations.jl right now.

@jpjones76
Copy link
Owner

I'm checking each package. Nothing currently in SeisIO requires interpolation, but a robust interpolation package will eventually be needed for many processing operations. I want to be certain that whatever we use is robust and well-maintained.

One thing that keeps me away from R, and which I feel Julia must actively prevent, is a huge ecosystem of packages with tremendous redundancy. By extension I feel like it falls on us to avoid that in SeisIO. (Julia is much better about this than R right now, but that's not setting the bar high enough; R's package ecosystem is a superfund site. Python machine learning might be even worse...)

@dylanmikesell
Copy link
Author

Any thoughts on which library? Interpolations.jl seemed pretty good to me after looking through the documentation. Last commit was 15 days ago so it looks like they are still maintaining this package.

@tclements
Copy link
Contributor

Just to chime in - I've used Interpolations.jl and it's had everything I've needed. Dierckx.jl wraps the same Fortran library that Scipy uses for interpolation. Interpolations.jl is probably easier to debug and the author is one of the most active contributors to Julia.

@dylanmikesell
Copy link
Author

I have a working version of the interpolation I need using Interpolations.jl. It is actually very slick in my opinion. I am going to keep using this package unless I hear otherwise.

@jpjones76
Copy link
Owner

Aha, yes, I think Interpolations.jl is best. Sorry for not responding sooner. I swear I replied to this yesterday, but there's no sign of my comment...??

@dylanmikesell
Copy link
Author

No worries. I got the basics finished today. See image that compares new Julia code to Matt's MATLAB output. This is an example that compares colocated broadband and short period sensors in the AV network. Station name is CLES, and this is what Matt gives an example data set when he shares the MATLAB code. You can see that in Julia we have exactly the same results. The BB and SP signals match from MATLAB or Julia, and the small discrepancies between the BHZ and SHZ channels are consistent in MATLAB or Julia.
decon_compare_zoom

I need to clean up a few more things and then I will share with @jpjones76 so you can help me figure out where to integrate into the SeisIO.jl base. I will start with a new file in src/Processing/Resp/.

@jpjones76
Copy link
Owner

Hi, I haven't heard back from you about this since my email of August 12. Haney's methods are elegant but prohibitively memory-intensive; can you please provide an example to demonstrate the need? I'm looking specifically for something where your code works, but ordinary instrument response translation fails: for example, if you can show me a case where detrend! ==> taper! ==> translate_resp! adds acausal artifacts that bias subsequent analysis, that's a pretty compelling argument.

@dylanmikesell
Copy link
Author

dylanmikesell commented Dec 2, 2020 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants