function [trl, event] = trialfun_ssvep(cfg)

% TRIALFUN_SSVEP defines the segments timelocked to the cue
%
% cfg.trialdef.pre  = time before the first trigger
% cfg.trialdef.post = time after the last trigger

% 0 to 4.5 s for ssveps
% -1 s to +1 s for ERPs (maybe this can be used for tf analysis)
% defined below!

global fp

hdr = ft_read_header(cfg.dataset);
event = ft_read_event(cfg.dataset);

eventlist = [{event.type}',{event.sample}',{event.value}'];
eventlist(strcmp(eventlist(:,1),'STATUS')==false,:)=[];
eventlist(:,1)=[];
eventlist=cell2mat(eventlist);

% % output raw event codes and their frequency
% test=eventlist(:,2);
% test2(:,1)=unique(test);
% [test2(:,2), nc]=histc(test,unique(test));
% test2

% the de2bi function is in the Mathworks communication toolbox
% this provides a replacement for it (if needed)
ft_hastoolbox('comm', 1);

% decimal to binary
% (if not all 16 channels are used, then the higher ones
% are missing. So, initialize a matrix with 16 columns.
bina=zeros(size(eventlist,1),16);
tmp=de2bi(eventlist(:,2));
bina(:,1:size(tmp,2))=tmp;
evars={'Sample','Code','Evt','Evt2','Light1','Light2','Light3','Light4',...
  'AudioL','AudioR','Micro','Evnt0','Evnt1','Evnt2','TTL1','TTL2','TTL3', ...
  'TTL4','TTL5','TTL6'};
eventlist=[eventlist,nan(size(eventlist,1),2),bina];
etable=array2table(eventlist,'VariableNames',evars);
% Evt: final codes for events of interest

% %
% Lightsensor 1 Trig. 1 / 1 (LSB)
% Lightsensor 2 Trig. 2 / 2
% Lightsensor 3 Trig. 3 / 4
% Lightsensor 4 Trig. 4 / 8
% Audio, left channel Trig. 5 / 16
% Audio, right channel Trig. 6 / 32
% Microphone Trig. 7 / 64
% Event marker bit 0 Trig. 8 / 128
% Event marker bit 1 Trig. 9 / 1 (MSB)
% Event marker bit 2 Trig. 10 / 2
% TTL input, line 1 Trig. 11 / 4
% TTL input, line 2 Trig. 12 / 8
% TTL input, line 3 Trig. 13 / 16
% TTL input, line 4 Trig. 14 / 32
% TTL input, line 5 Trig. 15 / 64
% TTL input, line 6 Trig. 16 / 128

etable(:,6:14)=[];
% delete irrelevant channels

% Define conditions:
% ==================
% SSVEP trial onset (4.5 s)
% low load has TTL1, 2, 3, or 4 for Rings 1 to 4 (small to large)
% high load has TTL1, 2, 3, or 4 and TTL5+6 for Rings 1 to 4 (small to large)
% square-count: TTL5
% square-no count: TTL6
% end of SSVEP trial: all TTLs on
% correct resp: TTL1+2
% incorrect resp: TTL3+4
for t=1:11   % loop through the different conditions
  switch t
    case 1
      indx=(etable.TTL1==1 & etable.TTL2==0 & etable.TTL3==0 & etable.TTL4==0 & etable.TTL5==0 & etable.TTL6==0);
    case 2
      indx=(etable.TTL1==0 & etable.TTL2==1 & etable.TTL3==0 & etable.TTL4==0 & etable.TTL5==0 & etable.TTL6==0);
    case 3
      indx=(etable.TTL1==0 & etable.TTL2==0 & etable.TTL3==1 & etable.TTL4==0 & etable.TTL5==0 & etable.TTL6==0);
    case 4
      indx=(etable.TTL1==0 & etable.TTL2==0 & etable.TTL3==0 & etable.TTL4==1 & etable.TTL5==0 & etable.TTL6==0);
    case 5
      indx=(etable.TTL1==1 & etable.TTL2==0 & etable.TTL3==0 & etable.TTL4==0 & etable.TTL5==1 & etable.TTL6==1);
    case 6
      indx=(etable.TTL1==0 & etable.TTL2==1 & etable.TTL3==0 & etable.TTL4==0 & etable.TTL5==1 & etable.TTL6==1);
    case 7
      indx=(etable.TTL1==0 & etable.TTL2==0 & etable.TTL3==1 & etable.TTL4==0 & etable.TTL5==1 & etable.TTL6==1);
    case 8
      indx=(etable.TTL1==0 & etable.TTL2==0 & etable.TTL3==0 & etable.TTL4==1 & etable.TTL5==1 & etable.TTL6==1);
    case 9  % end of trial
      indx=(etable.TTL1==1 & etable.TTL2==1 & etable.TTL3==1 & etable.TTL4==1 & etable.TTL5==1 & etable.TTL6==1);
    case 10 % correct response
      indx=(etable.TTL1==1 & etable.TTL2==1 & etable.TTL3==0 & etable.TTL4==0 & etable.TTL5==0 & etable.TTL6==0);
    case 11 % incorrect response
      indx=(etable.TTL1==0 & etable.TTL2==0 & etable.TTL3==1 & etable.TTL4==1 & etable.TTL5==0 & etable.TTL6==0);
  end
  etable.Evt(indx)=t;
end

if fp==15
  etable{etable.Sample==587697,[3,9]}=[7,0];
end

% detect the square events and code for target and nontarget
% ==========================================================
atar=etable;
indx=find(ismember(atar.Evt,1:8));
if length(indx) ~= 320
  disp(length(indx))
  error('... Not 320 trials!')
end
indx(:,2)=[indx(2:end,1)-1; size(atar.Evt,1)];
for t = 1:size(indx,1)
  atar.Trial(indx(t,1):indx(t,2),1)=t;
  atar.Cond(indx(t,1):indx(t,2),1)=atar.Evt(indx(t,1));
end
% target square
atar.Tar=nan(size(atar,1),1);
atar.Tar(atar.TTL1==0 & atar.TTL2==0 & atar.TTL3==0 & atar.TTL4==0 & atar.TTL5==1 & atar.TTL6==0,1)=1;
% nontarget square
atar.Tar(atar.TTL1==0 & atar.TTL2==0 & atar.TTL3==0 & atar.TTL4==0 & atar.TTL5==0 & atar.TTL6==1,1)=0;
atar=atar(atar.Tar>=0,[1,12:14]);
if size(atar,1) ~= 1280
  disp(size(atar,1))
  error('... Not 1280 squares!')
end

% extract ssvep trials
% ====================
atrial=etable;
% delete squares events (not necessary for ssveps, they mess up the latencies between diodes)
atrial((atrial.TTL1==0 & atrial.TTL2==0 & atrial.TTL3==0 & atrial.TTL4==0 & atrial.TTL5==1 & atrial.TTL6==0),:)=[];
atrial((atrial.TTL1==0 & atrial.TTL2==0 & atrial.TTL3==0 & atrial.TTL4==0 & atrial.TTL5==0 & atrial.TTL6==1),:)=[];
% find trial onsets and offsets
ind=find(ismember(atrial.Evt,1:8));
ind=length(ind);
if ind ~= 320
  disp(length(ind))
  error('... Not 320 trials!')
end
assvep=array2table(nan(ind,5),'VariableNames',{'Sample','Trial','Cond','tf','Resp'});
assvep.tf(:,1)=find(atrial.Evt<9);
assvep.Trial=(1:320)';
assvep.Cond=atrial.Evt(assvep.tf);
assvep.tf(:,2)=find(atrial.Evt==9);
% how many frames between trial end and next trial?
assvep.tf(:,3)=([assvep.tf(2:end,1); size(atrial,1)+1] -assvep.tf(:,2))-1;
if max(assvep.tf(:,3))>1
  error('More than 1 row for response row?')
  % this should not occur because the response occurs immediately after the
  % end of a trials
end
% set no response to 0
assvep.Resp(assvep.tf(:,3)==0)=0;
% find response
ind=find(assvep.tf(:,3)>0);
assvep.Resp(ind)=atrial.Evt(assvep.tf(ind,2)+1);
% correct
assvep.Resp(assvep.Resp==10)=1;
% incorrect
assvep.Resp(assvep.Resp==11)=2;

% add the info about correct trials to target squares (atar)
ind=find(assvep.Resp==1);
atar.Correct=ismember(atar.Trial,ind);

% add Sample info (take light diode, not parallel port)
if any(atrial.Light1(assvep.tf(:,1)+1)~=1)
  error('Next event after trial onset is NOT a light diode!')
end
assvep.Sample=atrial.Sample(assvep.tf(:,1)+1);

% add trial light diodes (there should be 38 cycles of 12 frames after trial onset)
assvep.Latency(:,1:39)=nan(size(assvep,1),39);
for t=1:size(assvep,1)
  % latency of all diode responses
  % 1st interval is that between parallel port and 1st diode
  % after that, list all diodes before end of trial code
  ind=assvep.tf(t,1):assvep.tf(t,2)-1;
  if any(atrial.Light1(ind(2:end))~=1)
    error('Not all rows in a trial are light diodes!')
  end
  tmp2=(atrial.Sample(ind)-atrial.Sample(ind(1)));
  % in ms
  %tmp2=tmp2/512*1000;
  assvep.Latency(t,1:length(tmp2))=tmp2;
end
assvep.tf=[];

% test=assvep{:,:};
% test=diff(assvep.Latency,1,2);
% test=test/512*1000;
% test(isnan(test))=[];
% range(test)
% hist(test)
% quantile(test,[0.025,0.975])
% quantile(test,[0.01,0.99])
% hist(test(:))

% create trl
% ==========
% for ssvep
ind = 320;
if size(assvep,1) ~= ind
  error(['Error: ',num2str(ind),' unequals 320 trials total!'])
end
assvep2=assvep{:,:};

% cfg.trialdef.prestim = 0; % in s <--- positive value!
% cfg.trialdef.poststim = 4.5; % in s
% begsample = samplenow - round(cfg.trialdef.prestim  * hdr.Fs);
% endsample = samplenow + round(cfg.trialdef.poststim * hdr.Fs);
% offset    = -round(cfg.trialdef.prestim  * hdr.Fs);
preSSVEP = 0.5;       % positive number for interval before onset
lengthSSVEP = 4.5;
begsample = assvep.Sample - round(preSSVEP * hdr.Fs);      % 0 s before
endsample = assvep.Sample + round(lengthSSVEP * hdr.Fs);    % 4.5 s after
offset = repmat(-round(preSSVEP * hdr.Fs),ind,1);
ssvep = ones(320,1);
trl = [begsample endsample offset ssvep assvep2(:,2:end)];

% append ERPs to squares
preERP = 1;       % positive number for interval before onset
lengthERP = 1;
ind = size(atar);
begsample = atar.Sample - round(preERP * hdr.Fs);
endsample = atar.Sample + round(lengthERP * hdr.Fs);
offset = repmat(-round(preERP * hdr.Fs),ind(1),1);
ssvep = zeros(ind(1),1);
trl2 = [begsample endsample offset ssvep atar{:,2:end}];
trl2 = [trl2, nan(ind(1),size(trl,2)-size(trl2,2))];
trl = vertcat(trl,trl2);

