[obspy-users] spectrum estimation - e.g., welch
megies at geophysik.uni-muenchen.de
Thu Jul 2 11:51:14 CEST 2015
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:
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.
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,
On 07/01/2015 09:51 PM, Mike Hagerty wrote:
> 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.
> for trace in stream: # Here stream is a multi-station mseed file
> vals, freq = psd(trace.data, NFFT=256, Fs=trace.stats.delta,
> detrend=mlab.detrend_mean, \
> noverlap=0, pad_to=None, sides='default', scale_by_freq=None)
> On 7/1/15 12:33 PM, Ringler, Adam wrote:
>> Hey Mike,
>> Does the following help you at all?
>> So something like the following (this is for the cross-power instead
>> of the power)?
>> On Wed, Jul 1, 2015 at 10:16 AM, Mike Hagerty <m.hagerty at isti.com
>> <mailto:m.hagerty at isti.com>> wrote:
>> 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 ?
>> -Mike Hagerty
>> obspy-users mailing list
>> obspy-users at lists.swapbytes.de <mailto:obspy-users at lists.swapbytes.de>
>> Follow USGS Seismic Instrumentation on Twitter: @USGS_Seismic
>> obspy-users mailing list
>> obspy-users at lists.swapbytes.de
> obspy-users mailing list
> obspy-users at lists.swapbytes.de
Dipl.-Geophys. Tobias Megies
Department für Geo- und Umweltwissenschaften
Tel: +49 (0) 89 2180-73981
+49 (0) 89 2180-4326
Mail: tobias.megies at geophysik.uni-muenchen.de
More information about the obspy-users