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