如何将两个时间序列与间隙和不同的时基相关联？答案

【问题标题】：How to correlate two time series with gaps and different time bases?如何将两个时间序列与间隙和不同的时基相关联？
【发布时间】：2011-07-05 02:05:29
【问题描述】：

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

我使用的加速度计是便宜的GCDC X250-2。我正在以最高增益运行加速度计，因此数据具有显着的本底噪声。

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

因此，数据不是原始的，甚至不是很漂亮。但是我的眼球检查显示它清楚地包含我感兴趣的信息。（如果需要，我可以发布情节。）

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

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

时间序列数据的性质使得 Python 的 numpy.correlate() 无法使用。我也看过 R 的 Zoo 包，但没有取得任何进展。我曾向信号分析的不同领域寻求帮助，但没有取得任何进展。

有人对我能做什么或我应该研究的方法有任何线索吗？

2011 年 2 月 28 日更新：添加了一些图表 here 显示数据示例。

【问题讨论】：

听起来您想（以编程方式）找到两个值：A 和 B 之间的时间偏移以及 A 和 B 之间的时间比例。仅此而已吗？间隔是否对记录的时间有任何影响（例如，时钟重新开始）还是只是缺少数据的时期？
@Paul：差距只是缺少数据，在两个设备之间是独立的。
好问题！但它可能更适合stats.stackexchange.com。
@Space_C0wb0y：好主意！我已经在那里发布了 this question。

标签： python r numpy signal-processing accelerometer

【解决方案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.
"""

【讨论】：

我会看看这个，但由于每个设备对每个外部事件的看法都不同（将一个设备称为“原始”，另一个称为“过滤”），我对匹配并不乐观工作。
@BobC 你试过了吗？我的怀疑是，在正确的平滑水平下，相关性仍然有效。
我还在玩它，但还没有得到可用的结果。
好吧，那太糟糕了。如果有一些直接的方法可以共享两个示例数据集，我会自己尝试一下。

【解决方案2】：

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

从那里你有更容易分析的小问题。

【讨论】：

这听起来像我想要的，但问题就变成了如何确定窗口。下一个问题是这两个数据集包含不同的数据，其中一个可以被认为是“原始的”，而另一个被监控的系统“过滤”。而且我不知道过滤器功能，尽管这是我想确定的。这些峰很少在单个事件中达成一致。此外，该系统有一个大约 180 Hz 的外部应用不稳定背景信号，在视觉上很容易区分，但变化太大而无法轻松去除。
从固定长度的重叠窗口开始，比如 300 长，每 100 个样本开始。会有一些冗余。您的其他问题属于他们自己单独的 DSP 问题。

【解决方案3】：

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

将绘图转换为图像，并将其粘贴到像 gimp 或 Photoshop 这样的体面的图像程序中
只要有间隙，就将图分成离散的图像
将第一组图放在一条水平线上
将第二个系列放在其正下方的水平线上
直观地识别第一个相关事件
如果两个事件没有垂直排列：
- 选择该行左侧的任何实例以及右侧的所有实例
- 将这些东西向右拖动，直到它们对齐

这几乎就是音频编辑器的工作方式，因此，如果您将其转换为简单的音频格式，例如未压缩的 WAV 文件，您可以直接使用 Audacity 之类的工具对其进行操作。（当然，这听起来很可怕，但您可以很容易地移动数据图。）

实际上，Audacity 也有一种称为 nyquist 的脚本语言，所以如果您不需要该程序检测相关性（或者您至少愿意暂时推迟该步骤是）您可能会使用大胆的标记和 nyquist 的某种组合来自动对齐并在标记相关点后以您选择的格式导出干净的数据。

【讨论】：

这通常是一个有趣的想法，但它增加了分析信号的复杂性。大多数走这条路的人通常会执行转换为功率谱密度并比较 PSD 图像 - 通常即使在维度增加的情况下域变化也会有净增益。
我只提出使用图像作为思想实验，但我认为大胆实际上是一个非常明智的选择。这正是音频编辑程序旨在解决的问题，当您添加 nyquist（基本上是音频/时间序列数据的方案）时，这是一个非常强大的组合。
我实际上正在使用 pyplot 和我发布的另一个问题的答案来做这样的事情，该问题允许我将绘图缩放在一起但单独平移它们。所以这就是我如何知道数据集中的第一个常见事件的时钟差约为 4.78 秒。

【解决方案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)

【讨论】：

我已经在 pyplot 中做了类似的事情。问题在于，对于数百个事件，手动执行该过程很痛苦。

【解决方案5】：

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

我会研究 SciPy 的双变量 optimize 函数。我会使用mask 或暂时将“间隙”上的数据（例如 Ax' 和 By）归零（假设可以通过编程方式确定间隙）。

为提高流程效率，请从 A 和 B 的粗略抽样开始，但在fmin（或您选择的任何优化器）中设置与您的抽样相称的精度。然后继续对整个数据集进行更精细的采样窗口，直到您的窗口变窄并且没有下采样。

编辑 - 匹配轴

关于尝试确定哪个轴与给定轴共线的问题，并且不知道数据的特征，我可以指出一个类似的问题。查看pHash 或this 帖子中列出的任何其他方法，以帮助识别类似的波形。

【讨论】：

我会研究您推荐的功能。如果可能的话，我还想自动确定跨设备的正确轴配对。
@BobC：听起来你知道 Axe 与 By 大致同轴？我理解错了吗？
@Paul：是的，加速度计安装时轴大致平行，可能在几度以内。但是 Ax 与 By 对齐，因此加速度计彼此之间有 2 个 90 度旋转。
@BobC：那我不知道你说的“我也想自动确定正确的轴配对”是什么意思
@Paul：给定来自一对未知方向但轴平行的 3 轴设备的数据，根据数据内容对轴进行配对。