detrending by fitting a linear function

[driel]

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

[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

[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" ;)

[driel]

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

[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.

[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.

[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.

[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')

[driel]

In [2967/obspy]:

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

[driel]

In [2968/obspy]:

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

[driel]

In [2969/obspy]:

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

[driel]

In [2970/obspy]:

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

[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...?

[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.

[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 ;)

[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?

[driel]

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

[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

[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).

[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?