多次重复实验得到的时序数据的对齐操作

一、问题的提出

进行多次重复实验会得到多次实验数据，如同世界上没有两片同样的树叶一样，实验过程所需的时间两两之间也不可能完全一样。如何将多次重复实验每个trial进行对齐是进行数据统计的关键步骤之一。本文介绍了一篇Neuron期刊文章上的方法。

上图所示是文章中多次重复实验和数据对齐的实现过程。图A进行正转，B进行反转，分别都进行了7次（指一次trial里转7次，从图片中看出这个示意图中将一个session里包含了18次trial）。图C是某个神经元发放的情形，图D是某块肌肉的情形，图E、F是scaling之后的，图G是图E的均值，曲线是SEM形式，背景颜色表示垂直方向上的速度（即亮表示速度曲线的顶点，暗表示速度曲线的低点）上方的小竖线表示一次转完成。

二、重复实验数据对齐实战

如何实现从图C到图E（或图D到图F）
原文如下：
To ensure response features were not lost to misalignment, we developed a technique to adaptively align trials within a condition. First, trials were aligned on movement onset. Individual trials were then scaled so that all trials had the same duration (set to be the median duration across trials). Because monkeys usually cycled at a consistent speed (within a given condition) this brought trials largely into alignment: e.g., the top of each cycle occurred at nearly the same time for each trial. The adaptive alignment procedure was used to correct any remaining slight misalignments. The time-base for each trial was scaled so that the position trace on that trial closely matched the average position of all trials. This involved a slight non-uniform stretching, and resulted in the timing of all key moments – such as when the hand passed the top of the cycle – being nearly identical across trials. This ensured that high-frequency temporal response features (e.g., the small peak in Figure 1G) were not lost to averaging.

简单讲，就是先对齐movement onset，然后计算duration的均值（这里是中间值），最后把所有trial都拉伸或压缩到均值。当然，在这之前要删除那些偏离太多的trial（这样的trial不能太多，只能是少数，为了达到这样的效果，动物训练要训好）。以下以模拟数据具体实现。

（1）脉冲电信号的对齐实战

计算均值就不说了，直接讲如何进行scaling：
主要的困惑是spike raster是怎么进行的？因为连续信号好弄，直接重采样，matlab里有resample命令，但是没弄过spike信号。

一步一步来：

%先生成一个raster，
X = randn(1,1050);
X(X>0)=0;
X(X<0)=1;
[~,index_spike] = find(X==1);
line([index_spike;index_spike], [0;0.5],'color','k')

得到：

如何把1050个点scaling到1000个点呢？
（大家最好看看index_spike里的值代表什么含义，没错，就是timestamp，而不是1）
思路：如果在1050的时刻点有个spike，那么scaling之后应该是在1000的时刻点有个与之对应的spike，在原来1000的时刻点有个spike，那么与之对应scaling之后应该是（1000/1050）*1000时刻点，即952.381有个spike，但是采样率是整数，所以近似952时刻点有个spike。看出来了吧，就是（1000/1050）乘以origin时刻点，然后取整，就得到scaling之后的时刻点，就能形成一个新的scaling的raster，我们来实现一下：

scaling_index_spike = (1000/1050)*index_spike;
int_scaling_index_spike = round(scaling_index_spike );
line([int_scaling_index_spike;int_scaling_index_spike], [0;0.5],'color','k')

把scaling前后的raster画在一起：

figure
line([index_spike;index_spike], [0;0.5],'color','k')
figure
line([int_scaling_index_spike;int_scaling_index_spike], [0;0.5],'color','k')

得到结果：

不要去看横坐标，只看raster图是不是很像啊，完成～

（2）连续信号的对齐实现

%首先生成两个正弦
Y1 = sin(0:0.01:5*pi);
Y2 = sin(0:0.0102:5*pi);
figure;
subplot(211)
plot(1:length(Y1),Y1)
subplot(212)
plot(1:length(Y2),Y2)

得到结果：

从结果可知，Y1由1571个数组成，Y2由1540个数组成
假设那么我们要将二者都scaling到1500个数，怎么操作呢？思路跟上面是一样的，只不过当index相同时，进行均值操作一下。

scaling_index_Y1 = (1500/1571)*(1:1571); %将1571个数与1500个数对应
int_scaling_index_Y1 = round(scaling_index_Y1);  %把小数点去掉
int_scaling_index_Y1 = unique(round(scaling_index_Y1));  %这步没啥特别，实质就是1:1500

for i = 1:length(unique_int_scaling_index_Y1)
    [m n] = find(int_scaling_index_Y1==unique_int_scaling_index_Y1(i));
    scaling_Y1(unique_int_scaling_index_Y1(i)) = mean(Y1(n));
end
%最后画出来瞅瞅
figure
plot(1:length(scaling_Y1),scaling_Y1)

上面完成了Y1的scaling，为了日后方便，我们把上面的操作包装成function：

function scaling_value = scaling_Y(Y, origin_index, scaling_index)
scaling_index_Y1 = (scaling_index/origin_index)*(1:origin_index); %将1571个数与1500个数对应
int_scaling_index_Y1 = round(scaling_index_Y1);  %把小数点去掉
unique_int_scaling_index_Y1 = unique(round(scaling_index_Y1));  %这步没啥特别，实质就是1:1500

for i = 1:length(unique_int_scaling_index_Y1)
    [m n] = find(int_scaling_index_Y1==unique_int_scaling_index_Y1(i));
    scaling_value(unique_int_scaling_index_Y1(i)) = mean(Y(n));
end
%最后画出来瞅瞅
figure
plot(1:length(scaling_value),scaling_value)

调用计算Y2的scaling：

scaling_Y2 = scaling_Y(Y2, length(Y2), 1500)

得到结果：

最后，为了直观，我们把Y1，Y2和scaling_Y1, scaling_Y2画在一起：

figure
subplot(511)
plot(1:length(Y1),Y1)
subplot(512)
plot(1:length(Y2),Y2)
subplot(513)
plot(1:length(scaling_Y1),scaling_Y1)
%set(gca,'xtick',0:1600:1600);
subplot(514)
plot(1:length(scaling_Y2),scaling_Y2)
%set(gca,'xtick',0:1600:1600);
subplot(515)
plot(1:length(Y1),Y1)
hold on
plot(1:length(Y2),Y2)
plot(1:length(scaling_Y1),scaling_Y1)
plot(1:length(scaling_Y2),scaling_Y2)

得到结果：