两个时间依赖的多维信号（信号向量）的相关性

Question

两个时间依赖的多维信号（信号向量）的相关性

pythonnumpy

5

我有一个矩阵M1，其中每一行都是一个时变信号。

我还有另一个矩阵M2，与M1具有相同的维度，其中每一行也是一个时变信号，作为“模板”，用于识别第一个矩阵中的信号形状。

我的目标是得到一个列向量v，其中v [i]是M1的第i行和M2的第i行之间的相关性。

我查看了numpy的corrcoef函数，并尝试了以下代码：

import numpy as np

M1 = np.array ([
    [1, 2, 3, 4],
    [2, 3, 1, 4]
])

M2 = np.array ([
    [10, 20, 30, 40],
    [20, 30, 10, 40]
])

print (np.corrcoef (M1, M2))

这会输出：

[[ 1.   0.4  1.   0.4]
 [ 0.4  1.   0.4  1. ]
 [ 1.   0.4  1.   0.4]
 [ 0.4  1.   0.4  1. ]]

我一直在阅读文档，但仍然不清楚我需要选择哪些矩阵条目作为我的向量v的条目。

有人可以帮忙吗？

（我已经学习了几个类似问题的S.O.答案，但还没有看到光...）

代码背景：

有256行（信号），我在“主信号”上运行一个窗口大小为200的滑动窗口，其长度为10k个样本。因此，M1和M2都是256行x200列。对不起，关于10k个样本的错误。那是总信号长度。通过使用滑动模板的相关性，我尝试找到最佳匹配模板的偏移量。实际上，我正在寻找256通道侵入性心电图（或者如医生所称的心电图）中的QRS波群。

    lg.info ('Processor: {}, time: {}, markers: {}'.format (self.key, dt.datetime.now ().time (), len (self.data.markers)))

    # Compute average signal shape over preexisting markers and uses that as a template to find the others.
    # All generated markers will have the width of the widest preexisting one.

    template = np.zeros ((self.data.samples.shape [0], self.bufferWidthSteps))

    # Add intervals that were marked in advance
    nrOfTerms = 0
    maxWidthSteps = 0
    newMarkers = []
    for marker in self.data.markers:
        if marker.key == self.markerKey:

            # Find start and stop sample index    
            startIndex = marker.tSteps - marker.stampWidthSteps // 2
            stopIndex = marker.tSteps + marker.stampWidthSteps // 2

            # Extract relevant slice from samples and add it to template
            template += np.hstack ((self.data.samples [ : , startIndex : stopIndex], np.zeros ((self.data.samples.shape [0], self.bufferWidthSteps - marker.stampWidthSteps))))

            # Adapt nr of added terms to facilitate averaging
            nrOfTerms += 1

            # Remember maximum width of previously marked QRS complexes
            maxWidthSteps = max (maxWidthSteps, marker.stampWidthSteps)
        else:
            # Preexisting markers with non-matching keys are just copied to the new marker list
            # Preexisting markers with a matching key are omitted from the new marker list
            newMarkers.append (marker)

    # Compute average of intervals that were marked in advance
    template = template [ : , 0 : maxWidthSteps] / nrOfTerms
    halfWidthSteps = maxWidthSteps // 2

    # Append markers of intervals that yield an above threshold correlation with the averaged marked intervals
    firstIndex = 0
    stopIndex = self.data.samples.shape [1] - maxWidthSteps
    while firstIndex < stopIndex:
        corr = np.corrcoef (
            template,
            self.data.samples [ : , firstIndex : firstIndex + maxWidthSteps]
        )

        diag = np.diagonal (
            corr,
            template.shape [0]
        )

        meanCorr = np.mean (diag)

        if meanCorr > self.correlationThreshold:
            newMarkers.append ([self.markerFactories [self.markerKey] .make (firstIndex + halfWidthSteps, maxWidthSteps)])

            # Prevent overlapping markers
            firstIndex += maxWidthSteps
        else:
            firstIndex += 5

    self.data.markers = newMarkers

    lg.info ('Processor: {}, time: {}, markers: {}'.format (self.key, dt.datetime.now ().time (), len (self.data.markers)))

- Jacques de Hooge

3个回答

0

我认为是这样的：（如果不正确，请纠正！）

import numpy as np

M1 = np.array ([
    [1, 2, 3, 4],
    [2, 3, 1, 4.5]
])

M2 = np.array ([
    [10, 20, 33, 40],
    [20, 35, 15, 40]
])

v = np.diagonal (np.corrcoef (M1, M2), M1.shape [0])

print (v)

输出什么：

[ 0.99411402  0.96131896]

因为它只有一个维度，我可以把它看作是一个列向量...

- Jacques de Hooge

0

由于我不太熟悉numpy数组的使用方法，所以我会挑选出行，将每个对传递给corrcoeff函数进行处理。

[np.corrcoef(i,j)[0][1] for i,j in zip(a,b)]

对于一个 np.array 列的输出

c, c.shape = np.array([np.corrcoef(i,j)[0][1] for i,j in zip(a,b)]), (a.shape[0], 1)

我相信使用numpy的广播/索引功能会更好

- f5r5e5d

我考虑过这个问题，尽管出于性能原因我不想在Python中使用循环。不过我认为性能下降应该是很小的，因为大部分时间都在计算corcoeff本身。 - Jacques de Hooge

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Divakar · Accepted Answer

基于此解决方案，用于查找两个2D数组之间的相关矩阵，我们可以有一个类似的方法来查找计算两个数组中对应行之间的相关向量。实现看起来会像这样 -

def corr2_coeff_rowwise(A,B):
    # Rowwise mean of input arrays & subtract from input arrays themeselves
    A_mA = A - A.mean(1)[:,None]
    B_mB = B - B.mean(1)[:,None]

    # Sum of squares across rows
    ssA = (A_mA**2).sum(1);
    ssB = (B_mB**2).sum(1);

    # Finally get corr coeff
    return np.einsum('ij,ij->i',A_mA,B_mB)/np.sqrt(ssA*ssB)

我们可以进一步优化这部分内容，通过在此引入einsum魔法来获得ssA和ssB！

def corr2_coeff_rowwise2(A,B):
    A_mA = A - A.mean(1)[:,None]
    B_mB = B - B.mean(1)[:,None]
    ssA = np.einsum('ij,ij->i',A_mA,A_mA)
    ssB = np.einsum('ij,ij->i',B_mB,B_mB)
    return np.einsum('ij,ij->i',A_mA,B_mB)/np.sqrt(ssA*ssB)

示例运行 -

In [164]: M1 = np.array ([
     ...:     [1, 2, 3, 4],
     ...:     [2, 3, 1, 4.5]
     ...: ])
     ...: 
     ...: M2 = np.array ([
     ...:     [10, 20, 33, 40],
     ...:     [20, 35, 15, 40]
     ...: ])
     ...: 

In [165]: corr2_coeff_rowwise(M1, M2)
Out[165]: array([ 0.99411402,  0.96131896])

In [166]: corr2_coeff_rowwise2(M1, M2)
Out[166]: array([ 0.99411402,  0.96131896])

运行时测试 -

In [97]: M1 = np.random.rand(256,200)
    ...: M2 = np.random.rand(256,200)
    ...: 

In [98]: out1 = np.diagonal (np.corrcoef (M1, M2), M1.shape [0])
    ...: out2 = corr2_coeff_rowwise(M1, M2)
    ...: out3 = corr2_coeff_rowwise2(M1, M2)
    ...: 

In [99]: np.allclose(out1, out2)
Out[99]: True

In [100]: np.allclose(out1, out3)
Out[100]: True

In [101]: %timeit np.diagonal (np.corrcoef (M1, M2), M1.shape [0])
     ...: %timeit corr2_coeff_rowwise(M1, M2)
     ...: %timeit corr2_coeff_rowwise2(M1, M2)
     ...: 
100 loops, best of 3: 9.5 ms per loop
1000 loops, best of 3: 554 µs per loop
1000 loops, best of 3: 430 µs per loop

einsum 相对于内置的 np.corrcoef，可以实现 20x+ 倍速提升！