import numpy as np
from numba import jit
@jit()
def astack(array, norm_p = 3, nsi = 1, roll_scale = np.array([-20,20])):
    nch, npts = array.shape
    ref_trace = np.zeros(npts)
    out_data  = np.zeros_like(array)
    for i in range(nsi):
        tau_ = np.zeros(nch)
        if i == 0:
            data_ = np.zeros_like(array)
            for ch, tr in enumerate(array):
                data_[ch] = tr/np.sqrt(np.sum(np.square(tr)))  #6s, 300個(gè)點(diǎn)做l2歸一化
                ref_trace += tr/np.sqrt(np.sum(np.square(tr))) #6s, 歸一化的reference trace
            ref_trace = ref_trace/nch
            tau = np.zeros(nch)
            
            ref_4s = ref_trace[50:-50]  #4s的reference trace
            data_4s = np.zeros((nch, 200))
            for n, tr in enumerate(data_):
                P_tr = np.zeros(roll_scale[1] - roll_scale[0])
                for j, shift_samp in enumerate(np.arange(roll_scale[0], roll_scale[1])):
                    tr_roll = np.roll(tr, shift_samp)[50:-50]
                    P_tr[j] = np.linalg.norm((ref_4s - tr_roll), norm_p)
                tau[n] = np.argmin(P_tr)
            
            for b, idx in enumerate(tau):
                tau_[b] = np.arange(roll_scale[0], roll_scale[1])[int(idx)]
#             np.arange(roll_scale[0], roll_scale[1])[]
            
#             tau_ = np.array([np.arange(roll_scale[0], roll_scale[1])[int(idx)] for idx in tau])
            
#     ref_4s = np.zeros(200)
#     for i, tr in enumerate(data_):
#         tr_roll = np.roll(tr, tau_[i])[50:-50]
#         ref_4s+=tr_roll
#     ref_4s = ref_4s/nch
            
    return tau_, ref_4s
            

還沒(méi)有寫多個(gè)循環(huán)，加入遞歸的操作。該算法的目的是把沒(méi)有對(duì)齊的波形自適應(yīng)拉齊。
參考自Rawlinson, N., & Kennett, B. L. (2004). Rapid estimation of relative and absolute delay times across a network by adaptive stacking. Geophysical Journal International, 157(1), 332-340.