function [c, A, fr, U, e, ccum] = multicoherence(y, X, npp, fs)

    %% compute multicoherence and coefficient estimates
    N = size(X,1);
    % compute cross-spectral-densities
    SX = zeros(N,N,npp/2+1);
    SY = zeros(npp/2+1,N);
    fprintf(2, 'Computing CSD matrix\n');
    for i=1:N
        fprintf(2, '\r%d/%d     ', i, N);
        SY(:,i) = cpsd(X(i,:), y, hanning(npp), npp/2, npp, fs);
        for j=i:N
            SX(i,j,:) = cpsd(X(j,:), X(i,:), hanning(npp), npp/2, npp, fs);
            SX(j,i,:) = conj(SX(i,j,:));
        end
    end
    [Syy,fr] = pwelch(y, hanning(npp), npp/2, npp, fs);

    % solve for coherence and coefficients
    fprintf(2, '\nComputing multicoherence\n');
    c = zeros(npp/2+1,1);
    A = zeros(N, npp/2+1);
    U = zeros(N,N,npp/2+1);
    e = zeros(N,npp/2+1);
    ccum = zeros(N,npp/2);
    for i=1:npp/2+1
        %A(:,i) = squeeze(SX(:,:,i)) \ SY(i,:)';
        A(:,i) = pinv(squeeze(SX(:,:,i))) * SY(i,:)';
        c(i) = (SY(i,:) * A(:,i)) / Syy(i);
        [U(:,:,i), ee] = eig(squeeze(SX(:,:,i)));
        e(:,i) = diag(ee);
        ccum(:,i) = abs(SY(i,:)*squeeze(U(:,:,i))).^2 ./ e(:,i).' / Syy(i);
    end
    % remove small imaginary part due to numerical errors
    c = real(c);
    
end
