GEH = cell(nbInstrument,nbPartials);

% Remodel Spectrogram
xHat = zeros(HFFTS,len);
for ins = 1:nbInstrument
    for harm = 1:nbPartials
        GEH{ins,harm} = sparse(zeros(HFFTS, len));
        for frame = first:curFrame
            GEH{ins,harm}(:,frame-first+1) = (filterBasis * filterCoefficients(:,ins)*gains(ins,frame)).*excitationSpectrums{ins,frame}(:,harm);
        end
        xHat = xHat + GEH{ins,harm}*partialsWeights(harm,ins);
    end
end

%Compute distance and r
distPartial(nbUpdates) = sum(sum( tmpMaskedSpectrogram .* log((tmpMaskedSpectrogram+eps) ./ (xHat+eps)) - tmpMaskedSpectrogram + xHat));
r = (currentMag+eps)./(xHat+eps);

%Update partialsWeights
updateGammaCoeff = ones(nbPartials,nbInstrument);
for ins=1:nbInstrument
    for harm = 1:nbPartials
        updateGammaCoeff(harm,ins) = sum(sum(GEH{ins,harm}.*r +eps,1))./sum(sum(GEH{ins,harm} +eps,1));
    end
end
clear GEH
partialsWeights = partialsWeights.*updateGammaCoeff;
for ins=1:nbInstrument
    maxValue = max(partialsWeights(:,ins));
    partialsWeights(:,ins) = partialsWeights(:,ins) / maxValue;
    gains(ins,:) = gains(ins,:) * maxValue;
end
%     figure;stem(partialsWeights), title(sprintf('Partial Weights (iter%d)',nbUpdates));
