function [RNew, phi, psi, alphaNew, gammaOld, gammaNew] = ...
        computeGammaMultipleSinus(...
        ROld, iOrder, crossCorrelationVectors, a, phi, psi, gammaOld, ...
        gammaNew, alphaOld, hankelMatrixIsAdded, dcIsIncluded)

    nPitches = length(a);
    RNew = computeRowsOfToeplitzHankelMatrix(iOrder, iOrder, ...
        crossCorrelationVectors, hankelMatrixIsAdded, dcIsIncluded);
    lambda = a-sum(ROld.*phi,1);
    mu = -sum(ROld.*psi,1);
    phi = [phi;zeros(1,nPitches)]+(ones(iOrder-1,1)*lambda).*gammaNew;
    psi = [psi;zeros(1,nPitches)]+(ones(iOrder-1,1)*mu).*gammaNew;
    alphaNew = sum(RNew(1:end-1,:).*gammaNew,1);
    b = (ones(iOrder-1,1)*(alphaOld-alphaNew)).*gammaNew+...
        [zeros(1,nPitches);gammaNew(1:iOrder-2,:)]+...
        [gammaNew(2:end,:);zeros(1,nPitches)]-...
        [gammaOld(1:iOrder-2,:);zeros(1,nPitches)]+...
        (ones(iOrder-1,1)*psi(end,:)).*phi-...
        (ones(iOrder-1,1)*phi(end,:)).*psi;
    nu = sum(RNew(1:end-1,:).*b)./gammaNew(end,:);
    gammaOld = gammaNew;
    gammaNew = nan(iOrder,nPitches);
    gammaNew(iOrder,:) = 1./(nu+RNew(iOrder,:));
    gammaNew(1:iOrder-1,:) = (ones(iOrder-1,1)*...
        (gammaNew(iOrder,:)./gammaOld(end,:))).*b;

end
