adding spectral methods to trace object

[megies]

Coming from the 'old school' interactive side, I dearly miss spectral methods attached directly to a trace object. While the spectrogram, which to me is a more sophisticated method is already available as a method, the direct estimation of psd/multi-taper and fft is not available. Certainly this can be done by modules of scipy or mlab BUT what I'm thinking of is to return a complex trace (resulting from a fft), manipulate it and possibly transfer it back to time domain. This means that all trace.stats information should be kept while the sampling or delta needs to be modified to frequency domain (i.e. 1/(time window). This would open up the world of trace manipulation in frequency domain, easy cross-spectral estimates, sub-sample shifts etc AND of course an easy and fast way of displaying the frequency content. Possibly this would result in a new 'trace object' and therefore needs a re-definition of the Trace object itself.
What is the opinion of the developers?

Cheers Jo

[claudiodsf]

👍 on this. I developed a code for source parameters estimation from spectral inversion, and I needed to build a custom spectral object to deal with spectra.

I don't feel like sharing the code here, since it is an horrible hack created from subclassing a Trace object.

I think that a proper implementation of spectral methods would benefit a significant class of applications.

[megies]

Talked with @jwassermann a bit about it.. while stats header entries would be pretty similar to a time domain Trace, it would probably be good to introduce a derived "frequency domain" Trace object. That way we could raise errors on e.g. mixed multiplications of "time domain Trace" and "frequency domain Trace". Furthermore we could prevent operations on one or the other that only make sense for exactly one type..

I'd propose adding a class BaseTrace and move everything from the current Trace class that applies to both types there. Then derive Trace and a new FrequencyTrace (or some other name) from BaseTrace.

Comments? There's a lot to be thought through properly here...

[krischer]

Good idea. The BaseTrace class will likely be pretty small as not much functionality can be shared but that is fine.

Just some ideas on how this could work:

import obspy

tr = obspy.read()[0]

# Go to frequency domain, maybe with other methods/options as well.
f_tr = tr.fft()

# data and freqs are stored here.
f_tr.data
f_tr.freqs

# Frequency domain taper, e.g. some filter.
f_tr.taper(...)
f_tr.plot()
...

# Go back to time domain. For that to be possible the starttime
# must also be stored for the frequeny domain trace.
tr = f_tr.ifft()

[claudiodsf]

f_tr.data would be a complex number.
Should we also store the same information as phase and amplitude? Or maybe provide methods to get them easily.

f_tr.freqs should probably be a function, like tr.times()

[claudiodsf]

Other thoughts: I would like to have trim() and slice() methods for spectral traces, as well (I'm using them in my custom implementation for the spectral inversion code).

What about streams? Could a Stream object mix Trace and FrequencyTrace objects? Should there be a FrequencyStream?

[krischer]

There should be FrequencyStream class in my opinion. I would also think that f_tr.data should be a complex number. We could provide f_tr.phase(), and f_tr.amplitude() or some related functions to calculate these things but the actual data value should be complex - many people will want to have access to that.

I agree with f_tr.freqs() being a function.

I don't really get the trim() and slice() operations in the frequency domain. Do you want to cut away frequency bins? Or have the same interpretation as in the time domain, e.g. time shift in the frequency domain and somehow reinterpolate/throw away some frequencies?

[claudiodsf]

trim() or slice() --or even spectral re-interpolation-- would be useful for people needing to stack spectra or to fit a spectral model using just a part of the spectrum (e.g., ignoring too low or too high frequencies).

[megies]

There should be FrequencyStream class in my opinion.

👍 we can't allow to mix both types.

I would also think that f_tr.data should be a complex number. We could provide f_tr.phase(), and f_tr.amplitude() or some related functions to calculate these things but the actual data value should be complex - many people will want to have access to that.

👍 for storing data as complex. talked with @jwassermann about it, we will need to define properties f_tr.amplitude, f_tr.phase that

• have getter methods that compute from complex data on the fly
• have setter methods that update the complex data (first getting the complementary other real data from the stored complex data)

I agree with f_tr.freqs() being a function.

absolutely. we should probably store this in stats like with `starttime/sampling_rate/npts

[QuLogic]

I'd propose adding a class BaseTrace

That proposal appears to coincide with some thoughts on #617, I think?

[megies]

That proposal appears to coincide with some thoughts on #617, I think?

Indeed, although I think there were more people speaking out against #617.

[petrrr]

Hi guys,
I was not aware of this discussion, but by incidence I proposed something quite similar earlier today (what a coincidence) in #841.

[megies]

2 I started to build up the frequency_domain_trace.py, but got problems with the circular imports in python, because I generate in the tr.fft Method an instance of the FrequencyDomainTrace class and in the tr_f.ifft Method an instance of the Trace class. Is there a nice way to do it?

If it's gonna be in separate files as @krischer proposed, the only way is to do the from obspy import Trace import not at top of frequency_domain_trace.py but inside the def ifft(self, ...): block (and vice versa with import of FrequencyDomainTrace in trace.py).

[fstettner]

I moved the FrequencyDomainTrace to the frequency_domain_trace.py like @krischer recommended. Additionally, I built in the numpy rfft, irfft and fftfreq , because of the return type of the scipy version. The fft is now attached to the TimeSeriesTrace with an differentiation of Trace and TimeSeriesTrace. I also started to build up the TimeSeriesFrequencyDomainTrace, that has currently just the ifft implemented.

[megies]

@fstettner, smoothing of spectra is missing here..

In the long run I would propose adding a method to FFTTrace, e.g. def smoothed(self, window, npts) (or width given in Hz).
For now you could just hardcode some fixed window from scipy, see http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window. For orientation, Poupinet 1984 used a five point triangular window for 100Hz sampled data.

CC @jwassermann

[fstettner]

I implemented the conjugate method. The coherence-method also should work now (with hanning).
I still have problems with the normalization of the cross_correlation-method. Are there some ideas?
Are the differentiation and integration method useful implemented?

[krischer]

This should just use the self.frequencies() method which would me much easier.

[ThomasLecocq]

👍 to get this done asap!

[megies]

This was done as a student project and likely needs more work.. :-/

[obspy-bot]

This pull request has been mentioned on ObsPy Forum. There might be relevant details there:

https://discourse.obspy.org/t/trace-remove-response-outputting-disp-or-acc/1259/6