我正在做一些研究，但一直停留在分析阶段（应该更多地关注我的统计讲座）。
我收集了两个同时的信号：体积和胸部扩张变化的综合流速。我想比较信号并最终希望从胸部扩张信号中得出体积。但首先我必须对齐/同步我的数据。
由于记录不是在同一时间开始的，并且胸部扩张被捕获的时间更长，我需要在胸部扩张数据集中找到与我的体积数据相对应的数据，并衡量它们的对齐程度。如果两个信号不是在完全相同的时间开始，或者在不同比例和不同分辨率的数据之间开始，我不太确定如何处理。
我附上了两个信号的示例（https://docs.google.com/spreadsheet/ccc?key=0As4oZTKp4RZ3dFRKaktYWEhZLXlFbFVKNmllbGVXNHc），如果我可以提供任何进一步的信息，请告诉我。
该问题询问当以规则但不同的间隔对序列进行采样时，如何找到一个时间序列（“扩展”）滞后于另一个（“体积”）的量。
在这种情况下，两个系列都表现出合理的连续行为，如图所示。这意味着 (1) 可能需要很少或不需要初始平滑，并且 (2) 重新采样可以像线性或二次插值一样简单。由于平滑度，二次方可能会稍微好一些。 重采样后，通过最大化互相关找到滞后，如线程所示，对于两个偏移采样的数据序列，它们之间的偏移量的最佳估计是多少？.
为了说明，我们可以使用问题中提供的数据，R用于伪代码。让我们从基本功能、互相关和重采样开始：
cor.cross <- function(x0, y0, i=0) {
#
# Sample autocorrelation at (integral) lag `i`:
# Positive `i` compares future values of `x` to present values of `y`';
# negative `i` compares past values of `x` to present values of `y`.
#
if (i < 0) {x<-y0; y<-x0; i<- -i}
else {x<-x0; y<-y0}
n <- length(x)
cor(x[(i+1):n], y[1:(n-i)], use="complete.obs")
}
这是一个粗略的算法：基于 FFT 的计算会更快。但是对于这些数据（涉及大约 4000 个值）来说已经足够了。
resample <- function(x,t) {
#
# Resample time series `x`, assumed to have unit time intervals, at time `t`.
# Uses quadratic interpolation.
#
n <- length(x)
if (n < 3) stop("First argument to resample is too short; need 3 elements.")
i <- median(c(2, floor(t+1/2), n-1)) # Clamp `i` to the range 2..n-1
u <- t-i
x[i-1]*u*(u-1)/2 - x[i]*(u+1)*(u-1) + x[i+1]*u*(u+1)/2
}
我将数据下载为逗号分隔的 CSV 文件并去掉了它的标题。（标题对 R 造成了一些我不想诊断的问题。）
data <- read.table("f:/temp/a.csv", header=FALSE, sep=",",
col.names=c("Sample","Time32Hz","Expansion","Time100Hz","Volume"))
注意 ：此解决方案假定每个系列的数据都按时间顺序排列，其中任何一个都没有间隙。 这允许它使用值中的索引作为时间的代理，并通过时间采样频率缩放这些索引以将它们转换为时间。
事实证明，这些仪器中的一种或两种会随着时间的推移而略有漂移。在继续之前消除这些趋势是很好的。此外，由于最后音量信号逐渐变细，我们应该将其剪掉。
n.clip <- 350 # Number of terminal volume values to eliminate
n <- length(data$Volume) - n.clip
indexes <- 1:n
v <- residuals(lm(data$Volume[indexes] ~ indexes))
expansion <- residuals(lm(data$Expansion[indexes] ~ indexes)
我重新采样频率较低的系列，以便从结果中获得最精确的结果。
e.frequency <- 32 # Herz
v.frequency <- 100 # Herz
e <- sapply(1:length(v), function(t) resample(expansion, e.frequency*t/v.frequency))
现在可以计算互相关——为了提高效率，我们只搜索一个合理的滞后窗口——并且可以识别找到最大值的滞后。
lag.max <- 5 # Seconds
lag.min <- -2 # Seconds (use 0 if expansion must lag volume)
time.range <- (lag.min*v.frequency):(lag.max*v.frequency)
data.cor <- sapply(time.range, function(i) cor.cross(e, v, i))
i <- time.range[which.max(data.cor)]
print(paste("Expansion lags volume by", i / v.frequency, "seconds."))
输出告诉我们扩展滞后了 1.85 秒。（如果最后 3.5 秒的数据没有被剪裁，输出将是 1.84 秒。）
以多种方式检查所有内容是个好主意，最好是目视检查。一、互相关函数：
plot(time.range * (1/v.frequency), data.cor, type="l", lwd=2,
xlab="Lag (seconds)", ylab="Correlation")
points(i * (1/v.frequency), max(data.cor), col="Red", cex=2.5)
接下来，让我们及时注册这两个系列，并将它们一起绘制在相同的轴上。
normalize <- function(x) {
#
# Normalize vector `x` to the range 0..1.
#
x.max <- max(x); x.min <- min(x); dx <- x.max - x.min
if (dx==0) dx <- 1
(x-x.min) / dx
}
times <- (1:(n-i))* (1/v.frequency)
plot(times, normalize(e)[(i+1):n], type="l", lwd=2,
xlab="Time of volume measurement, seconds", ylab="Normalized values (volume is red)")
lines(times, normalize(v)[1:(n-i)], col="Red", lwd=2)
看起来还不错！ 不过，我们可以通过散点图更好地了解配准质量。我按时间改变颜色以显示进度。
colors <- hsv(1:(n-i)/(n-i+1), .8, .8)
plot(e[(i+1):n], v[1:(n-i)], col=colors, cex = 0.7,
xlab="Expansion (lagged)", ylab="Volume")
我们正在寻找沿着一条线来回追踪的点：从这条线的变化反映了膨胀对体积的时滞响应中的非线性。虽然有一些变化，但它们非常小。然而，这些变化如何随时间变化可能具有一定的生理意义。统计学的美妙之处，尤其是它的探索性和视觉方面，是它如何倾向于创造好的问题和想法以及有用的答案。