search

ThomasLecocq / SeismoRMS

A simple Jupyter Notebook example for getting the RMS of a seismic signal (from PSDs)
European Union Public License 1.1
87 stars 38 forks source link

visual proof of RMS(PSD) == RMS(time domain) #7

Open ThomasLecocq opened 4 years ago

ThomasLecocq commented 4 years ago

@claudiodsf @FMassin @megies

should I include this ipynb that shows the correctness of the computations: ?? I mean we know it's Parsaval's identity, but still:

http://msnoise.org/tmp/Amplitude%20&%20Spectrum.html

And shows advantages:

  1. computing the PSD doesn't require pre_filtering the whole trace
  2. PSD doens't need to remove_response (ok this could be done by just tr.data / total_sensitivity, right?)
  3. RMS derived from PSD (calculated once) can be extracted for any frequency band, even wide/narrow

PS: in the example above, you'll see a ~33% noise reduction when it snows in Brussels (2 to 5cm max) and all the trafic is stopped.. Kinda matches what we see today :-)

megies commented 4 years ago

PSD doens't need to remove_response (ok this could be done by just tr.data / total_sensitivity, right?)

Only for the flat part of the response, but then yes. RaspiShake has ~5Hz geophones I think, though, so already there its not 100% valid if you look at 1-5Hz signal. Phase shifts also appear in the flat part but that should not matter too much in this case I guess

should I include this ipynb that shows the correctness of the computations: ?? I mean we know it's Parsaval's identity, but still:

I mean, can't hurt, I'd say

FMassin commented 4 years ago
  1. PSD doens't need to remove_response (ok this could be done by just tr.data / total_sensitivity, right?)

Well I believe this is taken care of (in frequency domain) by the PPSD init and add method (see https://github.com/obspy/obspy/blob/88687a146a7c3ca8f35c608db2243cf5fed6813c/obspy/signal/spectral_estimation.py#L973). It matters a lot so we can compare multiple stations.

ThomasLecocq commented 4 years ago
  1. PSD doens't need to remove_response (ok this could be done by just tr.data / total_sensitivity, right?)

Well I believe this is taken care of (in frequency domain) by the PPSD init and add method (see https://github.com/obspy/obspy/blob/88687a146a7c3ca8f35c608db2243cf5fed6813c/obspy/signal/spectral_estimation.py#L973). It matters a lot so we can compare multiple stations.

yes yes, sorry, thats what I meant, but it's not needed on the full 86400*sps samples like when doing the raw trace

ThomasLecocq commented 4 years ago

ok , will add some markdown to explain the different steps & why, later today