对时间戳数组使用互相关有意义吗?

Does it make sense to use cross-correlation on arrays of timestamps?

我想找到两个时间戳数组之间的偏移量。比方说,它们可以代表两个音轨中的哔哔声。

注意:任一音轨中可能有额外或缺失的开始。

我发现了一些看起来很有希望的关于互相关的信息(例如 https://dsp.stackexchange.com/questions/736/how-do-i-implement-cross-correlation-to-prove-two-audio-files-are-similar)。

我假设每个音轨的持续时间为 10 秒,并将蜂鸣声开始表示为 "square wave" 的峰值,采样率为 44.1 kHz:

import numpy as np

rfft = np.fft.rfft
irfft = np.fft.irfft

track_1 = np.array([..., 5.2, 5.5, 7.0, ...])
# The onset in track_2 at 8.0 is "extra," it has no
# corresponding onset in track_1
track_2 = np.array([..., 7.2, 7.45, 8.0, 9.0, ...])
frequency = 44100
num_samples = 10 * frequency
wave_1 = np.zeros(num_samples)
wave_1[(track_1 * frequency).astype(int)] = 1
wave_2 = np.zeros(num_samples)
wave_2[(track_2 * frequency).astype(int)] = 1
xcor = irfft(rfft(wave_1) * np.conj(rfft(wave_2)))
offset = xcor.argmax()

这种方法不是特别快,但即使频率很低,我也能得到相当一致的结果。但是……我不知道这是不是个好主意!有没有比互相关更好的方法来找到这个偏移量?

编辑:添加了关于缺失和额外发作的注释。

是的,它确实有道理。这通常在 matlab 中完成。这是一个类似应用程序的 link:

http://www.mathworks.com/help/signal/ug/cross-correlation-of-delayed-signal-in-noise.html

几个注意事项

当相关信号噪声过多时,通常会使用互相关。如果您不需要担心任何噪音,我会使用不同的方法。

如果track_1track_2是哔哔声的时间戳并且都捕捉到所有的哔哔声,那么就不需要构建波形并进行互相关。只需找到两个蜂鸣时间戳数组之间的平均延迟:

import numpy as np

frequency = 44100
num_samples = 10 * frequency
actual_time_delay = 0.020
timestamp_noise = 0.001

# Assumes track_1 and track_2 are timestamps of beeps, with a delay and noise added
track_1 = np.array([1,2,10000])
track_2 = np.array([1,2,10000]) + actual_time_delay*frequency + 
          frequency*np.random.uniform(-timestamp_noise,timestamp_noise,len(track_1))

calculated_time_delay = np.mean(track_2 - track_1) / frequency
print('Calculated time delay (s):',calculated_time_delay)

这会产生大约 0.020 秒,具体取决于引入的随机错误,并且速度会尽可能快。

编辑: 如果需要处理额外或缺失的时间戳,则可以按如下方式修改代码。本质上,如果时间戳随机性的误差是有界的,则对所有值之间的延迟执行统计 "mode" 函数。时间戳随机性范围内的任何事物都聚集在一起并被识别,然后对原始识别的延迟进行平均。

import numpy as np
from scipy.stats import mode

frequency = 44100
num_samples = 10 * frequency
actual_time_delay = 0.020

# Assumes track_1 and track_2 are timestamps of beeps, with a delay, noise added,
#   and extra/missing beeps
timestamp_noise = 0.001
timestamp_noise_threshold = 0.002
track_1 = np.array([1,2,5,10000,100000,200000])
track_2 = np.array([1,2,44,10000,30000]) + actual_time_delay*frequency 
track_2 = track_2 + frequency*np.random.uniform(-timestamp_noise,timestamp_noise,len(track_2))

deltas = []
for t1 in track_1:
    for t2 in track_2:
        deltas.append( t2 - t1)
sorted_deltas = np.sort(deltas)/frequency
truncated_deltas = np.array(sorted_deltas/(timestamp_noise_threshold)+0.5, 
dtype=int)*timestamp_noise_threshold

truncated_time_delay = min(mode(truncated_deltas).mode,key=abs)
calculated_time_delay = np.mean(sorted_deltas[truncated_deltas == truncated_time_delay])

print('Calculated time delay (s):',calculated_time_delay)

显然,如果缺少太多蜂鸣声或存在额外的蜂鸣声,那么代码会在某个时候开始失败,但通常它应该运行良好并且比尝试生成整个波形和执行相关要快得多。