
function [f, g] = refavg(u, ilayer, n_iterations)
% function [f, g] = refavg(u, ilayer, n_iterations)

[nd, nt] = size(u);

uu = u;
f = zeros(nd, 1);
g = zeros(1,nt);

for it = 1:n_iterations
   [ff, gg] = ref_est(uu, ilayer);
   %fprintf('%f %f %f %f\n', mean(ff), std(ff), mean(gg), std(gg))
   uu = uu - ff(:, ones(1, nt)) - gg(ones(1, nd), :);
   f = f + ff;
   g = g + gg;
end


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

function [f, g] = ref_est(u, ilayer);
nl = length(ilayer);
[nd, nt] = size(u);

[f, sf, nf] = mstdgap(u, 2);
f = f.' - mean(f(ilayer));

[g, sg, ng] = mstdgap(u(ilayer, :) - f(ilayer, ones(1, nt)));

%% This is slow; it can be sped up by using a simplified
%% local version of mstdgap, and probably eliminating the
%% function calls.  The gaps and numbers of points stay the
%% same for each iteration--no need to keep recalculating,
%% (Asking for the medians in mstdgap makes it much slower.)



