From 8035903fa5025955c822595f2c67423eaa93a77a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 29 Sep 2015 11:15:51 +0200 Subject: [PATCH 1/4] Zero crossings are calculated to derive only P pulse for calculating source spectrum. --- pylot/core/pick/autopick.py | 43 +++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 10b4f7d5..e3771d7a 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -313,7 +313,6 @@ def autopickstation(wfstream, pickparam): ############################################################## # get DC value (w0) and corner frequency (fc) of source spectrum # from P pulse - # restitute streams # initialize Data object data = Data() [corzdat, restflag] = data.restituteWFData(invdir, zdat) @@ -323,33 +322,45 @@ def autopickstation(wfstream, pickparam): # class needs stream object => build it z_copy = zdat.copy() z_copy[0].data = corintzdat - # calculate source spectrum and get w0 and fc - calcwin = 1 / bpz2[0] # largest detectable period == window length - # around P pulse for calculating source spectrum - specpara = DCfc(z_copy, mpickP, calcwin, iplot) - w0 = specpara.getw0() - fc = specpara.getfc() + # largest detectable period == window length + # after P pulse for calculating source spectrum + winzc = (1 / bpz2[0]) * z_copy[0].stats.sampling_rate + impickP = mpickP * z_copy[0].stats.sampling_rate + wfzc = z_copy[0].data[impickP : impickP + winzc] + # calculate spectrum using only first cycles of + # waveform after P onset! + zc = crossings_nonzero_all(wfzc) + if np.size(zc) == 0: + print ("Something is wrong with the waveform, " \ + "no zero crossings derived!") + print ("Cannot calculate source spectrum!") + else: + calcwin = (zc[3] - zc[0]) * z_copy[0].stats.delta + # calculate source spectrum and get w0 and fc + specpara = DCfc(z_copy, mpickP, calcwin, iplot) + w0 = specpara.getw0() + fc = specpara.getfc() - print 'autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, ' \ - 'Polarity: %s' % (Pweight, SNRP, SNRPdB, FM) + print ("autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, " \ + "Polarity: %s" % (Pweight, SNRP, SNRPdB, FM)) Sflag = 1 else: - print 'Bad initial (AIC) P-pick, skipping this onset!' + print ("Bad initial (AIC) P-pick, skipping this onset!") print 'AIC-SNR=', aicpick.getSNR(), 'AIC-Slope=', aicpick.getSlope(), 'counts/s' print '(min. AIC-SNR=', minAICPSNR, ', min. AIC-Slope=', minAICPslope, 'counts/s)' Sflag = 0 else: - print 'autopickstation: No vertical component data available!, ' \ - 'Skipping station!' + print ("autopickstation: No vertical component data available!, " \ + "Skipping station!") if edat is not None and ndat is not None and len(edat) > 0 and len( ndat) > 0 and Pweight < 4: - print 'Go on picking S onset ...' - print '##################################################' - print 'Working on S onset of station %s' % edat[0].stats.station - print 'Filtering horizontal traces ...' + print ("Go on picking S onset ...") + print ("##################################################") + print ("Working on S onset of station %s" % edat[0].stats.station) + print ("Filtering horizontal traces ...") # determine time window for calculating CF after P onset cuttimesh = [round(max([mpickP + sstart, 0])), From ce57f184e755201de5f480a879b9d051755a91bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 29 Sep 2015 11:17:47 +0200 Subject: [PATCH 2/4] In order to calculate DC value and corner frequency of source spectrum a synthetic spectrum is calculated and optimized using scipys curve_fit. --- pylot/core/analysis/magnitude.py | 66 +++++++++++++++++++++++++++++--- 1 file changed, 60 insertions(+), 6 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 4391d16b..17d3c16f 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -9,6 +9,7 @@ import matplotlib.pyplot as plt import numpy as np from obspy.core import Stream from pylot.core.pick.utils import getsignalwin +from scipy.optimize import curve_fit class Magnitude(object): ''' @@ -165,21 +166,74 @@ class DCfc(Magnitude): Y = abs(y[: N/2]) L = (N - 1) / tr.stats.sampling_rate f = np.arange(0, fny, 1/L) + + # remove zero-frequency and frequencies above + # corner frequency of seismometer (assumed + # to be 100 Hz) + fi = np.where((f >= 1) & (f < 100)) + F = f[fi] + YY = Y[fi] + # get plateau (DC value) and corner frequency + # initial guess of plateau + DCin = np.mean(YY[0:100]) + # initial guess of corner frequency + # where spectral level reached 50% of flat level + iin = np.where(YY >= 0.5 * DCin) + Fcin = F[iin[0][np.size(iin) - 1]] + fit = synthsourcespec(F, DCin, Fcin) + [optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [DCin, Fcin]) + self.w0 = optspecfit[0] + self.fc = optspecfit[1] + print ("DCfc: Determined DC-value: %f, \n" \ + "Determined corner frequency: %f" % (self.w0, self.fc)) - if self.getiplot() > 1: - f1 = plt.figure(1) + + #if self.getiplot() > 1: + iplot=2 + if iplot > 1: + f1 = plt.figure() plt.subplot(2,1,1) - plt.plot(t, np.multiply(tr, 1000), 'k') # show displacement in mm - plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') # show displacement in mm + # show displacement in mm + plt.plot(t, np.multiply(tr, 1000), 'k') + plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') plt.title('Seismogram and P pulse, station %s' % tr.stats.station) plt.xlabel('Time since %s' % tr.stats.starttime) plt.ylabel('Displacement [mm]') plt.subplot(2,1,2) - plt.semilogy(f, Y.real) - plt.title('Source Spectrum from P Pulse') + plt.semilogy(f, Y.real, 'k') + plt.semilogy(F, YY.real) + plt.semilogy(F, fit, 'g') + plt.title('Source Spectrum from P Pulse, DC=%f m/Hz, fc=%4.1f Hz' \ + % (self.w0, self.fc)) plt.xlabel('Frequency [Hz]') plt.ylabel('Amplitude [m/Hz]') plt.show() raw_input() plt.close(f1) + + +def synthsourcespec(f, omega0, fcorner): + ''' + Calculates synthetic source spectrum from given plateau and corner + frequency assuming Akis omega-square model. + + :param: f, frequencies + :type: array + + :param: omega0, DC-value (plateau) of source spectrum + :type: float + + :param: fcorner, corner frequency of source spectrum + :type: float + ''' + + #ssp = omega0 / (pow(2, (1 + f / fcorner))) + ssp = omega0 / (1 + pow(2, (f / fcorner))) + + #plt.plot(f, ssp) + #plt.show() + #raw_input() + + return ssp + From f4b905c2e63c3f8bec92a97f5050dbd899d25766 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 29 Sep 2015 11:19:25 +0200 Subject: [PATCH 3/4] Removed inserted plot command for debugging purposes. --- pylot/core/analysis/magnitude.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 17d3c16f..4fa02704 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -188,8 +188,7 @@ class DCfc(Magnitude): "Determined corner frequency: %f" % (self.w0, self.fc)) - #if self.getiplot() > 1: - iplot=2 + if self.getiplot() > 1: if iplot > 1: f1 = plt.figure() plt.subplot(2,1,1) From 708f0a1f1a096551a61316f0c51baacec859d33d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 29 Sep 2015 14:55:15 +0200 Subject: [PATCH 4/4] Some cosmetics on DCfc of class magnitude. --- pylot/core/analysis/magnitude.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 4fa02704..9d24392e 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -184,12 +184,11 @@ class DCfc(Magnitude): [optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [DCin, Fcin]) self.w0 = optspecfit[0] self.fc = optspecfit[1] - print ("DCfc: Determined DC-value: %f, \n" \ - "Determined corner frequency: %f" % (self.w0, self.fc)) + print ("DCfc: Determined DC-value: %e m/Hz, \n" \ + "Determined corner frequency: %f Hz" % (self.w0, self.fc)) if self.getiplot() > 1: - if iplot > 1: f1 = plt.figure() plt.subplot(2,1,1) # show displacement in mm @@ -200,13 +199,14 @@ class DCfc(Magnitude): plt.ylabel('Displacement [mm]') plt.subplot(2,1,2) - plt.semilogy(f, Y.real, 'k') - plt.semilogy(F, YY.real) - plt.semilogy(F, fit, 'g') - plt.title('Source Spectrum from P Pulse, DC=%f m/Hz, fc=%4.1f Hz' \ + plt.loglog(f, Y.real, 'k') + plt.loglog(F, YY.real) + plt.loglog(F, fit, 'g') + plt.title('Source Spectrum from P Pulse, DC=%e m/Hz, fc=%4.1f Hz' \ % (self.w0, self.fc)) plt.xlabel('Frequency [Hz]') plt.ylabel('Amplitude [m/Hz]') + plt.grid() plt.show() raw_input() plt.close(f1) @@ -230,9 +230,5 @@ def synthsourcespec(f, omega0, fcorner): #ssp = omega0 / (pow(2, (1 + f / fcorner))) ssp = omega0 / (1 + pow(2, (f / fcorner))) - #plt.plot(f, ssp) - #plt.show() - #raw_input() - return ssp