Can anyone explain this adaptive kernel density estimation method?

Question

I happened to need a reliable multivariete kernel density estimation method, and found it in this Matlab code, which turned out to work very well with my data. But I have trouble understanding its workings. Author provides no explanation (he links a paper, but unless I badly misunderstood it, it does not concern said method) nor have I been able to find it elsewhere.

I was hoping someone might recognize it / figure it out and explain it / direct me to resources explaining it.

What I do / don't understand

Once estimation parameters (set of kernel locations μ, weights w and squared bandwidths σ; for the sake of simplicity, I will stick with univariete case) are determined by iteratively evolving initial guess, estimate is provided simply as $$ \widehat{f}(x) = \sum_{i=1}^{\gamma} {\frac{w_i}{\sqrt{2 \pi \sigma_i}} e^{-\frac{(x - \mu_i )^2}{2 \sigma_i }}} $$ where $ \gamma $ is number of kernels used.

The procedure of updating parameters (univariete case, Matlab syntax) can be simplified without a loss of functionality as follows :

function [w,mu,Sig,bandwidth,ent]=regEM(w,mu,Sig,bandwidth,X)
[gam,d]=size(mu);[n,d]=size(X);
p=zeros(n,gam); psig=p;
for i=1:gam
    xRinv = (X - mu(i)).^2 / Sig(i);
    xSig = (xRinv/s)+eps;
    p(:,i) = exp(-0.5 * xRinv - 0.5 * bandwidth^2 / Sig(i)) * w(i) / sqrt(2 * pi * Sig(i));
    psig(:,i) = p(:,i) .* xSig;
end
density = sum(p,2); psigd = sum(psig,2);
p = bsxfun(@rdivide, p, density); % normalize classification prob.
ent = sum(log(density)); w = sum(p,1);
for i=find(w>0)
    mu(i) = p(:,i)'*X/w(i); %compute mu's
    Sig(i) = p(:,i)'*((X - mu(i)).^2)/w(i)+bandwidth^2; % compute sigmas
end
w=w/sum(w);
curv = mean(psigd ./ density);
bandwidth = 1/(4*n*(4*pi)^(1/2)*curv)^(1/3);
end

I see how this resembles finding an optimal bandwidth that minimizes mean integrated square error in basic KDE (as described in this paper). But I don't get how exactly does this work here (deriving expression for MISE is much harder in this case). I don't know how does the 'bandwidth' variable come to play here, how to interpret the additional apperances of it, etc.

Recapitulation

Could anyone help me understand motivation behind this method and explain it's workflow?