obspy.signal.cross_correlation.correlate_template

obsp.signal의 correlate 함수와는 다른 목적으로 쓰임. 

2026.09.09 - [Python/obpsy] - obspy.signal.cross_correlation.correlate

obspy.signal.cross_correlation.correlate

 

correlate(a,b, shift) ---- 두 비슷한 길이의 신호사이의 상대 shift 찾기 (b가 a에 비해 몇 sample 이동했나?)

correlate_template(data, template) ---- 긴 data안에서 짧은 template이 어디에 있는지 찾기 (이 template이 data에 어디에 있나?)

 

correlate_template은 shift를 지정하지 않고, template이 data의 전체위에서 움직인다. 

data
------------------------------------------------------------

template
██████

가능한 위치:

██████
  ██████
    ██████
      ██████
         ...
                              ██████

 

def _window_sum(data, window_len):
    """Rolling sum of data."""
    window_sum = np.cumsum(data)
    # in-place equivalent of
    # window_sum = window_sum[window_len:] - window_sum[:-window_len]
    # return window_sum
    np.subtract(window_sum[window_len:], window_sum[:-window_len],
                out=window_sum[:-window_len])
    return window_sum[:-window_len]

def correlate_template(data, template, mode='valid', normalize='full',
                       demean=True, method='auto'):
    # if we get Trace objects, use their data arrays
    if isinstance(data, Trace):
        data = data.data
    if isinstance(template, Trace):
        template = template.data
    data = np.asarray(data)
    template = np.asarray(template)
    lent = len(template)
    if len(data) < lent:
        raise ValueError('Data must not be shorter than template.')
    if demean:
        template = template - np.mean(template)
        if normalize != 'full':
            data = data - np.mean(data)
            
            
    cc = scipy.signal.correlate(data, template, mode=mode, method=method)
    if normalize is not None:
        tnorm = np.sum(template ** 2)
        if normalize == 'naive':
            norm = (tnorm * np.sum(data ** 2)) ** 0.5
            if norm <= np.finfo(float).eps:
                cc[:] = 0
            elif cc.dtype == float:
                cc /= norm
            else:
                cc = cc / norm
        elif normalize == 'full':
            pad = len(cc) - len(data) + lent
            if mode == 'same':
                pad1, pad2 = (pad + 2) // 2, (pad - 1) // 2
            else:
                pad1, pad2 = (pad + 1) // 2, pad // 2
            data = _pad_zeros(data, pad1, pad2)
            # in-place equivalent of
            # if demean:
            #     norm = ((_window_sum(data ** 2, lent) -
            #              _window_sum(data, lent) ** 2 / lent) * tnorm) ** 0.5
            # else:
            #      norm = (_window_sum(data ** 2, lent) * tnorm) ** 0.5
            # cc = cc / norm
            if demean:
                norm = _window_sum(data, lent) ** 2
                if norm.dtype == float:
                    norm /= lent
                else:
                    norm = norm / lent
                np.subtract(_window_sum(data ** 2, lent), norm, out=norm)
            else:
                norm = _window_sum(data ** 2, lent)
            norm *= tnorm
            if norm.dtype == float:
                np.sqrt(norm, out=norm)
            else:
                norm = np.sqrt(norm)
            mask = norm <= np.finfo(float).eps
            if cc.dtype == float:
                cc[~mask] /= norm[~mask]
            else:
                cc = cc / norm
            cc[mask] = 0
        else:
            msg = "normalize has to be one of (None, 'naive', 'full')"
            raise ValueError(msg)
    return cc

 

따라서 mode='valid' 인것이 합리적이고 당연하고 실제로 default 값이다. 당연히 template의 길이가 data 보다 짧아야한다.

예를 들어서.

data length     = 1000
template length = 100

#mode='valid'의 결과 길이는 1000-1000+1 = 901

cc[0]은 template와 data[0:100]의 비교 결과
cc[1]은 template와 data[1:101]의 비교 결과
cc[450]은 template와 data[450:550]의 비교결과

 

 

 

Normalization (full)

각 template 위치마다 data window의 normalization을 새로 계산 

 

그래서 full normalization은 위치마다 denominator가 다르다. 따라서 진폭이 달라도 shape가 같으면 cc=1이 될 수 있다.