首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >指定numpy.correlate的移位

指定numpy.correlate的移位
EN

Stack Overflow用户
提问于 2015-01-13 01:47:29
回答 2查看 2.8K关注 0票数 2

我想知道是否有可能为两个一维阵列的互相关指定由k变量表示的移位。因为通过将numpy.correlate函数及其mode参数设置为'full‘,我将获得所取数组的整个长度的每个k移位的互相关系数(假设两个数组大小相同)。让我在下面的例子中向你展示我的确切意思:

代码语言:javascript
复制
import numpy as np

# Define signal 1.
signal_1 = np.array([1, 2 ,3])

# Define signal 2.
signal_2 = np.array([1, 2, 3])

# Other definitions.
Xi = signal_1
Yi = signal_2
N = np.size(Xi)
k = 3
Xs = np.average(Xi)
Ys = np.average(Yi)

# Cross-covariance coefficient function.
def crossCovariance(Xi, Yi, N, k, Xs, Ys, forCorrelation = False):
    autoCov = 0
    for i in np.arange(0, N-k):
        autoCov += ((Xi[i+k])-Xs)*(Yi[i]-Ys)
    if forCorrelation == True:
        return autoCov/N
    else:
        return (1/(N-1))*autoCov

# Expected value function.
def E(X, P):
    expectedValue = 0
    for i in np.arange(0, np.size(X)):
        expectedValue += X[i] * (P[i] / np.size(X))
    return expectedValue 

# Cross-correlation coefficient function.
def crossCorrelation(Xi, Yi, k):
    # Calculate the covariance coefficient.
    cov = crossCovariance(Xi, Yi, N, k, Xs, Ys, forCorrelation = True)

    # Calculate standard deviations.
    EX = E(Xi, np.ones(np.size(Xi))) 
    SDX = (E((Xi - EX) ** 2, np.ones(np.size(Xi)))) ** (1/2)

    EY = E(Yi, np.ones(np.size(Yi)))
    SDY = (E((Yi - EY) ** 2, np.ones(np.size(Yi)))) ** (1/2)

    # Calculate correlation coefficient.
    return cov / (SDX * SDY)

# Express cross-covariance or cross-correlation function in a form of a 1D vector.
def array(k, norm = True):
    # If norm = True, return array of autocorrelation coefficients.
    # If norm = False, return array of autocovariance coefficients.
    vector = np.array([])
    shifts = np.abs(np.arange(-k, k+1, 1))
    for i in shifts:
        if norm == True:
            vector = np.append(crossCorrelation(Xi, Yi, i), vector)
        else:
            vector = np.append(crossCovariance(Xi, Yi, N, i, Xs, Ys), vector)
    return vector

在我的示例中,为k的不同值调用方法array(k, norm = True)将产生resuslt,如下所示:

代码语言:javascript
复制
k = 3, [ 0.  -0.5  0.   1.   0.  -0.5  0. ]
k = 2, [-0.5  0.   1.   0.  -0.5]
k = 1, [ 0.  1.  0.]
k = 0, [ 1.]

我的方法对于学习目的来说是很好的,但为了加快分析速度,我需要转移到本机numpy函数。如何在使用原生numpy.correlate函数的同时指定k移位值?PS k参数指定两个数组之间的“时间”移位。提前谢谢你。

EN

回答 2

Stack Overflow用户

发布于 2015-01-19 20:03:13

虽然我不知道有任何内置的函数来计算特定范围的信号滞后的互相关,但你可以通过向量化来加速你的版本,即在数组上执行操作,而不是数组中的单个元素。

此版本仅在lags上使用单个Python循环:

代码语言:javascript
复制
import numpy as np


def xcorr(x, y, k, normalize=True):

    n = x.shape[0]

    # initialize the output array
    out = np.empty((2 * k) + 1, dtype=np.double)
    lags = np.arange(-k, k + 1)

    # pre-compute E(x), E(y)
    mu_x = x.mean()
    mu_y = y.mean()

    # loop over lags
    for ii, lag in enumerate(lags):

        # use slice indexing to get 'shifted' views of the two input signals
        if lag < 0:
            xi = x[:lag]
            yi = y[-lag:]
        elif lag > 0:
            xi = x[:-lag]
            yi = y[lag:]
        else:
            xi = x
            yi = y

        # x - mu_x; y - mu_y
        xdiff = xi - mu_x
        ydiff = yi - mu_y

        # E[(x - mu_x) * (y - mu_y)]
        out[ii] = xdiff.dot(ydiff) / n

        # NB: xdiff.dot(ydiff) == (xdiff * ydiff).sum()

    if normalize:
        # E[(x - mu_x) * (y - mu_y)] / (sigma_x * sigma_y)
        out /=  np.std(x) * np.std(y)

    return lags, out

一些更一般的建议:

正如我在评论中提到的那样,你应该试着给你的函数起一个信息丰富的名字,并且不太可能与你的命名空间中的其他东西发生冲突(例如np.array).

  • It's
  • array ),这样你的函数才是自包含的。在您的版本中,NkXsYs是在main函数外部定义的。在这种情况下,您可能会意外地修改或覆盖这些变量中的一个,并且调试由这种thing.
  • Appending导致的错误可能会变得很棘手(例如,使用np.appendnp.concatenate) numpy数组的速度很慢,因此请尽可能避免这种情况。如果像本例一样,提前知道输出的大小,那么预先分配输出数组(例如使用np.emptynp.zeros),然后填充元素会快得多。如果您一定要进行连接,则通常更快的方法是附加到一个普通的Python数组,然后在最后将其转换为list数组。
票数 2
EN

Stack Overflow用户

发布于 2015-12-13 11:02:48

可以通过指定maxlags来使用它

代码语言:javascript
复制
import matplotlib.pyplot as plt
xcorr = plt.xcorr(signal_1, signal_2, maxlags=1)

文档可以在here上找到。这个实现是基于np.correlate的。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/27907866

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档