查找音频文件中摘录的开始位置：Python 中两个数组之间的互相关系数

Question

我一直在为一个问题烦恼，到目前为止我在 StackOverflow 上找到的所有答案都没有帮助 - 所以我请求你的帮助。

总体问题

我想创建一个函数来找到较大音频文件中音频摘录开始的确切时间戳。出于测试目的，我使用了 5 分钟的音频文件和 43 秒的摘录。下面，我在 Audacity 中对齐了两个音频文件：摘录正好从 00:01:55.554920 开始。

我还希望函数返回一个值，当且仅当它的“置信度值”超过某个阈值时，这实际上是函数的参数。我打算这样做的方法是检查两个对齐信号之间的相关系数是否超过给定的阈值。 换句话说，这是代码的简化版本：

find_excerpt_starting_sample(original_audio, excerpt, threshold):

    # Find the cross-correlation coefficients for each lag
    xcorr = cross_correlation(original_audio, excerpt)

    # Return the lag of the max correlation if it is over threshold
    if np.max(xcorr) > threshold:
        return np.argmax(xcorr)
    else:
        raise Exception("No correlation over threshold found.")

我在寻找正确的 

cross_correlation

函数时遇到了很多麻烦，因为我的所有尝试都没有返回一个介于 0 和 1 之间的数组。

问题，简化了

由于我对音频文件的尝试尚未得出结论，因此我尝试对两个数值数组执行相同的操作：

y1 = [2, 22, 14, 8, 0, 4, 8, 16, 26, 6, 12, 14, 16, 2, 6]
y2 = [4, 8, 16, 26, 6, 12]

这里，y2 包含 y1 的子集（从索引 5 开始）。为了确保该函数独立于幅度比例工作，我将 y2 的所有值减半：

y1 = [2, 22, 14, 8, 0, 4, 8, 16, 26, 6, 12, 14, 16, 2, 6]
y2 = [2, 4, 8, 13, 3, 6]

我想创建一个互相关函数，返回一个数组，其中滞后 5 处的值为 1。

到目前为止我的尝试

np.corrcoef

如果我们只是做一个简单的关联并沿着原始音频滑动摘录，它就会起作用：

import numpy as np
import matplotlib as plt

corr = np.zeros(len(y1) - len(y2))
for i in range(len(y1) - len(y2)):
    corr[i] = np.corrcoef(y1[i:i+len(y2)], y2)[0][1]

print(corr)
plt.plot(corr)
plt.show()

输出为：

[ 0.18961375 -0.71250433 -0.56075283 -0.08468414 0.21913077 1. -0.04179451 -0.46803451 -0.24815461]

问题是这种技术对于较长的文件来说真的非常非常不高效。

scipy.signal.correlate

现在，我不再重新发明轮子，而是开始使用 Stack Overflow 上找到的主要解决方案之一，即 scipy.signal 的关联函数。它返回找到适当滞后的值。然而，因为它执行的是卷积，所以无法量化相关性。

from scipy import signal

xcorr = signal.correlate(y1, y2, mode="full")
lags = signal.correlation_lags(len(y1), len(y2), mode="full")

print(xcorr)
print(lags)
plt.plot(lags, xcorr)
plt.show()

输出为：

[ 12 138 176 392 390 332 224 232 356 402 596 486 478 414 422 252 186  88 28  12]
[-5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14]

我看到了一些解决方案，但它们并没有按照我的预期工作。

首先，解决方案

here

建议使用此函数对系数进行归一化：

corr = signal.correlate(y1 / np.std(y1), y2 / np.std(y2), 'full') / min(len(y1), len(y2))
lags = signal.correlation_lags(len(y1), len(y2), mode="full")
print(c)
plt.plot(lags, c) 
plt.show() 

输出为：

[0.0736392 0.84685082 1.08004163 2.40554727 2.39327407 2.03735126 1.37459844 1.42369124 2.18462966 2.46691327 3.65741371 2.98238769 2.93329489 2.54055247 2.58964528 1.54642325 1.14140763 0.54002082 0.17182481 0.0736392 ]

如您所见，最大值不是 1，而是 3.65741371。

然后我尝试了

here

找到的另一种解决方案：

y1n = y1 / np.std(y1)
y2n = y2 / np.std(y2)
xcorr = signal.correlate(y1n, y2n, mode="full")
lags = signal.correlation_lags(len(y1), len(y2), mode="full")
print(xcorr)
plt.plot(lags, xcorr)
plt.show()

输出为：

[ 0.44183521 5.08110495 6.48024979 14.43328362 14.35964442 12.22410756 8.24759064 8.54214745 13.10777798 14.80147963 21.94448224 17.89432612 17.59976931 15.24331484 15.53787165 9.27853947 6.8484458 3.24012489 1.03094883 0.44183521]

再次强调，互相关的最大值不是 1，而是 21.94448224

呼救声

关于相关性，我有很多不了解的地方 - 我深入研究了它，但在深入研究之前，我问，你们中的一个人是否能够为我指出正确的方向，以及到目前为止我做错了什么.

非常感谢！

Answer 1

def cross_corr(x, y):
  x = np.array(x)
  y = np.array(y[::-1])
  yi = (y - y.mean())/ y.std() / np.sqrt(y.size)
  x_m = np.convolve(x, np.ones(yi.size), 'valid')**2/yi.size
  x_m2 = np.convolve(x**2, np.ones(yi.size), 'valid')
  return np.convolve(x, yi, 'valid')/np.sqrt(x_m2 - x_m)

cross_corr(y1,y2)
array([ 0.18961375, -0.71250433, -0.56075283, -0.08468414,  0.21913077,
        1.        , -0.04179451, -0.46803451, -0.24815461,  0.77716484])

这个函数比原来的解决方案快了很多数量级

Answer 2

signal.correlate()

替换卷积。这是我得到的：

import numpy as np

def cross_correlation(y1, y2):
    y2_normalized = (y2 - y2.mean()) / y2.std() / np.sqrt(y2.size)
    y1_m = signal.correlate(y1, np.ones(y2.size), 'valid') ** 2 / y2_normalized.size
    y1_m2 = signal.correlate(y1 ** 2, np.ones(y2.size), "valid")
    cross_correlation = signal.correlate(y1, y2_normalized, "valid") / np.sqrt(y1_m2 - y1_m)

交叉关联音频文件时，代码速度要快得多；对于我在问题开始时提到的音频剪辑，代码的执行时间不到 2 秒（分别为 5 分钟和 43 秒，以 44100 Hz 采样）。

这是输出：

峰值值为1.0，精确到毫秒。

注意：在执行互相关之前，我使用

y_env = np.abs(scipy.signal.hilbert(y))

获得了两个音频文件的包络，并使用

b, a = scipy.signal.butter(2, filter_over, "low", fs=44100)

和

y_filt = lfilter(b, a, y_env)

获得了 50 Hz 的低通滤波器。如果您的数据很长，您还可以随时执行下采样。