marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

In order to calculate DC value and corner frequency of source spectrum a synthetic spectrum is calculated and optimized using scipys curve_fit.

Browse Source
This commit is contained in:
Ludger Küperkoch 2015-09-29 11:17:47 +02:00
parent 8035903fa5
commit ce57f184e7
1 changed files with 60 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
pylot/core/analysis/magnitude.py
View File

@ -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