<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Delay Estimation
OpendTect has the Match Delta atttribute to measure time shifts between similar events in different seismic volumes. It works by finding the peaks in each volume and then tries to match them up. The algorithm is simple and fast but the results can be quite noisy. The following is a the result of running Match Delta between a depth section and itself shifted up by 13 metres. 
The purpose of this notebook is to test some other options that could then be used in an External Attribute script.
Simple Test Signal
Step1:  Local Normalised Cross-correlation
Here is a mainly numpy version of local normalised cross-correlation.
Step2:  Can we get a faster result using Numba?
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
%pylab inline
import test_signals as tst
def make_signals(nsamp,delay ):
    ref = np.random.rand(nsamp+abs(delay))*2-1
    wav = sig.ricker(80,5)
    filtered = np.convolve(ref, wav,'same')
    if delay < 0 :
        return filtered[0:nsamp], filtered[-delay:nsamp-delay]
    else:
        return filtered[delay:nsamp+delay], filtered[0:nsamp]
    
#res, shifted = make_signals(1000,10)
res, shifted = tst.make_delayed_signal_pair(1000,10)
fig = plt.figure(figsize=(12,1))
plt.plot(res)
plt.plot(shifted)
shifted.shape
Explanation: Delay Estimation
OpendTect has the Match Delta atttribute to measure time shifts between similar events in different seismic volumes. It works by finding the peaks in each volume and then tries to match them up. The algorithm is simple and fast but the results can be quite noisy. The following is a the result of running Match Delta between a depth section and itself shifted up by 13 metres. 
The purpose of this notebook is to test some other options that could then be used in an External Attribute script.
Simple Test Signal
End of explanation
#
# Rolling summ of squares
def rollingSSQ(inp, winlen):
    inpsq = np.square(inp)
    kernel = np.ones(winlen)
    return np.convolve(inpsq, kernel, 'same')
    
def localCorr_numpy(reference, match, winlen, nlag,lag,qual):
    hwin = winlen//2
    lags = 2*nlag+1
    ns = reference.shape[0]
    hxw = hwin-nlag
    cor = np.zeros(lags)
    refSSQ = np.sqrt(rollingSSQ(reference,2*hxw+1))
    matSSQ = np.sqrt(rollingSSQ(match,2*hxw+1))
    for ir in range(hwin,ns-hwin):
        rbeg = ir - hxw
        rend = ir + hxw + 1
        mbeg = rbeg - nlag
        mend = rend + nlag 
        cor = np.divide(np.correlate(match[mbeg:mend],reference[rbeg:rend],'valid'),matSSQ[ir-nlag:ir+nlag+1]*refSSQ[ir])
        pos = np.argmax(cor)
        if pos>0 and pos<lags-1:
            cp = (cor[pos-1]-cor[pos+1])/(2.*cor[pos-1]-4.*cor[pos]+2.*cor[pos+1])
            lag[ir] = pos-nlag+cp
            qual[ir] = cor[pos]
        else:
            lag[ir]=0.0
            qual[ir]=0.0
lag = np.zeros(res.shape)
qual = np.zeros(res.shape)
localCorr_numpy(res,shifted,51,15,lag,qual)
fig = plt.figure(figsize=(12,2))
plt.plot(lag)
plt.plot(qual)
%timeit -o localCorr_numpy(res,shifted,51,15,lag,qual)
Explanation: Local Normalised Cross-correlation
Here is a mainly numpy version of local normalised cross-correlation.
End of explanation
import sys,os
from numba import jit
sys.path.insert(0, os.path.join(sys.path[0], '..'))
import extnumba as xn
@jit(nopython=True)
def localCorr_numba(reference, match, winlen, nlag,lag,qual):
    hwin = winlen//2
    lags = 2*nlag+1
    ns = reference.shape[0]
    hxw = hwin-nlag
    cor = np.zeros(lags)
    refSSQ = np.zeros(ns)
    matSSQ = np.zeros(ns)
    xn.winSSQ(reference,2*hxw+1,refSSQ)
    xn.winSSQ(match,2*hxw+1,matSSQ)
    for ir in range(hwin,ns-hwin):
        rbeg = ir - hxw
        rend = ir + hxw + 1
        mbeg = rbeg - nlag
        mend = rend + nlag
        for il in range(lags):
            lbeg = rbeg + il - nlag
            lend = lbeg + 2 * hxw + 1
            sum = 0.0
            for iref,imat in zip(range(rbeg,rend),range(lbeg,lend)):
                sum += reference[iref]*match[imat]
            den = refSSQ[ir]*matSSQ[lbeg+hxw]
            if den== 0.0:
                cor[il] = 0.0
            else:
                cor[il] = sum/den
        pos = np.argmax(cor)
        if pos>0 and pos<lags-1:
            cp = (cor[pos-1]-cor[pos+1])/(2.*cor[pos-1]-4.*cor[pos]+2.*cor[pos+1])
            lag[ir] = pos-nlag+cp
            qual[ir] = cor[pos]
        else:
            lag[ir]=0.0
            qual[ir]=0.0
localCorr_numba(res,shifted,51,15,lag,qual)
fig = plt.figure(figsize=(12,2))
plt.plot(lag)
plt.plot(qual)
%timeit -o localCorr_numba(res,shifted,51,15,lag,qual)
Explanation: Can we get a faster result using Numba?
End of explanation
<END_TASK>