My apology if this is the wrong place to send this.
I noticed a potential artifact in using the PPSD function at very long periods. It appears that when you deconvolve the data, convert the data to acceleration and then estimate the psd you might be elevating the long period spectra.
Attached is an example of using the PPSD function from obspy.spectral_estimation (example.png). I have also included an example of first estimating the power spectra and then converting to acceleration by dividing out the response and using a factor of omega to switch to acceleration.
As you can see the second plot rides a few dB above the NLNM while (which I would expect) the first example shows a “hump” with elevated noise levels past a couple 100 seconds period (which I would not expect).
Below are more details.
Best,
Adam
I added the following code to a local version of spectral_estimation.py (line 663) to make this change
looking at the response
(http://service.iris.edu/irisws/evalresp/1/query?net=IU&sta=TUC&loc=00&cha=LHZ&time=2013-04-21T10:00:00&output=plot),
I suspect that the "Trace.simulate()" call inside PPSD is responsible.
It looks like it's using obspy.signal.invsim.seisSim's default values in
the PAZ/overall sensitivity/waterlevel deconvolution. seisSim is one of
signal's oldest routines and back then a default value of 600 for the
waterlevel got implemented, which is much too high from my experience
and might cause noise to blow up in your example. I'm using 10-30 usually.
Maybe you can try to check if hardcoding a more conservative value in
spectral_estimation
(https://github.com/obspy/obspy/blob/master/obspy/signal/spectral_estimation.py#L640-L641)
helps, or you could try to exchange that call with the "seedresp"-option
in simulate, providing a RESP file or Parser object (from dateless) and
going via evalresp (specifying a pre filtering).
Thanks for the help. I just confirmed the elevated vlp noise is coming from the water level. I switched it to 30 and it looked slightly elevated (however a definite improvement) and at 10 it went below the NLNM at about 1000 seconds period. To me this suggests it is probably better to remove the response from the power spectra and not on the time series for estimating pdfs at very long periods. Using paztofreqresp seems to solve this issue.
good to hear that it's cleared up.
If you think you have a way of processing that's better for the non-flat
parts or maybe if you even can help improve testing / cross checking the
PPSD routine, please don't hesitate to open an issue or pull request on
github.
