BML_TIMEWARP aligns and linearly warps slave to master Use as warpedcoords = bml_timealign(master, slave) warpedcoords = bml_timealign(cfg, master, slave) This function first aligns slave to master by calling BML_TIMEALIGN and then applies a constrained linear time-warp to slave, in order to maximize the correlation with master. The time-warping is done with the raws pre-processed by BML_TIMEALIGN. The time-warping is done by optimizing the cost function below using the 'simplex' method w(t) = wt0 + pivot_time + (t - pivot_time) ws1 s = f(w(t)) cost = -s*p/sqrt(s*s p*p) + (wt0/penalty_wt0)^2 + ((ws1-1)/penalty_ws1)^4 where pivot_time is the midpoint of the intersection of master and slave after initial alignment, f(t) is a interpolation function constructed from the slave's data and time points (after alingment) and '*' is the dot product . wt0 and ws1 are the time shifting and stretching paramenters, respectiveley. penalty_wt0 and penalty_ws1 are used to constrain their value by adding a penalty cost. penalty_wt0 is set to the max between penalty_wt0_min and abs(ws1-1)*duration to avoid the optimization algorithm to missalign the signals in the first iterations. Parameters ---------- cfg is a configuration structure see BML_TIMEALIGN parameters cfg.penalty_wt0_min - double: see details above. Defaults to 1e-6. cfg.penalty_ws1 - double: see details above. Defaults to 1e-3. cfg.timewarp - logical: defaults to true. If false no warping is done master - FT_DATATYPE_RAW continuous with single channel and trial slave - FT_DATATYPE_RAW continuous with single channel and trial Returns ------- warpedcoords structure with fields warpedcoords.s1 - integer: first coordinate sample warpedcoords.t1 - double: first coordinate time in seconds warpedcoords.s2 - integer: second coordinate sample warpedcoords.t2 - double: second coordinate time in seconds warpedcoords.wt0 - double: fitted parameter warpedcoords.ws1 - double: fitted parameter
0001 function warpedcoords = bml_timewarp(cfg, master, slave) 0002 0003 % BML_TIMEWARP aligns and linearly warps slave to master 0004 % 0005 % Use as 0006 % warpedcoords = bml_timealign(master, slave) 0007 % warpedcoords = bml_timealign(cfg, master, slave) 0008 % 0009 % This function first aligns slave to master by calling BML_TIMEALIGN and 0010 % then applies a constrained linear time-warp to slave, in order to 0011 % maximize the correlation with master. The time-warping is done with the 0012 % raws pre-processed by BML_TIMEALIGN. 0013 % 0014 % The time-warping is done by optimizing the cost function below using 0015 % the 'simplex' method 0016 % 0017 % w(t) = wt0 + pivot_time + (t - pivot_time) ws1 0018 % s = f(w(t)) 0019 % cost = -s*p/sqrt(s*s p*p) + (wt0/penalty_wt0)^2 + ((ws1-1)/penalty_ws1)^4 0020 % 0021 % where pivot_time is the midpoint of the intersection of master and 0022 % slave after initial alignment, f(t) is a interpolation function 0023 % constructed from the slave's data and time points (after alingment) and 0024 % '*' is the dot product . 0025 % wt0 and ws1 are the time shifting and stretching paramenters, 0026 % respectiveley. penalty_wt0 and penalty_ws1 are used to constrain their 0027 % value by adding a penalty cost. penalty_wt0 is set to the max between 0028 % penalty_wt0_min and abs(ws1-1)*duration to avoid the optimization 0029 % algorithm to missalign the signals in the first iterations. 0030 % 0031 % Parameters 0032 % ---------- 0033 % cfg is a configuration structure 0034 % see BML_TIMEALIGN parameters 0035 % cfg.penalty_wt0_min - double: see details above. Defaults to 1e-6. 0036 % cfg.penalty_ws1 - double: see details above. Defaults to 1e-3. 0037 % cfg.timewarp - logical: defaults to true. If false no warping is done 0038 % 0039 % master - FT_DATATYPE_RAW continuous with single channel and trial 0040 % slave - FT_DATATYPE_RAW continuous with single channel and trial 0041 % 0042 % Returns 0043 % ------- 0044 % warpedcoords structure with fields 0045 % warpedcoords.s1 - integer: first coordinate sample 0046 % warpedcoords.t1 - double: first coordinate time in seconds 0047 % warpedcoords.s2 - integer: second coordinate sample 0048 % warpedcoords.t2 - double: second coordinate time in seconds 0049 % warpedcoords.wt0 - double: fitted parameter 0050 % warpedcoords.ws1 - double: fitted parameter 0051 0052 if nargin == 2 0053 slave = master; 0054 master = cfg; 0055 cfg = []; 0056 elseif nargin ~=3 0057 error('incorrect number of arguments'); 0058 end 0059 0060 penalty_wt0_min = bml_getopt(cfg,'penalty_wt0_min', 60); 0061 penalty_ws1 = bml_getopt(cfg,'penalty_ws1', 1e-3); 0062 timewarp = bml_getopt(cfg,'timewarp', true); 0063 0064 mc=[]; sc=[]; %original master and slave time coordinates 0065 mc.s1=1; mc.s2=length(master.time{1}); 0066 mc.t1=master.time{1}(1); mc.t2=master.time{1}(end); 0067 sc.s1=1; sc.s2=length(slave.time{1}); 0068 sc.t1=slave.time{1}(1); sc.t2=slave.time{1}(end); 0069 0070 %aligning and transforming raws 0071 [slave_delta_t, ~, master, slave] = bml_timealign(cfg, master, slave); 0072 sc.t1 = sc.t1 + slave_delta_t; sc.t2 = sc.t2 + slave_delta_t; 0073 0074 %calculating pivot time 0075 ovlp = [max(mc.t1,sc.t1),min(sc.t2,sc.t2)]; 0076 pivot_time = mean(ovlp); 0077 0078 %saving cropped start and end indices of original slave raw 0079 crop_sc=[]; %cropped slave coordinates 0080 [crop_sc.s1, crop_sc.s2]=bml_crop_idx(sc,ovlp(1),ovlp(2)); 0081 crop_sc.t1=bml_idx2time(sc,crop_sc.s1); 0082 crop_sc.t2=bml_idx2time(sc,crop_sc.s2); 0083 0084 %cropping raws to overlap region 0085 master = bml_crop(master, ovlp(1), ovlp(2)); 0086 slave = bml_crop(slave, ovlp(1), ovlp(2)); 0087 ovlp_duration = ovlp(2) - ovlp(1); 0088 0089 %linear warping - wt0, ws1 0090 f = @(t) interp1(slave.time{1}(1,:),slave.trial{1}(1,:),t,'PCHIP',0); 0091 p = master.trial{1}(1,:); 0092 t = master.time{1}(1,:); 0093 s0 = f(t); 0094 dot0 = sqrt(dot(s0,s0) * dot(p,p)); 0095 0096 function cost = costfun(param) % [wt0, ws1] 0097 s = f(param(1) + pivot_time + (t - pivot_time) .* param(2)); 0098 penalty_wt0_dur = max([penalty_wt0_min, ovlp_duration*abs(param(2)-1)]); 0099 cost = -dot(s,p)/dot0 + (param(1)/penalty_wt0_dur)^2 + ((param(2)-1)/penalty_ws1)^4; 0100 end 0101 0102 wt0=0; ws1=1; 0103 if timewarp 0104 fited=fminsearch(@costfun,[wt0,ws1]); 0105 wt0=fited(1); ws1=fited(2); 0106 end 0107 0108 warpedcoords=[]; 0109 warpedcoords.t1 = pivot_time - wt0 - (1/ws1)*(crop_sc.t2 - crop_sc.t1)/2; 0110 warpedcoords.t2 = pivot_time - wt0 + (1/ws1)*(crop_sc.t2 - crop_sc.t1)/2; 0111 warpedcoords.s1 = crop_sc.s1; 0112 warpedcoords.s2 = crop_sc.s2; 0113 warpedcoords.wt0 = wt0; 0114 warpedcoords.ws1 = ws1; 0115 0116 end