function [PSRF] =  GRpsrf(X,T,P,Nch)
W = zeros(1,P); B = zeros(1,P); mn = mean(reshape(mean(X),P,Nch)');
for n=1:Nch
  xw = X(:,:,n) - repmat(mean(X(:,:,n)),T,1);   xb = mean(X(:,:,n)) - mn;
  W = W + sum(xw.*xw);  B = B + xb.*xb; end
W = W / ((T-1) * Nch); S = (T-1)/T * W + B/(Nch-1); 
PSRF = sqrt((Nch+1)/Nch * S ./ W - (T-1)/Nch/T)
end
