function [p,t,df]=pc1(Data,Data2)
% [p,t,df]=pc1g(Data);  or   [p,t,df]=pc1g(Data1, Data2);
% Data has any number of dimensions, but the last one (slowest moving index in MATLAB) in the number of cases
% Data is most often the case by case difference between two conditions that are repeated on each case) but can be any matrix of
% independently replicated data on which to test whether they contain signal or only noise
% If Data hase more than two dimensions, it is first made two dimensional by merging all but its last dimension, which is equivalent
% to testing, in a single statistical test, the null hypothesis that there is no signal, e.g., at any channel over the epoch provided
% (but including channels that do not express the difference reduces the power to detect a difference overall). 
% Returns the two-tailed probability p, a Student t value for paired data and the associated degrees of freedom.
% If two arguments are provided, their difference is calculated into Data1 and the test proceeds as if only one data matrix was presented
% The test has the null hypothesis that the data contain only noise (no systematic pattern that
% causes systematic similarity between independent replications of an ERP condition or of the difference between two conditions.
% If there is such a pattern, it could be evaluated as the mean of all replications, but using this as a template would insure that
% the mean similarity of individual cases with the mean is higher than 0, leading to a biased test.
% Instead, the first principal component in used as a template; the similarity with the template is quantified simply as
% the sum of products of each case with the template and a Student t test checks whether the mean similarity differs from 0.
% The algorithm, however, directly calculates scores proportional to these weights without extracting the first PC of the set of waves. 
%
% this function uses TCDF from the statistical toolbox
%
% reference:
% Achim, A.  Signal detection in averaged evoked potentials:  Monte Carlo comparison of the sensitivity of different methods.
%            Electroencephalography and Clinical Neurophysiology, 1995, 96: 574 584.
%
if nargin>1, Data=Data-Data2; end
k=ndims(Data);
n=size(Data,k);
if k>2, Data=reshape(Data,prod(size(Data))/n,n); end  %force into two diemnsional matrix
df=n-1;
V=extractPC1(Data.'*Data,1e-6);  %preferred to MATLAB's eigs which extract more components and uses a random element that can slightly modify the solution from call to call with same data
m=mean(V);
d=V-m(ones(1,n),:);
s2=d.'*d/(n*df);
if s2<=eps, error('The data provided as argument to PC1 have no variability across their last index'); end
t=sqrt(m*m/s2);
p=2 * tcdf(-t,df);


function [c]=extractPC1(A,converg)
% [c,eigv]=extractPC1(A,converg)
% A is a symmetrical (n,n) matrix and c is (n,1)
% converg (default 1e-6), if spercified is the maximal difference between scores (normalised for the initially largest element =1)
% between two successive iterations, below which iterations terminate
if nargin<2, converg=1e-6; end
[v,r]=max(diag(A));  % r is rank of maximal element in the idagonal
c=A(:,r)/v;
pre=zeros(size(c));
while max(abs(c-pre))>=converg,
   pre=c;
   c=A*c;
   c=c/c(r);
end
t=sqrt(c.'*c);
if mean(c)<0, t=-t; end
c=c/t;
if nargout>1
   pre=A*c;
   eigv=pre(r)/c(r);
end