如何将存在间隔和不同时间基准的两个时间序列相关联？

Question

我有两个3D加速度计的时间序列数据，它们具有不同的时间基准（时钟在不同时间开始，采样时间中有一些非常轻微的偏移），并且包含许多不同大小的间隙（由于写入不同闪存设备的延迟造成）。

我使用的加速度计是廉价的GCDC X250-2。我将加速度计的增益调到最高点，因此数据有显著的噪声水平。

每个时间序列约有2百万个数据点（每秒512个样本，持续一个小时），并且包含大约500个感兴趣的事件，其中典型事件跨越100-150个样本（每个200-300毫秒）。在闪存写入期间，许多事件受到数据中断的影响。

因此，数据不完美，甚至不太好看。但是我的肉眼检查清楚地显示了包含我感兴趣信息的数据。（如果需要，我可以发布图表。）

加速度计处于类似的环境中，但只是适度耦合，这意味着我可以通过肉眼判断哪些来自每个加速度计的事件相匹配，但是我目前在软件中还没有成功。由于物理限制，这些设备也安装在不同的方向上，其中轴不匹配，但它们尽可能靠近正交。因此，例如，对于3轴加速度计A和B，+ Ax映射到- By（上下），+ Az映射到- Bx（左右），+ Ay映射到- Bz（前后）。

我的最初目标是在垂直轴上相关冲击事件，尽管我最终想要a)自动发现轴映射，b)相关映射轴上的活动，以及c)提取两个加速度计之间的行为差异（例如扭曲或弯曲）。

由于时间序列数据的特性，Python的numpy.correlate()无法使用。我还尝试了R的Zoo包，但没有取得任何进展。我已经寻求不同领域的信号分析帮助，但没有取得任何进展。

有人有什么线索可以提供给我，或者我应该研究哪些方法？

2021年2月28日更新：在此处添加了一些图表here，显示数据的示例。

- BobC

---

听起来你想要（以编程方式）找到两个值：A和B之间的时间偏移量以及A和B之间的时间比例。就这些吗？这些间隙是否会对记录的时间产生影响（例如，时钟重新开始计时），还是它们只是缺失数据的时间段？ - Paul

@Paul：这些间隙只是缺失的数据，而且在两个设备之间是独立的。 - BobC

好问题！但它可能更适合于http://stats.stackexchange.com。 - Björn Pollex

@Space_C0wb0y：好主意！我已经在那里发布了**这个问题**。 - BobC

Answer 1

我对你的问题的理解是：给定两个非常长且嘈杂的时间序列，找到一个移位，使得一个信号中的大“突起”与另一个信号中的大突起匹配。

我的建议是：将数据插值为均匀间隔，矫正和平滑数据（假设快速振荡的相位不重要），并进行逐点互相关（假设小的移位将对齐数据）。

import numpy
from scipy.ndimage import gaussian_filter
"""
sig1 and sig 2 are assumed to be large, 1D numpy arrays
sig1 is sampled at times t1, sig2 is sampled at times t2
t_start, t_end, is your desired sampling interval
t_len is your desired number of measurements
"""

t = numpy.linspace(t_start, t_end, t_len)
sig1 = numpy.interp(t, t1, sig1)
sig2 = numpy.interp(t, t2, sig2)
#Now sig1 and sig2 are sampled at the same points.

"""
Rectify and smooth, so 'peaks' will stand out.
This makes big assumptions about your data;
these assumptions seem true-ish based on your plots.
"""
sigma = 10 #Tune this parameter to get the right smoothing
sig1, sig2 = abs(sig1), abs(sig2)
sig1, sig2 = gaussian_filter(sig1, sigma), gaussian_filter(sig2, sigma)

"""
Now sig1 and sig2 should look smoothly varying, with humps at each 'event'.
Hopefully we can search a small range of shifts to find the maximum of the 
cross-correlation. This assumes your data are *nearly* lined up already.
"""
max_xc = 0
best_shift = 0
for shift in range(-10, 10): #Tune this search range
    xc = (numpy.roll(sig1, shift) * sig2).sum()
    if xc > max_xc:
        max_xc = xc
        best_shift = shift
print 'Best shift:', best_shift
"""
If best_shift is at the edges of your search range,
you should expand the search range.
"""

- Andrew

---

我会看一下，但由于每个外部事件在每个设备上的表现不同（称一个设备为“原始”的，另一个设备为“过滤”的），我对匹配的成功持悲观态度。 - BobC

@BobC 你试过了吗？我怀疑只要平滑的程度适当，相关性仍然有效。 - Andrew

我还在尝试中，但是还没有得到可用的结果。 - BobC

好的，那太糟糕了。如果有一些简单的方法可以分享两个样本数据集，我会自己试试看。 - Andrew

Answer 2

如果数据包含未知大小的间隙，而且每个时间序列的间隙都不同，那么我会放弃尝试将整个序列相关联，而是尝试在每个时间序列上交叉相关配对短窗口对，例如重叠窗口长度为典型事件的两倍（300个样本长）。找到所有可能的高交叉相关匹配，并对潜在匹配施加顺序约束以获取匹配窗口的序列。

从那里开始，您就有了更容易分析的较小问题。

- hotpaw2

---

听起来像是我想要的，但问题接下来就变成了如何确定窗口。下一个问题是，这两个数据集包含不同的数据，其中一个可以被认为是“原始”的，另一个则是通过被监测系统进行了“过滤”的。而我不知道滤波函数是什么，虽然这是我想要确定的东西。峰值很少在个别事件之间达成一致。此外，该系统有一个大约为180赫兹的外部施加的杂乱背景信号，它在视觉上很容易区分，但变化太大，不能轻易地去除。 - BobC

从固定长度的重叠窗口开始，比如说每100个样本开始300个长。会有一些冗余。你提出的其他问题应该放在它们自己单独的DSP问题中。 - hotpaw2

Answer 3

这不是一个技术性的答案，但它可能会帮助你想出一个：

• 将图形转换为图像，并将其插入到像gimp或photoshop这样的良好图像程序中
• 每当有间隙时，将图形分成离散的图像
• 将第一系列图形放在水平线上
• 将第二个系列放在它正下方的水平线上
• 直观地识别第一个相关事件
• 如果两个事件没有垂直对齐：
  • 选择更靠左的实例和该行右侧的所有内容
  • 将这些内容向右拖动，直到它们对齐

这基本上就是音频编辑器的工作原理，因此，如果您将其转换为像未压缩的WAV文件这样的简单音频格式，您可以直接在Audacity之类的软件中操作它。（当然，它会听起来很糟糕，但您可以轻松地移动数据图形。）

实际上，Audacity也有一种名为Nyquist的脚本语言，因此如果您不需要程序来“检测”相关性（或者至少愿意暂时推迟该步骤），您可能可以使用Audacity标记和Nyquist的某些组合来自动对齐并导出干净的数据，一旦您标记了相关点，就可以选择以您喜欢的格式。

- tangentstorm

---

这通常是一个有趣的想法，但它增加了分析信号的复杂度。大多数采用这种方法的人通常会将其转换为功率谱密度并比较PSD图像-即使维度增加了，域更改通常也会带来净收益。 - Mark Elliot

我只是提出使用图像作为一个思维实验，但我认为Audacity实际上是一个相当明智的选择。这正是音频编辑程序旨在解决的问题类型，而当你加入Nyquist（它基本上是用于音频/时间序列数据的Scheme语言）时，这就是一个非常强大的组合。 - tangentstorm

我正在使用pyplot做类似的事情，并且在我发布的另一个问题的答案中，我可以让图形一起缩放但单独平移。因此，我知道数据集中第一个共同事件的时钟差约为4.78秒。 - BobC

Answer 4

我猜你需要手动构建一个偏移表，以对齐系列之间的“匹配项”。以下是获取这些匹配项的一种方法示例。思路是将数据左右移动直到对齐，然后调整比例直到“匹配”。试试吧。

library(rpanel)

#Generate the x1 and x2 data
n1 <- rnorm(500)
n2 <- rnorm(200)
x1 <- c(n1, rep(0,100), n2, rep(0,150))
x2 <- c(rep(0,50), 2*n1, rep(0,150), 3*n2, rep(0,50))

#Build the panel function that will draw/update the graph
lvm.draw <- function(panel) {
       plot(x=(1:length(panel$dat3))+panel$off, y=panel$dat3, ylim=panel$dat1, xlab="", ylab="y", main=paste("Alignment Graph   Offset = ", panel$off, "   Scale = ", panel$sca, sep=""), typ="l")
       lines(x=1:length(panel$dat3), y=panel$sca*panel$dat4, col="red")
       grid()
       panel
}

#Build the panel
xlimdat <- c(1, length(x1))
ylimdat <- c(-5, 5)
panel <- rp.control(title = "Eye-Ball-It", dat1=ylimdat, dat2=xlimdat, dat3=x1, dat4=x2, off=100, sca=1.0, size=c(300, 160))
rp.slider(panel, var=off, from=-500, to=500, action=lvm.draw, title="Offset", pos=c(5, 5, 290, 70), showvalue=TRUE)
rp.slider(panel, var=sca, from=0, to=2, action=lvm.draw, title="Scale", pos=c(5, 70, 290, 90), showvalue=TRUE)

- bill_080

---

我已经在pyplot中做了类似的事情。问题是，对于数百个事件，手动处理过程很痛苦。 - BobC

Answer 5

听起来你想要最小化函数(Ax'+By) + (Az'+Bx) + (Ay'+Bz)的一对值：即时间偏移量t₀和时间比例因子t_r，其中Ax' = t_r*(Ax + t₀)等。

我建议使用SciPy的双变量优化函数。并且我会在“间隙”（假设可以通过编程确定间隙）上使用掩码或暂时将数据（例如Ax'和By）置零。

为了使过程更加高效，从A和B的粗略采样开始，但设置与采样相称的fmin（或您选择的任何优化器）的精度。然后继续进行全数据集的逐渐细分窗口，直到您的窗口变窄且未被下采样。

编辑-匹配轴

关于尝试确定哪个轴与给定轴共线的问题，而不知道数据特征的情况下，我可以指向一个类似的问题。请查看pHash或本文中概述的其他方法，以帮助识别相似的波形。

- Paul

---

我会研究你推荐的函数。如果可能的话，我也想自动确定设备之间正确的轴配对。 - BobC

@BobC：听起来你知道 Ax 大致与 By 同轴？我理解错了吗？ - Paul

@Paul：是的，加速度计的轴大致平行，可能相差不到几度。但Ax与By对齐，因此加速度计相对于彼此有两个90度的旋转。 - BobC

@BobC：那么，“我还想自动确定正确的轴配对”这句话是什么意思，我就不知道了。 - Paul

@Paul：给定一对未知方向但轴平行的三轴装置数据，根据数据内容将轴配对。 - BobC

@BobC：明白了，听起来Ax是与-By配对的。所有排列组合都可以吗？（六面立方体，每面四个方向）总共24种？ - Paul