Skip to content

FEAT: support skip amplitude compensation from imaginary frequency #219

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Dengda98 merged 1 commit into from
Jan 17, 2026
Merged

FEAT: support skip amplitude compensation from imaginary frequency #219

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 7 additions & 3 deletions docs/source/Module/greenfn.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ greenfn
**grt greenfn**
|-M|\ *model*
|-D|\ *depsrc/deprcv*
|-N|\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**]
|-N|\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**][**+f**]
|-R|\ *file*\|\ *r1,r2,...*
|-O|\ *outdir*
[ |-H|\ *f1/f2* ]
Expand Down Expand Up @@ -56,7 +56,7 @@ greenfn

.. _-N:

**-N**\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**]
**-N**\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**][**+f**]
采样点数 *nt* 和采样时间间隔 *dt* (secs) ，这将决定计算的最高频。还可设置：

+ **+w**\ *zeta* - 虚频率系数 :math:`\zeta` [0.8]。 |bouchon1981| 提出在频率上添加一个虚频率偏移，
Expand All @@ -70,7 +70,11 @@ greenfn
默认情况下，程序会跳过非常低频的几个点以避免引入误差，
详见 :doc:`/Advanced/k_integ/drift/waveform_drift` 的介绍。

当时窗长度 nt\*dt 太小“包不住”有效信号，或时窗长度足够但时延(|-E|)不合适，输出的波形会发生混叠，
+ **+f** - 不进行由虚频率引起的振幅补偿。
默认情况下，程序会对时间域结果乘上 :math:`\exp(\frac{\zeta \pi}{T} t)` 进行振幅补偿，
但当时窗长度远超有效信号的长度时，这种补偿会放大尾部噪声。

**注：** 当时窗长度 nt\*dt 太小“包不住”有效信号，或时窗长度足够但时延(|-E|)不合适，输出的波形会发生混叠，
此时需调整 |-N| 和 |-E| 。

.. include:: explain_-R.rst_
Expand Down
30 changes: 21 additions & 9 deletions pygrt/C_extension/src/dynamic/grt_greenfn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@ typedef struct {
real_t wI; ///< 虚频率 zeta*PI/r
real_t *freqs;
size_t upsample_n; ///< 升采样倍数
bool keepAllFreq;
bool keepAllFreq; ///< 计算所有频率，不论频率多低
bool skipImagComps; ///< 跳过虚频率的补偿
} N;
/** 输出目录 */
struct {
Expand Down Expand Up @@ -316,7 +317,7 @@ printf("\n"
" <depsrc>: source depth (km).\n"
" <deprcv>: receiver depth (km).\n"
"\n"
" -N<nt>/<dt>[+w<zeta>][+n<fac>][+a] \n"
" -N<nt>/<dt>[+w<zeta>][+n<fac>][+a][+f] \n"
" <nt>: number of points. (NOT requires 2^n).\n"
" <dt>: time interval (secs). \n"
" +w<zeta>: define the coefficient of imaginary \n"
Expand All @@ -328,6 +329,8 @@ printf("\n"
" and calculated frequencies stay unchanged.\n"
" +a: All frequencies are calculated regardless of\n"
" how low the frequency is.\n"
" +f: skip the amplitude compensation from \n"
" imaginary frequency.\n"
"\n"
" -R<r1>,<r2>[,...]\n"
" Multiple epicentral distance (km), \n"
Expand Down Expand Up @@ -508,7 +511,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
}
break;

// 点数,采样间隔,虚频率 -Nnt/dt[+w<zeta>][+n<scale>][+a]
// 点数,采样间隔,虚频率 -Nnt/dt[+w<zeta>][+n<scale>][+a][+f]
case 'N':
Ctrl->N.active = true;
{
Expand Down Expand Up @@ -546,6 +549,10 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
case 'a':
Ctrl->N.keepAllFreq = true;
break;

case 'f':
Ctrl->N.skipImagComps = true;
break;

default:
GRTBadOptionError(N, "+%s is not supported.", token);
Expand Down Expand Up @@ -848,7 +855,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
static void write_one_to_sac(
const char *srcname, const char ch, GRT_FFTW_HOLDER *fh, const real_t wI,
SACTRACE *sac, const char *s_output_subdir, const char *s_prefix,
const int sgn, const cplx_t *grncplx)
const int sgn, bool skipImagComps, const cplx_t *grncplx)
{
snprintf(sac->hd.kcmpnm, sizeof(sac->hd.kcmpnm), "%s%s%c", s_prefix, srcname, ch);

Expand Down Expand Up @@ -879,12 +886,17 @@ static void write_one_to_sac(
// 归一化，并处理虚频
// 并转为 SAC 需要的单精度类型
real_t fac, coef;
coef = fh->df * exp(sac->hd.b * wI);
fac = exp(wI*fh->dt);
coef = fh->df;
fac = 1.0;
if (! skipImagComps) {
coef *= exp(sac->hd.b * wI);
fac = exp(wI*fh->dt);
}
for(size_t i = 0; i < fh->nt; ++i){
sac->data[i] = fh->w_t[i] * coef;
coef *= fac;
}


// 以sac文件保存到本地
grt_write_SACTRACE(s_outpath, sac);
Expand Down Expand Up @@ -1117,11 +1129,11 @@ int greenfn_main(int argc, char **argv) {

char ch = GRT_ZRT_CODES[c];

write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "", sgn, grn[ir][im][c]);
write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "", sgn, Ctrl->N.skipImagComps, grn[ir][im][c]);

if(Ctrl->e.active){
write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "z", sgn*(-1), grn_uiz[ir][im][c]);
write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "r", sgn, grn_uir[ir][im][c]);
write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "z", sgn*(-1), Ctrl->N.skipImagComps, grn_uiz[ir][im][c]);
write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "r", sgn , Ctrl->N.skipImagComps, grn_uir[ir][im][c]);
}
}

Expand Down
5 changes: 3 additions & 2 deletions pygrt/pygrn.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ def plot_response(self):
return fig, (ax1, ax2)


def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0):
def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0, skipImagComps:bool=False):
'''
Convert the Green's function from the frequency domain to the time domain
and return it in the form of :class:`obspy.Trace`
Expand Down Expand Up @@ -142,7 +142,8 @@ def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0):
# 实序列的傅里叶变换
data = irfft(cmlx_grn, nt, norm='backward') * (1/dt) # *(1/dt)和连续傅里叶变换幅值保持一致
# 抵消虚频率的影响
data *= np.exp((np.arange(0,nt)*dt + T0)*wI)
if not skipImagComps:
data *= np.exp((np.arange(0,nt)*dt + T0)*wI)

# 保存sac头段变量
sac.o = 0.0
Expand Down
8 changes: 5 additions & 3 deletions pygrt/pymod.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -410,6 +410,7 @@ def compute_grn(
converg_method:Union[str,None]=None,
delayT0:float=0.0,
delayV0:float=0.0,
skipImagComps:bool=False,
calc_upar:bool=False,
gf_source=['EX', 'VF', 'HF', 'DC'],
statsfile:Union[str,None]=None,
Expand Down Expand Up @@ -439,6 +440,7 @@ def compute_grn(
:param safilonTol: precision of Self-Adaptive Filon's Integration Method
:param filonCut: The splitting point of DWM and (SA)FIM, k*=<filonCut>/rmax, default is 0
:param converg_method: The method of explicit convergence, you can set "DCM", "PTAM" or "none". Default use "DCM" when abs(depsrc-deprcv) <= 1.0 km
:param skipImagComps: skip the amplitude compensation from imaginary frequency.
:param calc_upar: whether calculate the spatial derivatives of displacements.
:param gf_source: The source type to be calculated
:param statsfile: directory path for saving the statsfile during k integral, used to debug or observe the variations of :math:`F(k,\omega)` and :math:`F(k,\omega)J_m(kr)k`
Expand Down Expand Up @@ -513,10 +515,10 @@ def compute_grn(
continue

sgn = -1 if ZRTchs[c]=='Z'=='Z' else 1
stream.append(pygrnLst[ir][im][c].freq2time(delayT, travtP, travtS, sgn ))
stream.append(pygrnLst[ir][im][c].freq2time(delayT, travtP, travtS, sgn, skipImagComps))
if(calc_upar):
stream.append(pygrnLst_uiz[ir][im][c].freq2time(delayT, travtP, travtS, sgn*(-1) ))
stream.append(pygrnLst_uir[ir][im][c].freq2time(delayT, travtP, travtS, sgn ))
stream.append(pygrnLst_uiz[ir][im][c].freq2time(delayT, travtP, travtS, sgn*(-1), skipImagComps))
stream.append(pygrnLst_uir[ir][im][c].freq2time(delayT, travtP, travtS, sgn , skipImagComps))


# 在sac头段变量部分
Expand Down