diff --git a/.github/workflows/testbuild.yml b/.github/workflows/testbuild.yml index 6cfb2b8f..3cdfb94d 100644 --- a/.github/workflows/testbuild.yml +++ b/.github/workflows/testbuild.yml @@ -274,6 +274,13 @@ jobs: # python plot_compare_pygrt.py # continue-on-error: true # 即使失败，仍然标记为成功 + - name: Test signal convolution + working-directory: ${{ env.EXAMPLE_COPY_PATH }}/convolve_signals + run: | + chmod +x *.sh + ./run_grt.sh + python run_pygrt.py + - name: Remove the test files run: | rm -rf ${{ env.EXAMPLE_COPY_PATH }} diff --git a/example/convolve_signals/milrow b/example/convolve_signals/milrow new file mode 100644 index 00000000..a8a74b18 --- /dev/null +++ b/example/convolve_signals/milrow @@ -0,0 +1,9 @@ +0.2 3.4 1.7 2.3 9e10 9e10 +0.6 3.7 1.9 2.4 9e10 9e10 +0.5 4.2 2.1 2.4 9e10 9e10 +0.5 4.6 2.3 2.5 9e10 9e10 +0.7 4.9 2.8 2.6 9e10 9e10 +0.5 5.1 2.9 2.7 9e10 9e10 +6.0 5.9 3.3 2.7 9e10 9e10 +28. 6.9 4.0 2.8 9e10 9e10 +-0.1 8.2 4.7 3.2 9e10 9e10 diff --git a/example/convolve_signals/plot.sh b/example/convolve_signals/plot.sh new file mode 100755 index 00000000..0f5aa923 --- /dev/null +++ b/example/convolve_signals/plot.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +SAC_DISPLAY_COPYRIGHT=0 + +out="compare_pdf" +rm -rf $out +mkdir -p $out + + +for synpath in syn*; do +for cout in $synpath/cout*; do +fname=$(basename $cout) +echo $cout $(dirname $cout)/p${fname:1:} +sac << EOF +r $cout $(dirname $cout)/p${fname:1} +qdp off +color red inc list blue red +p1 +saveimg $out/compare_${synpath}_${fname}.pdf +q +EOF + +done +done + + +pdftk $out/* cat output out.pdf \ No newline at end of file diff --git a/example/convolve_signals/run_grt.sh b/example/convolve_signals/run_grt.sh new file mode 100755 index 00000000..8f463171 --- /dev/null +++ b/example/convolve_signals/run_grt.sh @@ -0,0 +1,40 @@ +#!/bin/bash + + +dist=10 +depsrc=2 +deprcv=0.5 + +nt=1024 +dt=0.01 + +modname="milrow" +out="GRN" + +grt -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${dist} + +# convolve different signals +G=$out/${modname}_${depsrc}_${deprcv}_${dist} +P=cout +S=1e24 +az=39.2 +grt.syn -G$G -Osyn_exp -P$P -A$az -S$S -Dt/0.2/0.2/0.4 + +fn=2 +fe=-1 +fz=4 +grt.syn -G$G -Osyn_sf -P$P -A$az -S$S -F$fn/$fe/$fz -Dt/0.1/0.3/0.6 + +stk=77 +dip=88 +rak=99 +grt.syn -G$G -Osyn_dc -P$P -A$az -S$S -M$stk/$dip/$rak -Dp/0.6 + +M11=1 +M12=-2 +M13=-5 +M22=0.5 +M23=3 +M33=1.2 +grt.syn -G$G -Osyn_mt -P$P -A$az -S$S -T$M11/$M12/$M13/$M22/$M23/$M33 -Dr/3 + diff --git a/example/convolve_signals/run_pygrt.py b/example/convolve_signals/run_pygrt.py new file mode 100644 index 00000000..1ef92f38 --- /dev/null +++ b/example/convolve_signals/run_pygrt.py @@ -0,0 +1,41 @@ +import numpy as np +import pygrt + +dist=10 +depsrc=2 +deprcv=0.5 + +nt=1024 +dt=0.01 + +modname="milrow" + +modarr = np.loadtxt(modname) + +pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) + +# compute green functions +st_grn = pymod.gen_gf_spectra(dist, nt, dt)[0] + +# synthetic +S=1e24 +az=39.2 +st = pygrt.utils.gen_syn_from_gf_EXP(st_grn, S, az) +sigs = pygrt.sigs.gen_triangle_wave(0.4, dt) +pygrt.utils.stream_convolve(st, sigs) +pygrt.utils.stream_write_sac(st, "syn_exp/pout") + +st = pygrt.utils.gen_syn_from_gf_SF(st_grn, S, 2, -1, 4, az) +sigs = pygrt.sigs.gen_trap_wave(0.1, 0.3, 0.6, dt) +pygrt.utils.stream_convolve(st, sigs) +pygrt.utils.stream_write_sac(st, "syn_sf/pout") + +st = pygrt.utils.gen_syn_from_gf_DC(st_grn, S, 77, 88, 99, az) +sigs = pygrt.sigs.gen_parabola_wave(0.6, dt) +pygrt.utils.stream_convolve(st, sigs) +pygrt.utils.stream_write_sac(st, "syn_dc/pout") + +st = pygrt.utils.gen_syn_from_gf_MT(st_grn, S, [1,-2,-5,0.5,3,1.2], az) +sigs = pygrt.sigs.gen_ricker_wave(3, dt) +pygrt.utils.stream_convolve(st, sigs) +pygrt.utils.stream_write_sac(st, "syn_mt/pout") \ No newline at end of file diff --git a/pygrt/C_extension/include/dynamic/signals.h b/pygrt/C_extension/include/dynamic/signals.h index ec512456..4f2ae649 100644 --- a/pygrt/C_extension/include/dynamic/signals.h +++ b/pygrt/C_extension/include/dynamic/signals.h @@ -167,4 +167,11 @@ float * get_ricker_wave(float dt, float f0, int *Nt); * * @return float指针 */ -float * get_custom_wave(int *Nt, const char *tfparams); \ No newline at end of file +float * get_custom_wave(int *Nt, const char *tfparams); + +/** + * 专用于在Python端释放C中申请的内存 + * + * @param pt 指针 + */ +void free1d(void *pt); \ No newline at end of file diff --git a/pygrt/C_extension/src/dynamic/signals.c b/pygrt/C_extension/src/dynamic/signals.c index f86a1868..226f2732 100644 --- a/pygrt/C_extension/src/dynamic/signals.c +++ b/pygrt/C_extension/src/dynamic/signals.c @@ -348,3 +348,8 @@ float * get_custom_wave(int *Nt, const char *tfparams){ *Nt = nt; return tfarr; } + + +void free1d(void *pt){ + free(pt); +} \ No newline at end of file diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 88e5d539..3a90a8c1 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -13,6 +13,9 @@ from .c_structures import * +FPOINTER = POINTER(c_float) +IPOINTER = POINTER(c_int) + c_PGRN = POINTER(c_GRN) @@ -81,4 +84,29 @@ def set_num_threads(n): C_compute_travt1d.argtypes = [ PREAL, PREAL, c_int, c_int, c_int, REAL -] \ No newline at end of file +] + + + +# ------------------------------------------------------------------- +# C函数定义的时间函数 +# ------------------------------------------------------------------- +C_free = libgrt.free1d +"""释放在C中申请的内存""" +C_free.restype = None +C_free.argtypes = [c_void_p] + +C_get_trap_wave = libgrt.get_trap_wave +"""梯形波""" +C_get_trap_wave.restype = FPOINTER +C_get_trap_wave.argtypes = [c_float, FPOINTER, FPOINTER, FPOINTER, IPOINTER] + +C_get_parabola_wave = libgrt.get_parabola_wave +"""抛物波""" +C_get_parabola_wave.restype = FPOINTER +C_get_parabola_wave.argtypes = [c_float, FPOINTER, IPOINTER] + +C_get_ricker_wave = libgrt.get_ricker_wave +"""雷克子波""" +C_get_ricker_wave.restype = FPOINTER +C_get_ricker_wave.argtypes = [c_float, c_float, IPOINTER] \ No newline at end of file diff --git a/pygrt/pygrn.py b/pygrt/pygrn.py index 634eae7f..cf894217 100755 --- a/pygrt/pygrn.py +++ b/pygrt/pygrn.py @@ -153,7 +153,7 @@ def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0): sac = self.SACTrace nt = sac.npts dt = sac.delta - wI = self.wI + wI = sac.user0 T = nt*dt if not np.isclose(T*df, 1.0): diff --git a/pygrt/signals.py b/pygrt/signals.py index 57523b21..ed4975e9 100755 --- a/pygrt/signals.py +++ b/pygrt/signals.py @@ -10,187 +10,92 @@ import numpy as np +import numpy.ctypeslib as npct + +from .c_interfaces import * __all__ = [ "gen_triangle_wave", - "gen_square_wave", "gen_parabola_wave", - "gen_dirac_wave", - "gen_step_wave", - "gen_ramp_wave", - "gen_smooth_ramp_wave", + "gen_trap_wave", "gen_ricker_wave", ] -def gen_triangle_wave(nlen, vlen, dt): +def gen_triangle_wave(vlen, dt): ''' 生成三角信号 - :param nlen: 总时长(s) :param vlen: 信号时长(s) :param dt: 采样间隔(s) :return: - - **t** - 时间序列 - **wave** - 波形幅值序列 ''' - t = np.arange(0, nlen+dt, dt) - triangle_wave = np.zeros_like(t) - - vnt = int(vlen/dt)+1 - phase = np.linspace(0, 1, vnt) - triangle_wave[0:vnt] = (1 - 2 * np.abs(phase - 0.5)) - - return t, triangle_wave - -def gen_square_wave(nlen, vlen, dt): - ''' - 生成矩形信号 - - :param nlen: 总时长(s) - :param vlen: 信号时长(s) - :param dt: 采样间隔(s) + return gen_trap_wave(vlen/2.0, vlen/2.0, vlen, dt) - :return: - - **t** - 时间序列 - - **wave** - 波形幅值序列 - ''' - t = np.arange(0, nlen+dt, dt) - square_wave = np.zeros_like(t) - - square_wave[0:int(vlen/dt)+1] = 1 - - return t, square_wave - -def gen_parabola_wave(nlen, vlen, dt): + +def gen_parabola_wave(vlen, dt): ''' 生成抛物线信号 - :param nlen: 总时长(s) :param vlen: 信号时长(s) :param dt: 采样间隔(s) :return: - - **t** - 时间序列 - **wave** - 波形幅值序列 ''' - t = np.arange(0, nlen+dt, dt) - parabola_wave = np.zeros_like(t) - - vnt = int(vlen/dt)+1 - phase = np.linspace(0, 1, vnt) - parabola_wave[0:vnt] = 1 - (phase - 0.5) ** 2 * 4 - - return t, parabola_wave - -def gen_dirac_wave(nlen, dt): - ''' - 生成脉冲信号 + ct1 = c_float(vlen) + cnt = c_int(0) - :param nlen: 总时长(s) - :param dt: 采样间隔(s) - - :return: - - **t** - 时间序列 - - **wave** - 波形幅值序列 - ''' - t = np.arange(0, nlen+dt, dt) - wave = np.zeros_like(t) - wave[0] = 1 - return t, wave + carr = C_get_parabola_wave(dt, byref(ct1), byref(cnt)) + arr = npct.as_array(carr, shape=(cnt.value,)).copy() -def gen_step_wave(nlen, dt): - ''' - 生成阶跃信号(全1信号) + C_free(carr) - :param nlen: 总时长(s) - :param dt: 采样间隔(s) + return arr - :return: - - **t** - 时间序列 - - **wave** - 波形幅值序列 +def gen_trap_wave(t1, t2, t3, dt): ''' - t = np.arange(0, nlen+dt, dt) - step_wave = np.ones_like(t) - return t, step_wave - -def gen_ramp_wave(nlen, t0, dt): - ''' - 生成斜坡信号 + 生成梯形信号 - :param nlen: 总时长(s) - :param t0: 上坡时长(s) + :param t1: 上坡截止时刻(s) + :param t2: 平台截止时刻(s) + :param t3: 下坡截止时刻(s) :param dt: 采样间隔(s) :return: - - **t** - 时间序列 - **wave** - 波形幅值序列 ''' - t = np.arange(0, nlen+dt, dt) - ramp_wave = np.ones_like(t) - - vnt = int(t0/dt)+1 - phase = np.linspace(0, 1, vnt) - ramp_wave[:vnt] = phase - return t, ramp_wave - -def gen_smooth_ramp_wave(nlen, t0, dt): - ''' - 生成光滑斜坡信号 + ct1 = c_float(t1) + ct2 = c_float(t2) + ct3 = c_float(t3) + cnt = c_int(0) - :param nlen: 总时长(s) - :param t0: 上坡(1-1/e)时长(s) - :param dt: 采样间隔(s) + carr = C_get_trap_wave(dt, byref(ct1), byref(ct2), byref(ct3), byref(cnt)) + arr = npct.as_array(carr, shape=(cnt.value,)).copy() - :return: - - **t** - 时间序列 - - **wave** - 波形幅值序列 - ''' - t = np.arange(0, nlen+dt, dt) - wave = (1 - np.exp(-t/t0)) - return t, wave + C_free(carr) -def gen_trap_wave(nlen, t1, t2, t3, dt): - ''' - 生成梯形信号 + return arr - :param nlen: 总时长(s) - :param t1: 上坡截止时刻(s) - :param t2: 平台截止时刻(s) - :param t3: 下坡截止时刻(s) - :param dt: 采样间隔(s) - :return: - - **t** - 时间序列 - - **wave** - 波形幅值序列 - ''' - t = np.arange(0, nlen+dt, dt) - wave = np.piecewise( - t, - [(t > 0 )* (t < t1), (t >= t1) * (t < t2), (t >= t2) * (t < t3), t >= t3], - [lambda t: t / t1, - 1, - lambda t: (t3 - t)/(t3 - t2), - 0] - ) - return t, wave - -def gen_ricker_wave(nlen, f0:float, dt:float): +def gen_ricker_wave(f0:float, dt:float): ''' 生成Ricker子波 - :param nlen: 总时长(s) :param f0: 中心频率(Hz) :param dt: 采样间隔(s) :return: - - **t** - 时间序列 - **wave** - 波形幅值序列 ''' + cnt = c_int(0) + + carr = C_get_ricker_wave(dt, f0, byref(cnt)) + if cast(carr, c_void_p).value is None: + raise ValueError("NULL pointer") + arr = npct.as_array(carr, shape=(cnt.value,)).copy() - t = np.arange(0, nlen+dt, dt) - t0 = 1.0/f0 - a = np.pi**2 * f0**2 * (t-t0)**2 - wave = (1 - 2*a ) * np.exp(-a) + C_free(carr) - return t, wave + return arr \ No newline at end of file diff --git a/pygrt/utils.py b/pygrt/utils.py index bd2ae423..5ad08f77 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -24,6 +24,8 @@ import math import os +from numpy.typing import ArrayLike + from .c_structures import * @@ -284,7 +286,7 @@ def gen_syn_from_gf_EXP(st:Stream, M0:float, az:float, calc_upar:bool=False): return gen_syn_from_gf(st, calc_upar, "COMPUTE_EXP", M0, az) -def gen_syn_from_gf_MT(st:Stream, M0:float, MT:np.ndarray, az:float, calc_upar:bool=False): +def gen_syn_from_gf_MT(st:Stream, M0:float, MT:ArrayLike, az:float, calc_upar:bool=False): ''' 矩张量源，单位为dyne*cm @@ -364,22 +366,41 @@ def __check_trace_attr_sac(tr:Trace, **kwargs): # return st -def stream_convolve(st0:Stream, signal:np.ndarray, inplace=True): +def stream_convolve(st0:Stream, signal0:np.ndarray, inplace=True): ''' 对stream中每一道信号做线性卷积 :param st0: 记录多个Trace的 :class:`obspy.Stream` 类型 - :param signal: 卷积信号 + :param signal0: 卷积信号 :param inplace: 是否做原地更改 :return: - **stream** - 处理后的结果, :class:`obspy.Stream` 类型 ''' st = st0 if inplace else deepcopy(st0) + signal = deepcopy(signal0) for tr in st: data = tr.data - data[:] = oaconvolve(data, signal, mode='full')[:data.shape[0]] + dt = tr.stats.delta + + fac = None + user_wI = hasattr(tr.stats, "sac") and "user0" in tr.stats.sac + # 使用虚频率先压制 + if user_wI: + npts = tr.stats.npts + wI = tr.stats.sac['user0'] + fac = np.exp(np.arange(0, npts)*dt*wI) + signal = deepcopy(signal0) + + signal[:] /= fac[:len(signal)] + data[:] /= fac + + data1 = np.pad(data, (len(signal)-1, 0), mode='wrap') # 强制循环卷 + data[:] = oaconvolve(data1, signal, mode='valid')[:data.shape[0]] * dt # dt是连续卷积的系数 + + if user_wI: + data[:] *= fac return st @@ -433,7 +454,9 @@ def stream_write_sac(st:Stream, path:str): 将一系列Trace以SAC形式保存到本地，以发震时刻作为参考0时刻 :param st: 记录多个Trace的 :class:`obspy.Stream` 类型 - :param path: 保存文件名，不要加后缀 + :param path: 保存文件名，不要加后缀 + 例如path="GRN/res"表明sac文件将保存在GRN文件中， + 文件名分别为resR.sac, resT.sac, resZ.sac ''' # 新建对应文件夹 @@ -443,7 +466,7 @@ def stream_write_sac(st:Stream, path:str): # 每一道的保存路径为path+{channel} for tr in st: - tr.write(".".join([path, tr.stats.channel]) + '.sac', format='SAC') + tr.write(path+tr.stats.channel[-1] + '.sac', format='SAC')