The following is an example analysis pipeline that was used for the FieldTrip workshop at CUNY, New York by Stephen and Saskia in 2011.

% load biosemi160 into q (which has sensor locations in polar coordinates
% add fiducials then save

% cfg=[];
% cfg.layout='biosemi160.lay';
% q=ft_prepare_layout(cfg);

zeroline = 92/40*50;

q{161,1} = 'nasion';
q{161,2} = zeroline; %ten percent under FPz
q{161,3} = q{104,3};

q{162,1} = 'inion';
q{162,2} = zeroline; %ten percent under Oz
q{162,3} = q{23,3};

q{163,1} = 'left';
q{163,2} = -zeroline; %ten percent under T7
q{163,3} = q{137,3};

q{164,1} = 'right';
q{164,2} = zeroline; %ten percent under T8
q{164,3} = q{64,3};

save('LayoutNY_fiducials',q);



% construct 3D electrode positions

load LayoutNY_fiducials.mat

ph = cell2mat(q(:,2));
th = cell2mat(q(:,3));

x = sin(ph*pi/180) .* cos(th*pi/180);
y = sin(ph*pi/180) .* sin(th*pi/180);
z = cos(ph*pi/180);

plot3(x, y, z, '.');
elec.label = q(:,1);
elec.pnt = [x y z];

% scale to get into mm
elec.pnt = 100*elec.pnt;

vol = ft_read_vol('headmodel/standard_vol.mat');

%realign electrodes to headmodel
cfg = [];
cfg.method = 'interactive';
cfg.elec = elec;
cfg.headshape = vol.bnd(1); %1 = skin
elec = ft_electroderealign(cfg);

%save new electrodes
save elec160.mat elec

clear all

cond = {'Pitch20', 'Pitch40', 'Timbre490', 'Timbre510'};

for c = 1:4 % loop over the 4 conditions
  for block = 1:4 % loop over the 4 blocks within each condition
    
    % create the trial definition
    cfg=[];
    cfg.filename              = ['data/RON_', cond{c}, '_Block', num2str(block), '_Sub03.bdf'];
    cfg.trialfun              = 'trialfunNY';
    cfg.trialdef.eventtype    = 'STATUS';
    cfg.trialdef.eventvalue   = [121 122 103 104 111 112]; % stimulus triggers (standard before (2), deviant (2), standard after(2))
    cfg.trialdef.eventcorrect = [1   2   1   2   1   2  ]; % correct response triggers
    cfg.trialdef.prestim      = 0.2; % latency in seconds
    cfg.trialdef.poststim     = 1;   % latency in seconds
    
    cfg = ft_definetrial(cfg);
    
    % remember condition and block
    trl = cfg.trl;
    trl(:,6) = c;
    trl(:,8) = block;
    
    % define the relevant trials
    for i=1:length(trl)
      if ismember(trl(i,4), [121 122])
        trl(i,7) = 1; % standard before
      elseif ismember(trl(i,4), [103 104])
        trl(i,7) = 2; % deviant
      elseif ismember(trl(i,4), [111 112])
        trl(i,7) = 3; % standard after
      end
    end
    
    % read and preprocess the data using the trial definition
    cfg=[];
    cfg.dataset       = ['data/RON_', cond{c}, '_Block', num2str(block), '_Sub03.bdf'];
    cfg.trl           = trl;
    cfg.reref         = 'yes';
    cfg.refchannel    = 'EXG5';
    
    data = ft_preprocessing(cfg);
    
    % save the data to disk
    save(['analysis/prepro/data_', cond{c}, '_block', num2str(block)], 'data')
  end
end

clear all
cond = {'Pitch20' 'Pitch40', 'Timbre490', 'Timbre510'};

% concatenate all trials/blocks
ipart=1;
for icond = 1:4
  for iblock = 1:4
    disp(['Loading analysis/prepro/data_', cond{icond}, '_block', num2str(iblock)]);
    load(['analysis/prepro/data_', cond{icond}, '_block', num2str(iblock)], 'data');
    datapart(ipart) = data;
    ipart=ipart+1;
  end
end

cfg=[];
data = ft_appenddata(cfg, datapart(1), datapart(2),datapart(3),datapart(4),datapart(5),datapart(6),datapart(7),datapart(8),datapart(9),datapart(10),datapart(11),datapart(12),datapart(13),datapart(14),datapart(15),datapart(16));
clear datapart*

% detect eog artifacts using ICA
cfg=[];
cfg.method  = 'runica';
cfg.channel = 1:160; % EEG channels only
datacomp = ft_componentanalysis(cfg, data);
save('analysis/ica/datacomp', 'datacomp')

% plot the components to detect the artifacts
figure
k=1; f=1;
for icomp=1:length(datacomp.topo)
  if k>20
    k=1;
    figure
  end
  cfg=[];
  cfg.layout = 'biosemi160.lay';
  cfg.xlim   = [icomp icomp];
  
  subplot(4,5,k);
  ft_topoplotER(cfg, datacomp);
  title(icomp);
  
  k = k+1;
end


% remove components that reflect eog artifacts
cfg=[];
cfg.component = [12 49]; % note the exact numbers will vary per run
data = ft_rejectcomponent(cfg, datacomp);
save('analysis/data_clean', 'data')



%% manual artifact rejection

% remove the mean
cfg=[];
cfg.demean = 'yes';
data = ft_preprocessing(cfg, data);

% shuffle the trials for non-biased artifact rejection
shuffle = randperm(length(data.trial));
datashuff = data;
for i=1:length(data.trial)
  datashuff.trial{i} = data.trial{shuffle(i)};
  datashuff.time{i}  = data.time{shuffle(i)};
  datashuff.trialinfo(i,:) = data.trialinfo(shuffle(i),:);
end

% browse for artifact rejection
cfg = [];
cfg.channel = 'EEG';
cfg.continuous = 'no';
cfg = ft_databrowser(cfg, data);
data = ft_rejectartifact(cfg,data);

% visual artifact rejection in summary mode
cfg = [];
cfg.method = 'summary';
data = ft_rejectvisual(cfg, data);

% the following trials were removed: 43, 87, 174, 261, 307, 333, 343, 398, 557, 583, 593, 606, 666, 670, 674, 741, 742, 743, 744, 745, 899
% the following channels were removed: A7, A13, A14, A25, A26, A27, A28, B8, B9, B12, B19, B32, D1, D21, D23, D25, D26, D31, D32, E1, E17, E21, E32


save('analysis/data_clean', 'data')

clear all

cond = {'Pitch20' 'Pitch40', 'Timbre490', 'Timbre510'};

load('analysis/data_clean', 'data')

for icond = 1:4 % loop over the 4 conditions    
  for istim = 1:3 % loop over standard-before, deviant, standard-after
      
    % compute timelocked averages for each condition
    cfg=[];
    cfg.lpfilter   = 'yes';
    cfg.lpfreq     = 40;
    cfg.trials   = find(data.trialinfo(:,3) == icond & data.trialinfo(:,4) == istim);
    timelock{icond,istim} = ft_timelockanalysis(cfg, data);
      
    % baseline correct
    cfg=[];
    cfg.baseline = [-0.2 0];
    timelock{icond,istim} = ft_timelockbaseline(cfg, timelock{icond,istim});
  end  
end

save('analysis/timelock/timelock_all','timelock');

% plot the ERPs over all sensors
figure
for i = 1:4
  subplot(2,2,i);
  ft_singleplotER([],timelock{i,1},timelock{i,2},timelock{i,3});
  title(cond{i});
  if i==1
    legend('standard-before', 'deviant', 'standard-after', 'Location', 'NorthWest')
  end
end
print(gcf, '-dpng', 'figures/fig1_ERP')

% plot ERP in interactive mode, only for standard-before
cfg = [];
cfg.layout      = 'biosemi160lay.mat';
cfg.interactive = 'yes';
figure; ft_multiplotER(cfg,timelock{1,1});

% compute the contrasts
for icond = 1:4
  % standard before vs deviant
  stb_vs_dev{icond} = timelock{1,1};
  stb_vs_dev{icond}.avg = timelock{icond,2}.avg - timelock{icond,1}.avg;
  
  % standard before vs standard after
  stb_vs_sta{icond} = timelock{1,1};
  stb_vs_sta{icond}.avg = timelock{icond,3}.avg - timelock{icond,1 }.avg;
  
end

% plot the contrasts
figure
for icond = 1:4
  cfg=[];
  cfg.xlim        = [0.5 7];
  cfg.zlim        = 'maxabs';
  cfg.interactive = 'yes';
  cfg.layout      = 'biosemi160lay.mat';
  cfg.colorbar    = 'yes';
  
  subplot(2,2,icond);
  ft_topoplotER(cfg, stb_vs_dev{icond});
end
print(gcf, '-dpng', 'figures/fig2_ERP')

figure
for icond = 1:4
  cfg=[];
  cfg.xlim        = [0.5 7];
  cfg.zlim        = 'maxabs';
  cfg.interactive = 'yes';
  cfg.layout      = 'biosemi160lay.mat';
  cfg.colorbar    = 'yes';
  
  subplot(2,2,icond);
  ft_topoplotER(cfg, stb_vs_sta{icond});
end
print(gcf, '-dpng', 'figures/fig3_ERP')



%% now collapse

% compute averages for pitch and timbre
for istim = 1:3
  cfg=[];
  cfg.lpfilter   = 'yes';
  cfg.lpfreq     = 40;
  cfg.keeptrials = 'yes';
  cfg.trials     = find((data.trialinfo(:,3) == 1 | data.trialinfo(:,3) == 2 ) & data.trialinfo(:,4) == istim);
  timelock_pitch{istim}  = ft_timelockanalysis(cfg, data);
  
  cfg.trials     = find((data.trialinfo(:,3) == 3 | data.trialinfo(:,3) == 4 ) & data.trialinfo(:,4) == istim);
  timelock_timbre{istim} = ft_timelockanalysis(cfg, data);
  
  cfg=[];
  cfg.baseline   = [-0.2 0];
  timelock_pitch{istim}  = ft_timelockbaseline(cfg, timelock_pitch{istim});
  timelock_timbre{istim} = ft_timelockbaseline(cfg, timelock_timbre{istim});
end
save('analysis/timelock/timelock_avg','timelock_*');


% plot contrasts
figure
subplot(1,2,1);
ft_singleplotER([],timelock_pitch{1},timelock_pitch{2},timelock_pitch{3});
title('Pitch');

subplot(1,2,2);
ft_singleplotER([],timelock_timbre{1},timelock_timbre{2},timelock_timbre{3});
title('Timbre');

print(gcf, '-dpng', 'figures/fig4_ERP')

cfg=[];
cfg.layout = 'biosemi160lay.mat'; %in meters
%cfg.layout        = 'elec160.mat'; %in mm
cfg.neighbourdist = .1; 
cfg.neighbours    = ft_neighbourselection(cfg,timelock_pitch{1});

cfg.latency       = [0 1];
cfg.parameter     = 'trial';
cfg.method        = 'montecarlo';
cfg.design        = [1:size(timelock_pitch{2}.trial,1) 1:size(timelock_pitch{1}.trial,1);
  ones(1,size(timelock_pitch{2}.trial,1)), ones(1,size(timelock_pitch{1}.trial,1))*2];

cfg.numrandomization = 1000;
cfg.correctm      = 'cluster';
cfg.correcttail   = 'prob';
cfg.ivar          = 2;
cfg.uvar          = 1;

cfg.statistic     = 'indepsamplesT';

stat = ft_timelockstatistics(cfg, timelock_pitch{2}, timelock_pitch{1});

save('analysis/timelock/stat','stat');



%% plot

% cfg=[];
% cfg.layout = 'biosemi160lay.mat';
% ft_clusterplot(cfg, stat)


% find relevant clusters
ipos = find([stat.posclusters.prob]<=0.05);
ineg = find([stat.negclusters.prob]<=0.05);

% loop over all sig positive clusters
for i=ipos
  
  cfg=[];
  cfg.highlight = 'on';
  cfg.zparam    = 'stat';
  cfg.layout    = 'biosemi160lay.mat';
  cfg.style     = 'straight';
  cfg.gridscale = 500;
  
  % find the significant time range for this cluster
  tmp=[];
  for t = 1:length(stat.time)
    if ~isempty(find(any(stat.posclusterslabelmat(:,t)==ipos)))
      tmp = [tmp t];
    end
  end
  cfg.xlim      = [stat.time(tmp(1)) stat.time(tmp(end))];
  
  % find the channels belonging to this cluster
  cfg.highlightchannel = [];
  for c = 1:length(stat.label)
    if ~isempty(find(any(stat.posclusterslabelmat(c,:)==ipos)))
      cfg.highlightchannel = [cfg.highlightchannel c];
    end
  end
  
  figure
  ft_topoplotER(cfg, stat);
  title('positive cluster')
  print(gcf, '-dpng', ['figures/fig5_STAT_pos', num2str(i)])
end

% loop over all sig negative clusters
for i=ineg
  
  cfg=[];
  cfg.highlight = 'on';
  cfg.zparam    = 'stat';
  cfg.layout    = 'biosemi160lay.mat';
  cfg.style     = 'straight';
  cfg.gridscale = 500;
  
  % find the significant time range for this cluster
  tmp=[];
  for t = 1:length(stat.time)
    if ~isempty(find(any(stat.negclusterslabelmat(:,t)==ineg)))
      tmp = [tmp t];
    end
  end
  cfg.xlim      = [stat.time(tmp(1)) stat.time(tmp(end))];
  
  % find the channels belonging to this cluster
  cfg.highlightchannel = [];
  for c = 1:length(stat.label)
    if ~isempty(find(any(stat.negclusterslabelmat(c,:)==ineg)))
      cfg.highlightchannel = [cfg.highlightchannel c];
    end
  end
  
  figure
  ft_topoplotER(cfg, stat);
  title('negative cluster')
  print(gcf, '-dpng', ['figures/fig6_STAT_neg', num2str(i)])
end



load analysis/timelock/timelock_all.mat
load headmodel/standard_mri.mat

elec = ft_read_sens('elec160.mat');
vol  = ft_read_vol('headmodel/standard_vol.mat');

  vol.r = [86 88 92 100];
  vol.o = [0 0 40];
  figure, ft_plot_vol(vol)
ft_plot_vol([],vol);

% prepare data
cfg = [];
cfg.demean = 'yes';
cfg.baselinewindow = [-inf 0];
cfg.reref = 'yes';
cfg.refchannel = 'all'; % thereby averaging systematic leadfield error over sensors
for icond = 1 : 4
    for istim = 1 : 3
        timelock{icond,istim} = ft_preprocessing(cfg, timelock{icond,istim});
    end
end

% do some averages
before = timelock{1,1};
before.trial{1} = (timelock{1,1}.trial{1} + timelock{2,1}.trial{1} + timelock{3,1}.trial{1} + timelock{4,1}.trial{1} ) ./ 4;
deviant = timelock{1,1};
deviant.trial{1} = (timelock{1,2}.trial{1} + timelock{2,2}.trial{1} + timelock{3,2}.trial{1} + timelock{4,2}.trial{1} ) ./ 4;
after = timelock{1,1};
after.trial{1} = (timelock{1,3}.trial{1} + timelock{2,3}.trial{1} + timelock{3,3}.trial{1} + timelock{4,3}.trial{1} ) ./ 4;

% plot ERP's
cfg = [];
cfg.layout = 'biosemi160lay.mat';
cfg.interactive = 'yes';
figure; ft_multiplotER(cfg, before, deviant, after);
legend;

% dipole fitting
cfg = [];
cfg.vol = vol;
cfg.elec = elec;
cfg.model = 'regional';
cfg.numdipoles = 2;
cfg.latency = [0.080 0.1]; %n100
cfg.gridsearch = 'yes';
cfg.resolution = 20; % mm
cfg.nonlinear = 'yes';
cfg.symmetry = 'x';
source = ft_dipolefitting(cfg, before);

% topoplot ERP
cfg = [];
cfg.layout = 'biosemi160lay.mat';
cfg.zparam = 'Vdata';
figure; ft_topoplotER(cfg, source);

% topoplot dipole model
cfg.zparam = 'Vmodel';
figure; ft_topoplotER(cfg, source);

% topoplot difference data with model
source.Vdifference = source.Vdata ./ source.Vmodel;
cfg.zparam = 'Vdifference';
cfg.zlim = [0 2];
figure; ft_topoplotER(cfg, source);

% plot first dipole position on mni
cfg = [];
cfg.method = 'ortho';
cfg.location = source.dip.pos(1,:);
figure; ft_sourceplot(cfg, mri);

% plot second dipole position on mni
cfg.location = source.dip.pos(2,:);
figure; ft_sourceplot(cfg, mri);

% use dipole positions to extract timecourse of whole trial
cfg = source.cfg;
cfg.dip.pos = source.dip.pos;
cfg.gridsearch = 'no';
cfg.nonlinear = 'no';
cfg.latency = [-inf inf];
cfg.vol = vol;
cfg.elec = elec;
source2 = ft_dipolefitting(cfg, before);

% plot timecourses of both dipoles
figure; hold; 
plot(source2.time, source2.dip.mom(1:3,:), '-')
plot(source2.time, source2.dip.mom(4:6,:), ':')

% use Singular Value Decomposition to extract component from all
% three orientations of first dipole. u = identity matrix, s = scaling, v =
% signal strength
[u1, s1, v1] = svd(source2.dip.mom(1:3,:));
[u2, s2, v2] = svd(source2.dip.mom(4:6,:));
figure; hold;
plot(source2.time, v1(:,1)); 
plot(source2.time, v2(:,1),':'); 

% use pythagoras to make one dimensional dipole from three orientations
figure; hold;
plot(source2.time, sqrt(sum(source2.dip.mom(1:3,:).^2)));
plot(source2.time, sqrt(sum(source2.dip.mom(4:6,:).^2)),':');