数据集有11719 个点。可以像这样检测连续段
d = load('summerRunningMean.csv');
% detect continuous parts
s = sqrt(sum(diff(d)' .^ 2));
ind = [0, find(s > 100 * mean(s)), size(d, 1)];
这导致 6 个段的长度从 534 到 4255 点。
绘制在空间中叠加的轨迹
for i = 1 : 6
x = d(ind(i) + 1 : ind(i + 1), 1);
y = d(ind(i) + 1 : ind(i + 1), 2);
plot(x, y)
hold all
end
xlabel x
ylabel y
axis equal
给出这个结果
并显示轨道具有大致相同的形状，并且除了其中一个之外，所有轨道都从同一个位置开始。
然而，时间序列的叠加表明沿该形状的运动以不同的速度发生：
因此，对这些数据进行统计的问题是没有共同的参考点来将来自不同轨道的数据点相互分配。
在下文中，我排除了轨道#6，因为它的起点与其他轨道不同。
定义这种公共参考的一种方法是计算沿轨道的曲线长度：
l = [0 ; cumsum(sqrt(diff(x) .^ 2 + diff(y) .^ 2))];
绘制这个长度而不是时间使得 x 和 y 坐标更相似，因此具有可比性：
根据这个结果，可以计算平均值和其他统计数据。要真正做到这一点，数据必须通过插值重新参数化：
li = 0 :0.001: l(end);
xi = interp1(l, x, li);
yi = interp1(l, y, li);
现在，对于每个轨道，我们都有一个共同的参考，并且可以将转换后的数据存储在共同的数据矩阵中：
n = numel(li);
xp(1 : n, i) = xi;
yp(1 : n, i) = yi;
i轨道索引在哪里。矩阵xp和yp必须初始化为 NaN，因为轨道具有不同的长度。然后，可以计算统计数据，例如平均值：
xm = nanmean(xp, 2);
ym = nanmean(yp, 2);
生成的平均轨迹与原始轨迹一起：
为了进一步提高轨道之间的一致性，可以平滑它们以减少随机变化（span = 5000似乎效果很好）：
xs = smooth(xi, span);
ys = smooth(yi, span);
之后，必须重新计算曲线长度参数：
l = [0 ; cumsum(sqrt(diff(xs) .^ 2 + diff(ys) .^ 2))];
结果：
共同参考所依据的协议明显得到了改进。数据必须再次重新参数化
li = 0 :0.001: l(end);
xsi = interp1(l, xs, li);
ysi = interp1(l, ys, li);
存储在公共数据矩阵中的结果
n = numel(li);
xp(1 : n, i) = xsi;
yp(1 : n, i) = ysi;
和计算的平均值
xm = nanmean(xp, 2);
ym = nanmean(yp, 2);
生成的平均轨迹与原始轨迹一起：
结果看起来更加平滑。然而，平滑减少了两个轨道末端的大环的周长，因此“平均轨道”不再位于原始轨道之间。平滑度和接近原始轨迹之间的权衡可以通过 的值来调节span。
生成最后一个数字的完整代码包含在这里；为了重现未平滑的版本，设置span = 1;.
d = load('summerRunningMean.csv');
% detect tracks as continuous segments
s = sqrt(sum(diff(d)' .^ 2));
ind = [0, find(s > 100 * mean(s)), size(d, 1)];
% remove 6th track because itis an outlier
ind = ind(1 : end - 1);
N = size(ind, 2) - 1;
% smoothing parameter
span = 5000;
xp = nan(13000, N); % I know, hardcoded, not nice.
yp = nan(13000, N);
for i = 1 : N
% extract data
x = d(ind(i) + 1 : ind(i + 1), 1);
y = d(ind(i) + 1 : ind(i + 1), 2);
% determine length along curve
% to use as curve parameter
l = [0 ; cumsum(sqrt(diff(x) .^ 2 + diff(y) .^ 2))];
% reparametrize by interpolation
li = 0 :0.001: l(end);
xi = interp1(l, x, li);
yi = interp1(l, y, li);
% smooth to remove small deviations
xs = smooth(xi, span);
ys = smooth(yi, span);
% determine length along smoothed curve
% as improved curve parameter
l = [0 ; cumsum(sqrt(diff(xs) .^ 2 + diff(ys) .^ 2))];
% again, reparametrize by interpolation
li = 0 :0.001: l(end);
xsi = interp1(l, xs, li);
ysi = interp1(l, ys, li);
% store
n = numel(li);
xp(1 : n, i) = xsi;
yp(1 : n, i) = ysi;
plot(x, y)
axis equal
hold all
end
% compute mean
xm = nanmean(xp, 2);
ym = nanmean(yp, 2);
% plot mean
plot(xm, ym, 'k', 'LineWidth', 2)
xlabel x
ylabel y