% POP_COMPERP - Compute the grand average ERP waveforms of multiple datasets
%                 currently loaded into EEGLAB, with optional ERP difference-wave 
%                 plotting and t-tests. Creates a plotting figure.
% Usage:
%       >> pop_comperp( ALLEEG, flag );  % pop-up window, interactive mode
%       >> [erp1 erp2 erpsub time sig] = pop_comperp( ALLEEG, flag, ...
%                                         datadd, datsub, 'key', 'val', ...);
% Inputs:
%   ALLEEG  - Array of loaded EEGLAB EEG structure datasets
%   flag    - [0|1] 0 -> Use ICA components; 1 -> use data channels {default: 1}
%   datadd  - [integer array] List of ALLEEG dataset indices to average to make 
%             an ERP grand average and optionally to compare with 'datsub' datasets.
%
% Optional inputs:
%   datsub  - [integer array] List of ALLEEG dataset indices to average and then 
%             subtract from the 'datadd' result to make an ERP grand mean difference. 
%             Together, 'datadd' and 'datsub' may be used to plot and compare grand mean 
%             responses across subjects or conditions. Both arrays must contain the same 
%             number of dataset indices and entries must be matched pairwise (Ex:
%             'datadd' indexes condition A datasets from subjects 1:n, and 'datsub',
%             condition B datasets from the same subjects 1:n). {default: []}
%   'alpha'    - [0 < float < 1] Apply two-tailed t-tests for p < alpha. If 'datsub' is
%                not empty, perform t-tests at each latency. If 'datasub' is empty, 
%                perform two-tailed t-tests against a 0 mean dataset with same variance. 
%                Significant time regions are highlighted in the plotted data. 
%   'chans'    - [integer array] Vector of chans. or comps. to use {default: all}
%   'geom'     - ['scalp'|'array'] Plot erps in a scalp array (PLOTTOPO)
%                or as a rectangular array (PLOTDATA). Note: Only channels
%                (see 'chans' above) can be plotted in a 'scalp' array.
%   'tlim'     - [min max] Time window (ms) to plot data {default: whole time range}
%   'title'    - [string] Plot title {default: none}
%   'ylim'     - [min max] y-axis limits {default: auto from data limits}
%   'mode'     - ['ave'|'rms'] Plotting mode. Plot  either grand average or RMS 
%                (root mean square) time course(s) {default: 'ave' -> grand average}.
%   'std'      - ['on'|'off'|'none'] 'on' -> plot std. devs.; 'none' -> do not 
%                interact with other options {default:'none'}
%
% Visualisation options:
%   'addavg'   - ['on'|'off'] Plot grand average (or RMS) of 'datadd' datasets
%                {default: 'on' if 'datsub' empty, otherwise 'off'}
%   'subavg'   - ['on'|'off'] Plot grand average (or RMS) of 'datsub' datasets 
%                {default:'off'}
%   'diffavg'  - ['on'|'off'] Plot grand average (or RMS) difference 
%   'addall'   - ['on'|'off'] Plot the ERPs for all 'dataadd' datasets only {default:'off'}
%   'suball'   - ['on'|'off'] Plot the ERPs for all 'datasub datasets only {default:'off'}
%   'diffall'  - ['on'|'off'] Plot all the 'datadd'-'datsub' ERP differences {default:'off'}
%   'addstd'   - ['on'|'off'] Plot std. dev. for 'datadd' datasets only 
%                {default: 'on' if 'datsub' empty, otherwise 'off'}
%   'substd'   - ['on'|'off'] Plot std. dev. of 'datsub' datasets only {default:'off'}
%   'diffstd'  - ['on'|'off'] Plot std. dev. of 'datadd'-'datsub' differences {default:'on'}
%   'diffonly' - ['on'|'off'|'none'] 'on' -> plot difference only; 'none' -> do not affect 
%                other options {default:'none'}
%   'allerps'  - ['on'|'off'|'none'] 'on' -> show ERPs for all conditions; 
%                'none' -> do not affect other options {default:'none'}
%   'tplotopt' - [cell array] Pass 'key', val' plotting options to PLOTTOPO
%
% Output:
%   erp1   - Grand average (or rms) of the 'datadd' datasets
%   erp2   - Grand average (or rms) of the 'datsub' datasets
%   erpsub - Grand average (or rms) 'datadd' minus 'datsub' difference
%   times  - Vector of epoch time indices
%   sig    - T-test significance values (chans,times). 
%
% Author: Arnaud Delorme, CNL / Salk Institute, 15 March 2003
%
% Note: t-test functions were adapted for matrix preprocessing from C functions 
%       by Press et al. See the description in the PTTEST code below 
%       for more information.
%
% See also: EEGLAB, PLOTTOPO

% Copyright (C) 15 March 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [erp1, erp2, erpsub, times, pvalues] = pop_comperp( ALLEEG, flag, datadd, datsub, varargin);

erp1 = '';
if nargin < 1
   help pop_comperp;
   return;
end
if isempty(ALLEEG)
    error('pop_comperp: cannot process empty sets of data');
end
if nargin < 2
   flag = 1;
end
allcolors = { 'b' 'r' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' ...
              'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm'};
erp1 = '';
if nargin < 3
    gtmp = [1.1 0.8 .21 .21 .21 0.1]; gtmp2 = [1.48 1.03 1];
    uigeom = { [2.6 0.95] gtmp gtmp gtmp [1] gtmp2 gtmp2 [1.48 0.25 1.75] gtmp2 gtmp2 };
    commulcomp= ['if get(gcbo, ''value''),' ...
                 '    set(findobj(gcbf, ''tag'', ''multcomp''), ''enable'', ''on'');' ...
                 'else,' ...
                 '    set(findobj(gcbf, ''tag'', ''multcomp''), ''enable'', ''off'');' ...
                 'end;'];
	uilist = { { } ...
               { 'style' 'text' 'string' 'avg.        std.      all ERPs' } ...
               { 'style' 'text' 'string' 'Datasets to average (ex: 1 3 4):' } ...
               { 'style' 'edit' 'string' '' } ...
               { 'style' 'checkbox' 'string' '' 'value' 1 } ...
               { 'style' 'checkbox' 'string' '' } ...
               { 'style' 'checkbox' 'string' '' } { } ...
	           { 'style' 'text' 'string' 'Datasets to average and subtract (ex: 5 6 7):' } ...
               { 'style' 'edit' 'string' '' } ...
               { 'style' 'checkbox' 'string' '' 'value' 1 } ...
               { 'style' 'checkbox' 'string' '' } ...
               { 'style' 'checkbox' 'string' '' } { } ...
	           { 'style' 'text' 'string' 'Plot difference' } { } ...
               { 'style' 'checkbox' 'string' '' 'value' 1 } ...
               { 'style' 'checkbox' 'string' '' } ...
               { 'style' 'checkbox' 'string' '' } { } ...
               { } ...
	           { 'style' 'text' 'string' fastif(flag, 'Channels subset ([]=all):', ...
                                                  'Components subset ([]=all):') } ...
               { 'style' 'edit' 'string' '' } { } ...
	           { 'style' 'text' 'string' 'Highlight significant regions (.01 -> p=.01)' } ...
               { 'style' 'edit' 'string' '' } { } ...
	           { 'style' 'text' 'string' 'Use RMS instead of average (check):' } { 'style' 'checkbox' 'string' '' } { } ...
	           { 'style' 'text' 'string' 'Low pass (Hz) (for display only)' } ...
               { 'style' 'edit' 'string' '' } { } ...
               { 'style' 'text' 'string' 'Plottopo options (''key'', ''val''):' } ...
               { 'style' 'edit' 'string' '''ydir'', -1' } ...
               { 'style' 'pushbutton' 'string' 'Help' 'callback', 'pophelp(''plottopo'')' } ...
               };
    
    % remove geometry textbox for ICA components
    result = inputgui( uigeom, uilist, 'pophelp(''pop_comperp'')', 'ERP grand average/RMS - pop_comperp()');
    if length(result) == 0, return; end

    %decode parameters list
    options = {};
    datadd = eval( [ '[' result{1} ']' ]);
    if result{2},  options = { options{:} 'addavg'  'on' }; else, options = { options{:} 'addavg'  'off' }; end
    if result{3},  options = { options{:} 'addstd'  'on' }; else, options = { options{:} 'addstd'  'off' }; end
    if result{4},  options = { options{:} 'addall'  'on' }; end
    datsub = eval( [ '[' result{5} ']' ]);
    if result{6},  options = { options{:} 'subavg'  'on' }; end
    if result{7},  options = { options{:} 'substd'  'on' }; end
    if result{8},  options = { options{:} 'suball'  'on' }; end
    if result{9},  options = { options{:} 'diffavg' 'on' }; else, options = { options{:} 'diffavg' 'off' }; end
    if result{10}, options = { options{:} 'diffstd' 'on' }; else, options = { options{:} 'diffstd' 'off' }; end
    if result{11}, options = { options{:} 'diffall' 'on' }; end
    
    if result{12},           options = { options{:} 'chans' eval( [ '[' result{12} ']' ]) }; end
    if ~isempty(result{13}), options = { options{:} 'alpha' str2num(result{13}) }; end
    if result{14},           options = { options{:} 'mode' 'rms' }; end
    if ~isempty(result{15}), options = { options{:} 'lowpass' str2num(result{15}) }; end
    if ~isempty(result{16}), options = { options{:} 'tplotopt' eval([ '{ ' result{16} ' }' ]) }; end; 
else 
    options = varargin;
end

if nargin == 3
  datsub = []; % default 
end

% decode inputs
% -------------
if isempty(datadd), error('First edit box (datasets to add) can not be empty'); end
g = finputcheck( options, ... 
                 { 'chans'    'integer'  []               [1:ALLEEG(datadd(1)).nbchan];
                   'title'    'string'   []               '';
                   'alpha'    'float'    []               [];
                   'geom'     'string'  {'scalp';'array'} fastif(flag, 'scalp', 'array');
                   'addstd'   'string'  {'on';'off'}      fastif(isempty(datsub), 'on', 'off');
                   'substd'   'string'  {'on';'off'}     'off';
                   'diffstd'  'string'  {'on';'off'}     'on';
                   'addavg'   'string'  {'on';'off'}     fastif(isempty(datsub), 'on', 'off');
                   'subavg'   'string'  {'on';'off'}     'off';
                   'diffavg'  'string'  {'on';'off'}     'on';
                   'addall'   'string'  {'on';'off'}     'off';
                   'suball'   'string'  {'on';'off'}     'off';
                   'diffall'  'string'  {'on';'off'}     'off';
                   'std'      'string'  {'on';'off';'none'}     'none';
                   'diffonly' 'string'  {'on';'off';'none'}     'none';
                   'allerps'  'string'  {'on';'off';'none'}     'none';
                   'lowpass'  'float'    [0 Inf]         [];
                   'tlim'     'float'    []              [];
                   'ylim'     'float'    []              [];
                   'tplotopt' 'cell'     []              {};
                   'mode'     'string'  {'ave';'rms'}    'ave';
                   'multcmp'  'integer'  [0 Inf]         [] });
if ischar(g), error(g); end
if length(datadd) == 1
    disp('Cannot perform statistics using only 1 dataset');
    g.alpha = [];
end

h = figure; axcopy

try, icadefs; set(h, 'color', BACKCOLOR); axis off; catch, end

% backward compatibility of param
% -------------------------------
if ~strcmpi(g.diffonly, 'none')
    if strcmpi(g.diffonly, 'off'), g.addavg = 'on'; g.subavg = 'on'; end
end
if ~strcmpi(g.allerps, 'none')
    if isempty(datsub)
         g.addall  = g.allerps;
    else g.diffall = g.allerps;
    end
end
if ~strcmpi(g.std, 'none')
    if isempty(datsub)
         g.addstd  = g.std;
    else g.diffstd = g.std;
    end
end

% check consistency
% -----------------
if length(datsub) > 0 && length(datadd) ~= length(datsub)
    error('The number of component to subtract must be the same as the number of components to add');
end

% if only 2 dataset entered, toggle average to single trial
% ---------------------------------------------------------
if length(datadd) == 1 &strcmpi(g.addavg, 'on')
    g.addavg  = 'off';
    g.addall = 'on';
end
if length(datsub) == 1 &strcmpi(g.subavg, 'on')
    g.subavg  = 'off';
    g.suball = 'on';
end
if length(datsub) == 1 && length(datadd) == 1 &strcmpi(g.diffavg, 'on')
    g.diffavg  = 'off';
    g.diffall = 'on';
end;      

regions = {};
pnts   = ALLEEG(datadd(1)).pnts;
srate  = ALLEEG(datadd(1)).srate;
xmin   = ALLEEG(datadd(1)).xmin;
xmax   = ALLEEG(datadd(1)).xmax;
nbchan = ALLEEG(datadd(1)).nbchan;
chanlocs = ALLEEG(datadd(1)).chanlocs;
unionIndices = union_bc(datadd, datsub);
for index = unionIndices(:)'
    if ALLEEG(index).pnts ~= pnts,     error(['Dataset '  int2str(index) ' does not have the same number of points as others']); end
    if ALLEEG(index).xmin ~= xmin,     error(['Dataset '  int2str(index) ' does not have the same xmin as others']); end
    if ALLEEG(index).xmax ~= xmax,     error(['Dataset '  int2str(index) ' does not have the same xmax as others']); end
    if ALLEEG(index).nbchan ~= nbchan, error(['Dataset '  int2str(index) ' does not have the same number of channels as others']); end
end
    
if ~isempty(g.alpha) && length(datadd) == 1
    error([ 'T-tests require more than one ''' datadd ''' dataset' ]);
end

% compute ERPs for add
% --------------------
for index = 1:length(datadd)
    TMPEEG = eeg_checkset(ALLEEG(datadd(index)),'loaddata');
    if flag == 1, erp1ind(:,:,index)  = mean(TMPEEG.data,3);
    else          erp1ind(:,:,index)  = mean(eeg_getdatact(TMPEEG, 'component', [1:size(TMPEEG.icaweights,1)]),3);
    end
    addnames{index} = [ '#' int2str(datadd(index)) ' ' TMPEEG.setname ' (n=' int2str(TMPEEG.trials) ')' ];
    clear TMPEEG;
end

% optional: subtract
% ------------------
colors = {}; % color aspect for curves
allcolors = { 'b' 'r' 'g' 'c' 'm' 'y' [0 0.5 0] [0.5 0 0] [0 0 0.5] [0.5 0.5 0] [0 0.5 0.5] [0.5 0 0.5] [0.5 0.5 0.5] };
allcolors = { allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} };
allcolors = { allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} allcolors{:} };
if length(datsub) > 0 % dataset to subtract

    % compute ERPs for sub
    % --------------------
    for index = 1:length(datsub)
        TMPEEG = eeg_checkset(ALLEEG(datsub(index)),'loaddata');
        if flag == 1, erp2ind(:,:,index)  = mean(TMPEEG.data,3);
        else          erp2ind(:,:,index)  = mean(eeg_getdatact(TMPEEG, 'component', [1:size(TMPEEG.icaweights,1)]),3);
        end
        subnames{index} = [ '#' int2str(datsub(index)) ' ' TMPEEG.setname '(n=' int2str(TMPEEG.trials) ')' ];
        clear TMPEEG
    end
    
    l1 = size(erp1ind,3);
    l2 = size(erp2ind,3);
    allcolors1 = allcolors(3:l1+2);
    allcolors2 = allcolors(l1+3:l1+l2+3);
    allcolors3 = allcolors(l1+l2+3:end);
    [erps1, erpstd1, colors1, colstd1, legend1] = preparedata( erp1ind        , g.addavg , g.addstd , g.addall , g.mode, 'Add ' , addnames, 'b', allcolors1 );
    [erps2, erpstd2, colors2, colstd2, legend2] = preparedata( erp2ind        , g.subavg , g.substd , g.suball , g.mode, 'Sub ' , subnames, 'r', allcolors2 );
    [erps3, erpstd3, colors3, colstd3, legend3] = preparedata( erp1ind-erp2ind, g.diffavg, g.diffstd, g.diffall, g.mode, 'Diff ', ...
                                                      { addnames subnames }, 'k', allcolors3 );
    
    % handle special case of std
    % --------------------------
    erptoplot  = [ erps1 erps2 erps3 erpstd1 erpstd2 erpstd3 ];
    colors     = { colors1{:} colors2{:} colors3{:} colstd1{:} colstd2{:} colstd3{:}};
    legend     = { legend1{:} legend2{:} legend3{:} };
    
    % highlight significant regions
    % -----------------------------
    if ~isempty(g.alpha)
        pvalues = pttest(erp1ind(g.chans,:,:), erp2ind(g.chans,:,:), 3);
        regions = p2regions(pvalues, g.alpha, [xmin xmax]*1000);
    else 
        pvalues= [];
    end
    
else
    [erptoplot, erpstd, colors, colstd, legend] = preparedata( erp1ind, g.addavg, g.addstd, g.addall, g.mode, '', addnames, 'k', allcolors);
    erptoplot = [ erptoplot erpstd ];
    colors    = { colors{:} colstd{:} };
    
    % highlight significant regions
    % -----------------------------
    if ~isempty(g.alpha)
        pvalues = ttest(erp1ind, 0, 3);
        regions = p2regions(pvalues, g.alpha, [xmin xmax]*1000);
    else 
        pvalues= [];
    end
end
    
% lowpass data
% ------------
if ~isempty(g.lowpass)
    if exist('filtfilt') == 2
        erptoplot = eegfilt(erptoplot, srate, 0, g.lowpass);
    else
        erptoplot = eegfiltfft(erptoplot, srate, 0, g.lowpass);
    end
end
if strcmpi(g.geom, 'array') || flag == 0, chanlocs = []; end
if ~isfield(chanlocs, 'theta'), chanlocs = []; end

% select time range
% -----------------
if ~isempty(g.tlim)
    pointrange = round(eeg_lat2point(g.tlim/1000, [1 1], srate, [xmin xmax]));
    g.tlim     = eeg_point2lat(pointrange, [1 1], srate, [xmin xmax]);
    erptoplot  = reshape(erptoplot, size(erptoplot,1), pnts, size(erptoplot,2)/pnts);
    erptoplot  = erptoplot(:,[pointrange(1):pointrange(2)],:);
    pnts       = size(erptoplot,2);
    erptoplot  = reshape(erptoplot, size(erptoplot,1), pnts*size(erptoplot,3));    
    xmin = g.tlim(1);
    xmax = g.tlim(2);
end

% plot data
% ---------
set(0, 'CurrentFigure', h);
plottopo( erptoplot, 'chanlocs', chanlocs, 'frames', pnts, ...
          'limits', [xmin xmax 0 0]*1000, 'title', g.title, 'colors', colors, ...
          'chans', g.chans, 'legend', legend, 'regions', regions, 'ylim', g.ylim, g.tplotopt{:});

% outputs
% -------
times  = linspace(xmin, xmax, pnts);
erp1   = mean(erp1ind,3);
if length(datsub) > 0 % dataset to subtract
    erp2   = mean(erp2ind,3);
    erpsub = erp1-erp2;
else
    erp2 = [];
    erpsub = [];
end

if nargin < 3 && nargout == 1
    erp1 = sprintf('pop_comperp( ALLEEG, %d, %s);', flag, vararg2str({ datadd datsub options{:} }) );
end
return;

% convert significance values to alpha
% ------------------------------------
function regions = p2regions( pvalues, alpha, limits);
    
    for index = 1:size(pvalues,1)
        signif = diff([1 pvalues(index,:) 1] < alpha);
        pos    = find([signif] > 0);
        pos    = pos/length(pvalues)*(limits(2) - limits(1))+limits(1);
        neg    = find([signif(2:end)] < 0);
        neg    = neg/length(pvalues)*(limits(2) - limits(1))+limits(1);
        if length(pos) ~= length(neg), signif, pos, neg, error('Region error'); end
        regions{index} = [neg;pos];
    end
    
% process data
% ------------
function [erptoplot, erpstd, colors, colstd, legend] = preparedata( erpind, plotavg, plotstd, plotall, mode, tag, dataset, coloravg, allcolors);

    colors    = {};
    legend    = {};
    erptoplot = [];
    erpstd    = [];
    colstd    = {};

    % plot individual differences
    % ---------------------------
    if strcmpi(plotall, 'on')
        erptoplot = [ erptoplot erpind(:,:) ];
        for index=1:size(erpind,3)
            if iscell(dataset)
                if strcmpi(tag, 'Diff ')
                    legend = { legend{:} [ dataset{1}{index} ' - ' dataset{2}{index} ] };
                else
                    legend = { legend{:} dataset{index} };
                end
            else
                legend = { legend{:} [ 'Dataset ' int2str(dataset(index)) ] };
            end
            colors = { colors{:}  allcolors{index} };
        end
    end
    
    % plot average
    % ------------
    if strcmpi( plotavg, 'on')
        if strcmpi(mode, 'ave')
             granderp    = mean(erpind,3);
             legend      = { legend{:} [ tag 'Average' ] };
        else granderp    = sqrt(mean(erpind.^2,3));
             legend      = { legend{:} [ tag 'RMS' ] };
        end
        colors    = { colors{:}  {coloravg;'linewidth';2 }};
        erptoplot = [ erptoplot granderp];
    end

    % plot standard deviation
    % -----------------------
    if strcmpi(plotstd, 'on')
        if strcmpi(plotavg, 'on')
            std1      = std(erpind, [], 3);
            erptoplot = [ erptoplot granderp+std1 ];
            erpstd    = granderp-std1;
            legend    = { legend{:} [ tag 'Standard dev.' ] };
            colors    = { colors{:} { coloravg;'linestyle';':' } };
            colstd    = { { coloravg 'linestyle' ':' } };
        else 
            disp('Warning: cannot show standard deviation without showing average');
        end
    end

% ------------------------------------------------------------------
    
function [p, t, df] = pttest(d1, d2, dim)
%PTTEST Student's paired t-test.
%       PTTEST(X1, X2) gives the probability that Student's t
%       calculated on paired data X1 and X2 is higher than
%       observed, i.e. the "significance" level. This is used
%       to test whether two paired samples have significantly
%       different means.
%       [P, T] = PTTEST(X1, X2) gives this probability P and the
%       value of Student's t in T. The smaller P is, the more
%       significant the difference between the means.
%       E.g. if P = 0.05 or 0.01, it is very likely that the
%       two sets are sampled from distributions with different
%       means.
%
%       This works for PAIRED SAMPLES, i.e. when elements of X1
%       and X2 correspond one-on-one somehow.
%       E.g. residuals of two models on the same data.
%       Ref: Press et al. 1992. Numerical recipes in C. 14.2, Cambridge.

if size(d1,dim) ~= size(d2, dim)
   error('PTTEST: paired samples must have the same number of elements !')
end
if size(d1,dim) == 1
    close; error('Cannot compute paired t-test for a single ERP difference')
end; 
a1 = mean(d1, dim);
a2 = mean(d2, dim);
v1 = std(d1, [], dim).^2;
v2 = std(d2, [], dim).^2;
n1 = size(d1,dim);
df = n1 - 1;
disp(['Computing t-values, df:' int2str(df) ]);

d1 = d1-repmat(a1, [ones(1,dim-1) size(d1,3)]);
d2 = d2-repmat(a2, [ones(1,dim-1) size(d2,3)]);
%cab = (x1 - a1)' * (x2 - a2) / (n1 - 1);
cab = sum(d1.*d2,3)/(n1-1);
% use abs to avoid numerical errors for very similar data
% for which v1+v2-2cab may be close to 0.
t = (a1 - a2) ./ sqrt(abs(v1 + v2 - 2 * cab) / n1) ;
p = betainc( df ./ (df + t.*t), df/2, 0.5) ;

% ------------------------------------------------------------------

function [p, t] = ttest(d1, d2, dim)
%TTEST Student's t-test for equal variances.
%       TTEST(X1, X2) gives the probability that Student's t
%       calculated on data X1 and X2, sampled from distributions
%       with the same variance, is higher than observed, i.e.
%       the "significance" level. This is used to test whether
%       two sample have significantly different means.
%       [P, T] = TTEST(X1, X2) gives this probability P and the
%       value of Student's t in T. The smaller P is, the more
%       significant the difference between the means.
%       E.g. if P = 0.05 or 0.01, it is very likely that the
%       two sets are sampled from distributions with different
%       means.
%
%       This works if the samples are drawn from distributions with
%       the SAME VARIANCE. Otherwise, use UTTEST.
%
%See also: UTTEST, PTTEST.
if size(d1,dim) == 1
    close; error('Cannot compute t-test for a single ERP')
end; 
a1 = mean(d1, dim);
v1 = std(d1, [], dim).^2;
n1 = size(d1,dim);
if length(d2) == 1 && d2 == 0
    a2 = 0;
    n2 = n1;
    df = n1 + n2 - 2;
    pvar = (2*(n1 - 1) * v1) / df ;
else 
    a2 = mean(d2, dim);
    v2 = std(d2, [], dim).^2;
    n2 = size(d2,dim);
    df = n1 + n2 - 2;
    pvar = ((n1 - 1) * v1 + (n2 - 1) * v2) / df ;     
end
disp(['Computing t-values, df:' int2str(df) ]);

t = (a1 - a2) ./ sqrt( pvar * (1/n1 + 1/n2)) ;
p = betainc( df ./ (df + t.*t), df/2, 0.5) ;


