detrending by fitting a linear function

Question

Before adding this to obspy, I decided put it here for discussion, maybe there is some problem I did not anticipate:

For now, detrending is done using the first and last and last value of the trace in obspy.signal.invsim.detrend. This fails with a lot of noise or when the time window of the trace contains only part of the signal. Wouldn't it be better to fit a linear function to the trace and remove that from the data? E.g. like this:

{{{ for tr in st: dt = tr.stats.delta t = np.linspace(0., (tr.stats.npts - 1) * dt, tr.stats.npts)

A = np.vstack([t, np.ones(len(t))]).T
m, c = np.linalg.lstsq(A, tr.data)[0]

tr.data = tr.data - m*t - c

}}}

Works fine for me.

Cheers, Martin

该提问来源于开源项目：obspy/obspy

Answer 1

[barsch] I can't answer from a scientific view as I'm not really into in signal processing. However as far as I know its always going down to the question what you actually want to to with your data. Some prefer method x to do something while some prefer something else.

Therefore I would suggest to include "your" detrending function as an addition into obspy.signal - either by introducing a new keyword "method=..."to the current detrend function or giving it a completely new name - pretty much the same we did with xcorr, merge etc.

just my 2 cents Robert

Answer 2

[barsch] oh most important part is missing in the last paragraph above: "as long you have a proper documentation for the new method it in the docstring" ;)

Answer 3

[driel] Sure, would keep the old version for compatibility reasons. Also, for very long data, the fitting might be slow and memory intensive.

Answer 4

[anonymous] Accidently answered with an email that does not show up here... So once again:

I have always used scipy's detrend routine which, like your solution, just subtracts a linear least square fit:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html

It also has on optional 'method' argument so one could also just remove the mean value.

Maybe just wrap that routine? I also agree that simply removing a line through the first and last sample makes little sense.

Answer 5

[megies] I think using the scipy routine is the way to go. I would suggest implementing a detrend() method on Trace/Stream to make this easier to handle for users. We can leave the current detrend function in obspy.signal.invsim untouched and use the scipy routine directly at Trace/Stream. probably also with an kwarg switch like: {{{ trace.detrend(method="LS")

method == "LS": use scipy least squares
method == "simple": use method from `signal.invsim` with first/last sample

}}} With the least squares linear regression approach it might be a good idea to also include tapering (controlled with another kwarg maybe) for a few samples at the start. Otherwise it is not guaranteed that the data starts at 0 which is important for any filtering coming after detrending.

Answer 6

[driel] I didn't know the function proposed by Lion, thanks for that. I agree with Tobias as for adding the detrend to Stream / Trace and I suggest when someone starts with wrapping functions in Stream and Trace, it wouldn't be much more then copy and paste to also do sth. about #265 (also using a method argument and using simple diff/cumptrapz as a start).

I'll have a look at it, but might need some help with the tests at some point.

Answer 7

[krischer] I would also suggest to add some routines to Trace/Stream. I think its cleaner to add seperate detrend and taper convenience methods and maybe even give the optional possibility to detrend and taper right in the filter method, e.g. {{{ stream.filter('lowpass', freq=1.0, corners=2, zerophase=True, detrend='LS', taper='hanning') }}}

Answer 8

[driel] In [2967/obspy]: {{{

!CommitTicketReference repository="obspy" revision="2967"

• adding detrend to obspy.core.trace
• see #304 }}}

Answer 9

[driel] In [2968/obspy]: {{{

!CommitTicketReference repository="obspy" revision="2968"

• addin detrend to the stream
• see #304
• can someone please check wether I sticked to all conventions and if tests are running? }}}

Answer 10

[driel] In [2969/obspy]: {{{

!CommitTicketReference repository="obspy" revision="2969"

• adding convinience wrapping of differentiation to Stream / Trace
• keyword method included for future implementation of fourier domain differentiation
• see #265 and #304 }}}

Answer 11

[driel] In [2970/obspy]: {{{

!CommitTicketReference repository="obspy" revision="2970"

• adding convinience wrapping for integration (cumtrapz for now) to stream and trace
• see #304 and #265 }}}

Answer 12

[driel] Okay, I implemented part of what we discussed above. I would ask one of the core developers to quickly double check. How can I check the docstrings format? docs.obspy.org is not updated after svn commit I guess...?

Answer 13

[driel] Right, thanks Robert for cleaning up. The only thing I dislike is the detrend(type='constant'), I'd rather suggest to call it demean(). As it does not remove any 'trend', the name detrend is not intuitive to me.

Answer 14

[barsch] Its taken from http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html - change it if you don't like it or just add another keyword 'demean' doing exactly the same thing ;)

Answer 15

[megies] During filtering, data arrays get converted to float.. should we do this for detrend ({{{tr.data = tr.data.astype('float')}}} before detrending), too, what do you think?

Answer 16

[driel] I don't see any reason not to, but a couple of good ones to do so...

Answer 17

[krischer] They definitely have to be converted to floats. The detrending array should be a float array anyway and therefore the resulting detrended data array should also be float.

Assuming tr.data is an integer array and detrend a float array, then

{{{ tr.data -= detrend }}} results in an integer array which is wrong. {{{ tr.data = tr.data - detrend }}} on the other hand would result in a float array but I would prefer the following: {{{ tr.data = np.require(tr.data, 'float32') tr.data -= detrend }}}

Answer 18

[megies] .astype()/require() is the only option since the rest is done inside the scipy routine that is called (and no type conversion done there).

Answer 19

[driel] I kind of lost track, what's the status here? Can this ticket be closed or is there some issue with the float stuff discussed above?