首页 > 解决方案 > 为什么在填充输入之一时python scipy互相关不起作用

问题描述

scipy 互相关函数根本不适用于特定的一维数组,我不知道为什么。下面的代码演示了这个问题,只需使用一个跟踪而不是另一个跟踪即可。

这个问题有点与互相关Python互相关没有返回正确的移位有关

#!/usr/bin/python3

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

def _main():
    """
    trace = np.array([0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, # down the step
                      0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, # up the step
                      0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002]) # down the step
    """
    trace = np.array([0.51231204949426460, 0.47472182808002383, 0.48806029762272723, 0.51352464310119930, 0.58506742537603330, 0.62993314829830390, 0.57657927012749040, 0.55369158834668990, 0.56255864527226200, 0.61576098682569510,
                      0.62955418648769630, 0.64236215760241170, 0.69063835641941580, 0.75073729780384960, 0.86896478361172370, 0.92216712516515690, 0.91329988783884970, 0.92807831604813670, 0.99113300320800610, 0.99999999999999990, 0.91527040506699960, 
                      0.80098377331469030, 0.71723934679539750, 0.68275634764039450, 0.65812563395824950, 0.63250963159524040, 0.59999708953480900, 0.55172083058422660, 0.54975037348965490, 0.57011178351142090, 0.52807534544936740])


    left_padded_trace = np.pad(trace, (10, 0), mode='constant', constant_values=trace.min())
    center_padded_trace = np.pad(trace, (5, 5), mode='constant', constant_values=trace.min())
    right_padded_trace = np.pad(trace, (0, 10), mode='constant', constant_values=trace.min())

    correlation1 = signal.correlate(center_padded_trace, left_padded_trace, mode='full', method='fft')
    correlation2 = signal.correlate(center_padded_trace, center_padded_trace, mode='full', method='fft')
    correlation3 = signal.correlate(center_padded_trace, right_padded_trace, mode='full', method='fft')

    corr_peak_index1 = np.argmax(correlation1)
    corr_max1 = np.max(correlation1)

    corr_peak_index2 = np.argmax(correlation2)
    corr_max2 = np.max(correlation2)

    corr_peak_index3 = np.argmax(correlation3)
    corr_max3 = np.max(correlation3)

    offset1 = corr_peak_index1-(center_padded_trace.size-1)
    offset2 = corr_peak_index2-(center_padded_trace.size-1)
    offset3 = corr_peak_index3-(center_padded_trace.size-1)

    print("Corr1: {}, Corr2: {}, Corr3: {}".format(corr_peak_index1, corr_peak_index2, corr_peak_index3))
    print("Offset1: {}, Offset2: {}, Offset3: {}".format(offset1, offset2, offset3))

    plt.figure(1)

    plt.subplot(311)
    plt.plot(range(0, center_padded_trace.size), center_padded_trace, 'r-',
            range(offset1, left_padded_trace.size+offset1), left_padded_trace, 'b--',
            range(0, correlation1.size), correlation1/corr_max1, 'g-',
            [corr_peak_index1], [1], 'k+')

    plt.subplot(312)
    plt.plot(range(0, center_padded_trace.size), center_padded_trace, 'r-',
            range(offset2, center_padded_trace.size+offset2), center_padded_trace, 'b--',
            range(0, correlation2.size), correlation2/corr_max2, 'g-',
            [corr_peak_index2], [1], 'k+')

    plt.subplot(313)
    plt.plot(range(0, center_padded_trace.size), center_padded_trace, 'r-',
            range(offset3, right_padded_trace.size+offset3), right_padded_trace, 'b--',
            range(0, correlation3.size), correlation3/corr_max3, 'g-',
            [corr_peak_index3], [1], 'k+')

    plt.show()

由于填充添加的移位是相同的,唯一的区别是输入轨迹上的变化,因此相关性的移位和对齐结果应该是相同的,但事实并非如此。

对于第一个跟踪(更多合成步骤),相关性和偏移量是:(1 是左侧填充,2 是居中,3 是右侧填充)

对于第二个轨迹(更自然),

跟随情节:

解决方案

请参阅下面 Paul Panzer 的回答和评论。

问题在于原始代码与非零填充有关。

当用非零值移动数组时,互相关值越来越高,峰值受到影响。以下代码和图像演示了这种效果:

    trace = np.array([0.51231204949426460, 0.47472182808002383, 0.48806029762272723, 0.51352464310119930, 0.58506742537603330, 0.62993314829830390, 0.57657927012749040, 0.55369158834668990, 0.56255864527226200, 0.61576098682569510, 0.62955418648769630, 0.64236215760241170, 0.69063835641941580, 0.75073729780384960, 0.86896478361172370, 0.92216712516515690, 0.91329988783884970, 0.92807831604813670, 0.99113300320800610, 0.99999999999999990, 0.91527040506699960, 0.80098377331469030, 0.71723934679539750, 0.68275634764039450, 0.65812563395824950, 0.63250963159524040, 0.59999708953480900, 0.55172083058422660, 0.54975037348965490, 0.57011178351142090, 0.52807534544936740])

    for padding_value in np.arange(0, trace.min(), trace.min()/10):
        left_padded_trace = np.pad(trace, (10, 0), mode='constant', constant_values=padding_value)
        center_padded_trace = np.pad(trace, (5, 5), mode='constant', constant_values=padding_value)

        correlation = signal.correlate(center_padded_trace, left_padded_trace, mode='full', method='fft')

        corr_peak_index = np.argmax(correlation)

        plt.figure(2)
        plt.subplot(211)
        plt.title('Left Padded Trace')
        plt.xticks([])
        plt.plot(left_padded_trace)
        plt.subplot(212)
        plt.title('Centered Padded Trace')
        plt.plot(center_padded_trace)

        plt.figure(3)
        plt.plot(range(0, correlation.size), correlation)
        plt.plot([corr_peak_index], [correlation[corr_peak_index]], 'k+')

    plt.show()

结果如下所示。可以看出,随着填充值的增加,相关峰值向中心移动。

标签: pythonnumpyscipycross-correlation

解决方案


差异可以通过以下事实来解释,即在第二次跟踪的情况下,您正在填充不为零的最小值。因此,您不能期望峰值会随着偏移量而移动。相反,您会得到移动的峰值曲线加上一个与最小值成比例的三角形。

概括

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

def _main(offset=0, trace_idx=0):
    trace = [np.array([0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, # down the step
                      0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, 0.99999999999999999, 0.99999999999999998, # up the step
                       0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002, 0.00000000000000001, 0.00000000000000002]), # down the step
    np.array([0.51231204949426460, 0.47472182808002383, 0.48806029762272723, 0.51352464310119930, 0.58506742537603330, 0.62993314829830390, 0.57657927012749040, 0.55369158834668990, 0.56255864527226200, 0.61576098682569510,
                      0.62955418648769630, 0.64236215760241170, 0.69063835641941580, 0.75073729780384960, 0.86896478361172370, 0.92216712516515690, 0.91329988783884970, 0.92807831604813670, 0.99113300320800610, 0.99999999999999990, 0.91527040506699960, 
                      0.80098377331469030, 0.71723934679539750, 0.68275634764039450, 0.65812563395824950, 0.63250963159524040, 0.59999708953480900, 0.55172083058422660, 0.54975037348965490, 0.57011178351142090, 0.52807534544936740])][trace_idx]

    trace += offset - trace.min()

    left_padded_trace = np.pad(trace, (10, 0), mode='constant', constant_values=trace.min())
    center_padded_trace = np.pad(trace, (5, 5), mode='constant', constant_values=trace.min())
    right_padded_trace = np.pad(trace, (0, 10), mode='constant', constant_values=trace.min())

    correlation1 = signal.correlate(center_padded_trace, left_padded_trace, mode='full', method='fft')
    correlation2 = signal.correlate(center_padded_trace, center_padded_trace, mode='full', method='fft')
    correlation3 = signal.correlate(center_padded_trace, right_padded_trace, mode='full', method='fft')

    corr_peak_index1 = np.argmax(correlation1)
    corr_max1 = np.max(correlation1)

    corr_peak_index2 = np.argmax(correlation2)
    corr_max2 = np.max(correlation2)

    corr_peak_index3 = np.argmax(correlation3)
    corr_max3 = np.max(correlation3)

    offset1 = corr_peak_index1-(center_padded_trace.size-1)
    offset2 = corr_peak_index2-(center_padded_trace.size-1)
    offset3 = corr_peak_index3-(center_padded_trace.size-1)

    return offset1, offset2, offset3

    print("Corr1: {}, Corr2: {}, Corr3: {}".format(corr_peak_index1, corr_peak_index2, corr_peak_index3))
    print("Offset1: {}, Offset2: {}, Offset3: {}".format(offset1, offset2, offset3))

    plt.figure(1)

    plt.subplot(311)
    plt.plot(range(0, center_padded_trace.size), center_padded_trace, 'r-',
            range(offset1, left_padded_trace.size+offset1), left_padded_trace, 'b--',
            range(0, correlation1.size), correlation1/corr_max1, 'g-',
            [corr_peak_index1], [1], 'k+')

    plt.subplot(312)
    plt.plot(range(0, center_padded_trace.size), center_padded_trace, 'r-',
            range(offset2, center_padded_trace.size+offset2), center_padded_trace, 'b--',
            range(0, correlation2.size), correlation2/corr_max2, 'g-',
            [corr_peak_index2], [1], 'k+')

    plt.subplot(313)
    plt.plot(range(0, center_padded_trace.size), center_padded_trace, 'r-',
            range(offset3, right_padded_trace.size+offset3), right_padded_trace, 'b--',
            range(0, correlation3.size), correlation3/corr_max3, 'g-',
            [corr_peak_index3], [1], 'k+')

    plt.show()


x = np.arange(200)*0.01
y1 = np.array([*map(_main, x)])

y2 = np.array([*map(_main, x, np.ones(x.size,int))])

plt.figure(1)
plt.subplot(211)
plt.title('synthetic')
plt.plot(x,y1)
plt.legend(('left-shifted input', 'centered input', 'right-shifted input'))
plt.subplot(212)
plt.title('natural')
plt.plot(x,y2)
plt.ylabel('x-offset of result')
plt.xlabel('y-offset')
plt.savefig("summary.png")

推荐阅读