bml_timealign

 BML_TIMEALIGN aligns slave to master and returns the slave's delta t

 Use as
   slave_delta_t = bml_timealign(master, slave)
   slave_delta_t = bml_timealign(cfg, master, slave)
   [slave_delta_t, max_corr] = bml_timealign(cfg, master, slave)
   [slave_delta_t, max_corr, master, slave] = bml_timealign(cfg, master, slave)


 cfg is a configuration structure with fields:
 cfg.resample_freq - double: frequency to resample and aligned master and
            slave (Hz). Defaults to 10000.
 cfg.method - string: method use for preprocessing master and slave
            'env' or 'envelope' (default) - see BML_ENVELOPE_BINABS
            'lpf' or 'low-pass-filter' - see ft_preproc_lowpassfilter
 cfg.env_freq - double: frequency of the envelope. Defaults to 100Hz
 cfg.lpf_freq - double: cut-frequency of the low-pass filter. Defaults to
            4000Hz. 
 cfg.scan - double: time window in which to scan for a maximal correlation
            if a scalar is given the time window is [-scan, scan]
            if a length 2 vector is given it should be [-a, b], where 'a'
            and 'b' are positive numbers in seconds. 
 cfg.freqreltol - double: frequency relative tolerance. defaults to 1e-5
 cfg.penalty_tau - double: penalty time use to weight the correlation.
            This allows to bound slave_delta_t (as with cfg.scan) but in
            a continuous way. If empty (default) no penalty is applied.
 cfg.penalty_n - double: penalty 'hill-coefficient' use to weight the
            correlation. Defines how abrupt is the penalty increase when
            slave_delta_t > cfg.penalty_tau. Defaults to 2. 
 cfg.ft_feedback - string: default to 'no'. Defines verbosity of fieldtrip
            functions 

 master - FT_DATATYPE_RAW continuous with single channel and trial
 slave - FT_DATATYPE_RAW continuous with single channel and trial

 returns
 slave_delta_t - double: time in seconds by which to shift the slave to 
           align it to master
 max_corr - double: maximum correlation coefficient achieved for the shift
           slave_delta_t
 master - FT_DATATYPE_RAW: master raw after applying the preprocessing
 slave - FT_DATATYPE_RAW: slave raw after applying the preprocessing

0001 function [slave_delta_t, max_corr, master, slave] = bml_timealign(cfg, master, slave)
0002 
0003   % BML_TIMEALIGN aligns slave to master and returns the slave's delta t
0004   %
0005   % Use as
0006   %   slave_delta_t = bml_timealign(master, slave)
0007   %   slave_delta_t = bml_timealign(cfg, master, slave)
0008   %   [slave_delta_t, max_corr] = bml_timealign(cfg, master, slave)
0009   %   [slave_delta_t, max_corr, master, slave] = bml_timealign(cfg, master, slave)
0010   %
0011   %
0012   % cfg is a configuration structure with fields:
0013   % cfg.resample_freq - double: frequency to resample and aligned master and
0014   %            slave (Hz). Defaults to 10000.
0015   % cfg.method - string: method use for preprocessing master and slave
0016   %            'env' or 'envelope' (default) - see BML_ENVELOPE_BINABS
0017   %            'lpf' or 'low-pass-filter' - see ft_preproc_lowpassfilter
0018   % cfg.env_freq - double: frequency of the envelope. Defaults to 100Hz
0019   % cfg.lpf_freq - double: cut-frequency of the low-pass filter. Defaults to
0020   %            4000Hz.
0021   % cfg.scan - double: time window in which to scan for a maximal correlation
0022   %            if a scalar is given the time window is [-scan, scan]
0023   %            if a length 2 vector is given it should be [-a, b], where 'a'
0024   %            and 'b' are positive numbers in seconds.
0025   % cfg.freqreltol - double: frequency relative tolerance. defaults to 1e-5
0026   % cfg.penalty_tau - double: penalty time use to weight the correlation.
0027   %            This allows to bound slave_delta_t (as with cfg.scan) but in
0028   %            a continuous way. If empty (default) no penalty is applied.
0029   % cfg.penalty_n - double: penalty 'hill-coefficient' use to weight the
0030   %            correlation. Defines how abrupt is the penalty increase when
0031   %            slave_delta_t > cfg.penalty_tau. Defaults to 2.
0032   % cfg.ft_feedback - string: default to 'no'. Defines verbosity of fieldtrip
0033   %            functions
0034   %
0035   % master - FT_DATATYPE_RAW continuous with single channel and trial
0036   % slave - FT_DATATYPE_RAW continuous with single channel and trial
0037   %
0038   % returns
0039   % slave_delta_t - double: time in seconds by which to shift the slave to
0040   %           align it to master
0041   % max_corr - double: maximum correlation coefficient achieved for the shift
0042   %           slave_delta_t
0043   % master - FT_DATATYPE_RAW: master raw after applying the preprocessing
0044   % slave - FT_DATATYPE_RAW: slave raw after applying the preprocessing
0045   
0046   if nargin == 2
0047     slave = master;
0048     master = cfg;
0049     cfg = [];
0050   elseif nargin ~=3
0051     error('incorrect number of arguments in call');
0052   end
0053   
0054   resample_freq     = bml_getopt(cfg,'resample_freq', 10000);
0055   scan              = bml_getopt(cfg, 'scan');
0056   freqreltol        = bml_getopt(cfg, 'freqreltol', 1e-5);
0057   method            = string(bml_getopt(cfg, 'method', 'envelope'));
0058   env_freq          = bml_getopt(cfg,'env_freq', 100);
0059   lpf_freq          = bml_getopt(cfg,'lpf_freq', 4000);
0060   penalty_tau       = bml_getopt(cfg,'penalty_tau');
0061   penalty_n         = bml_getopt(cfg,'penalty_n', 2); 
0062   ft_feedback       = bml_getopt(cfg,'ft_feedback','no');
0063   ft_feedback       = ft_feedback{1};
0064 
0065   %assert single trial and channel
0066   if numel(master.trial) > 1; error('master should be single trial raw'); end
0067   if numel(slave.trial) > 1; error('slave should be single trial raw'); end
0068   if numel(master.label) > 1; error('master should be single channel raw'); end
0069   if numel(slave.label) > 1; error('slave should be single channel raw'); end
0070 
0071   %calculating scan range
0072   mc=[]; sc=[]; %master and slave time coordinates
0073   mc.s1=1; mc.s2=length(master.time{1});
0074   mc.t1=master.time{1}(1); mc.t2=master.time{1}(end);
0075   sc.s1=1; sc.s2=length(slave.time{1});
0076   sc.t1=slave.time{1}(1); sc.t2=slave.time{1}(end);
0077   max_scan_range = [mc.t1 - sc.t2, mc.t2 - sc.t1];
0078   if prod(max_scan_range) > 0 % if files do not overlap
0079     slave_delta_t=nan;
0080     max_corr=nan;
0081     return
0082   end
0083   if isempty(scan)
0084     scan = max_scan_range;
0085   elseif length(scan)==1
0086     scan=[max(-scan,max_scan_range(1)), min(scan,max_scan_range(2))];
0087   elseif length(scan)==2
0088     scan=[max(scan(1),max_scan_range(1)), min(scan(2),max_scan_range(2))];  
0089   else
0090     error('invalid use of cfg.scan argument');
0091   end
0092 
0093   %robust estimation of mean and std
0094   master_median = median(master.trial{1});
0095   slave_median = median(slave.trial{1});
0096   master_std = robust_std(master.trial{1});
0097   slave_std = robust_std(slave.trial{1});
0098   if master_std==0; error('master can''t be correlated'); end
0099   if slave_std==0; error('slave can''t be correlated'); end
0100 
0101   %cropping and padding to correlation time window
0102   ctw = [sc.t1+scan(1), sc.t2+scan(2)];
0103   master = bml_pad(master, ctw(1), ctw(2), 0);
0104   slave = bml_pad(slave, master.time{1}(1), master.time{1}(end), 0);
0105 
0106   %common resample frequency
0107   cfg=[]; cfg.feedback=ft_feedback;
0108   cfg.resamplefs=resample_freq;
0109   master = ft_resampledata(cfg, master);
0110   cfg=[]; cfg.feedback=ft_feedback;
0111   cfg.time=master.time; cfg.method='linear';
0112   slave = ft_resampledata(cfg, slave);
0113 
0114   %checking slave resampling
0115   if abs(slave.fsample/master.fsample-1) < freqreltol
0116     slave.fsample = master.fsample;
0117   else
0118     error('failed to resample slave to master''s time');
0119   end
0120   is_nan=isnan(slave.trial{1});
0121   if sum(is_nan)>0
0122     master.trial{1} = master.trial{1}(:,~is_nan);
0123     master.time{1} = master.time{1}(:,~is_nan);
0124     slave.trial{1} = slave.trial{1}(:,~is_nan);
0125     slave.time{1} = slave.time{1}(:,~is_nan);
0126   end
0127 
0128   %methods
0129   if ismember(lower(method),{'env','envelope'})  %envelope correlation
0130     cfg=[]; cfg.freq = env_freq; %calculating envelops
0131     master = bml_envelope_binabs(cfg,master);
0132     slave = bml_envelope_binabs(cfg,slave);
0133     try_polarity=false;
0134 
0135   elseif ismember(lower(method),{'lpf','low-pass-filter'})  %low-pass-filter
0136     master.trial{1} = ft_preproc_lowpassfilter(master.trial{1},...
0137                       master.fsample, lpf_freq, 4, 'but', 'twopass');
0138     slave.trial{1} = ft_preproc_lowpassfilter(slave.trial{1},...
0139                       slave.fsample, lpf_freq, 4, 'but', 'twopass');              
0140     try_polarity=true;                
0141 
0142   else
0143     error('unknown method');
0144   end
0145 
0146   %normalizing data
0147   master.trial{1} = (master.trial{1} - master_median) / master_std;
0148   slave.trial{1} = (slave.trial{1} - slave_median) / slave_std;
0149 
0150   %correlation
0151   [corr,lag]=xcorr(master.trial{1}(1,:), slave.trial{1}(1,:),...
0152                     floor(max(abs(scan))*master.fsample),'coeff');
0153   [max_corr_idx,max_corr] = find_delta_corr(corr,lag,penalty_tau,penalty_n);  
0154 
0155   if try_polarity                              
0156     [corr_m,lag_m]=xcorr(master.trial{1}(1,:), (-1).* slave.trial{1}(1,:),...
0157                       floor(max(abs(scan))*master.fsample),'coeff');                 
0158     [max_corr_idx_m,max_corr_m] = find_delta_corr(corr_m,lag_m,penalty_tau,penalty_n);   
0159     if max_corr_m > max_corr
0160       slave.trial{1}=(-1).* slave.trial{1};
0161       max_corr_idx = max_corr_idx_m;
0162       max_corr = max_corr_m;
0163       lag=lag_m;
0164     end
0165   end
0166 
0167   slave_delta_t = lag(max_corr_idx)/master.fsample;
0168   slave.time{1} = slave.time{1} + slave_delta_t;
0169 
0170 end
0171 
0172 % private function
0173 function [max_corr_idx,max_corr] = find_delta_corr(corr,lag,penalty_tau,penalty_n)
0174   if ~isempty(penalty_tau)
0175     [~,max_corr_idx]=max(corr./((penalty_tau*master.fsample)^penalty_n + abs(lag).^penalty_n));
0176     max_corr=corr(max_corr_idx);
0177   else
0178     [max_corr,max_corr_idx]=max(corr);
0179   end
0180 
0181 end