spectrum estimation - e.g., welch

Hello.

I'm trying to do a quick spectrum estimation in python/obspy.

My choices seem to be:

1. obspy.signal.spectral_estimation.psd() = wrapper for matplotlib.mlab.psd() - no documentation on return type(s),
     either on obspy or matplotlib pages. Seems to return a matplotlib GUI plot of psd by default (how to change this ?)

2. obspy.signal.freqattributes.welch(data, win, Nfft, L=0, over=0)
     data = numpy.ndarray
     win = ? numpy.ndarray ?

I'm doing:
     spec = welch(data, welch_window(NFFT), NFFT)

     where len(data) = multiple of NFFT

e.g., if NFFT=256 and len(data)=512, welch gives:
     Caught e= operands could not be broadcast together with shapes (512,) (256,)

     So clearly, I am not dimensioning this right somehow (?)
     Are there any working examples of this module ?

     My understanding is that welch should be dividing the data up into segments of length=NFFT, multiply each of
     these with the window, compute the periodogram, and average the segments together.

     Also, the obspy doc for welch() says that the window is optional, but I don't see how to use this feature.
     e.g., welch(data, 1, NFFT) ?? How do I call it without a window array ?

Cheers!
-Mike Hagerty

Hey Mike,

Does the following help you at all?

http://matplotlib.org/api/mlab_api.html#matplotlib.mlab.psd

So something like the following (this is for the cross-power instead of the power)?

cpval,fre = csd(tr1.data,tr2.data,NFFT=lenfft,Fs=sr,noverlap=lenol,scale_by_freq=True)

Hi Adam,

your syntax is similar to mine.

I think there is an issue with how obspy wraps mlab.psd.

In the code snippet below, if I comment out “trace.plot()”, then it loops through the stations
as expected (no extra psd plots).
If I uncomment trace.plot(), then the first plot I see is the trace plot of the first channel as expected.
However, just before I see the trace plot of the second channel, psd throws up ‘Fig 1’ of the
first channel’s psd.

So, either obspy is generating the plots (?) or is passing them forward from mlab by mistake (?)

I don’t see anything in the mlab docs that suggest it would autogenerate the psd plots, so I think
maybe obspy is doing this.

Cheers!
-Mike

Hi Mike,

I just got forwarded some code pieces related to psd and plotting. I
didn't really review it myself but dropped it onto github for
visibility, could be it's interesting/helpful for you.
check it out in this not yet polished pull request:
https://github.com/obspy/obspy/pull/1088/files

Also the psd wrapper in obspy will get deprecated in the next major
release as it's not really necessary anymore with ancient matplotlib
versions dying out, so you might consider using mlabs psd directly.
https://github.com/obspy/obspy/commit/687705f5912854e9dd341b7e66cd3a86a063bd1a

The opening of (empty?) figures by psd is unwanted I would think,
although we'd need some self contained example code to reproduce. Note
the difference of using matplotlib.pyplot.psd vs. matplotlib.mlab.psd

hope it helps,
Tobias