function [MPSRF] = GRmpsrf(X,T,P,Nch)
W = zeros(P); B = zeros(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 + xw'*xw;  B = B + xb'*xb;
end
W = W / ((T-1) * Nch); B = B / (Nch-1); V = sort(abs(eig(W\B)));
MPSRF = sqrt( (T-1)/T + V(end) * (Nch+1)/Nch); end

