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

Limit number of fft bins.

Browse Source
This commit is contained in:
Ludger Küperkoch 2020-09-21 17:01:06 +02:00
parent f93499da7d
commit 37c8858096
1 changed files with 2 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -522,7 +522,6 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
Fc = None Fc = None
w0 = None w0 = None
zdat = select_for_phase(wfstream, "P") zdat = select_for_phase(wfstream, "P")
if len(zdat) == 0: if len(zdat) == 0:
@ -537,7 +536,6 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
# trim traces to common range (for rotation) # trim traces to common range (for rotation)
trstart, trend = common_range(wfstream) trstart, trend = common_range(wfstream)
wfstream.trim(trstart, trend) wfstream.trim(trstart, trend)
# rotate into LQT (ray-coordindate-) system using Obspy's rotate # rotate into LQT (ray-coordindate-) system using Obspy's rotate
# L: P-wave direction # L: P-wave direction
# Q: SV-wave direction # Q: SV-wave direction
@ -591,7 +589,8 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
#n = freq * l #n = freq * l
# find next power of 2 of data length # find next power of 2 of data length
m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2))) m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2)))
N = int(np.power(m, 2)) N = min(int(np.power(m, 2)), 16384)
#N = int(np.power(m, 2))
y = dt * np.fft.fft(xdat, N) y = dt * np.fft.fft(xdat, N)
Y = abs(y[: N / 2]) Y = abs(y[: N / 2])
L = (N - 1) / freq L = (N - 1) / freq