function [ A_est, S_est, time_Proposed, H_est, Y_normalized] = ProSPA(R_whiten,K,B,X)
%Input:
%R_whiten: Local Covariances (whitened)
%K: number of sources
%L: frame length
%B: prewhiten matrix
%X: data
%output:
% A_est: channel estimation
% S_est: source estimation
% H_est = KR(conj(A_est),A_est)
% Y_normalized: for initializing VolMin-AO


% copyrights: Xiao Fu
% xfu@umn.edu
% See: Xiao Fu, Wing-Kin Ma, Kejun Huang, and Nicholas Sidiropoulos, 
% "Blind Separation of Quasi-Stationary Sources: Exploiting Convex Geometry
% in Covariance Domain" submitted to IEEE Trans. Signal Process., 2014



%--- ------------------ ProSPA-------------------------

tic;
[K,K,M]=size(R_whiten);
z = R_whiten(:);
Y_hat_whiten_unnorm = reshape(z,K^2,M);
Runnorm = (1/M)*Y_hat_whiten_unnorm*Y_hat_whiten_unnorm';

[Q_unnorm E_unnorm] = eig(0.5*(Runnorm +Runnorm'));
Us_unnorm = Q_unnorm(:,end-K+1:end);
%-------------
Y_hat_whiten_unnorm = Us_unnorm*Us_unnorm'*Y_hat_whiten_unnorm;
Ind = bsxfun(@plus,((1:K)+K*(0:K-1))',K*K*(0:M-1));
TR = sum( R_whiten(Ind));
Y_hat_whiten_unnorm = bsxfun(@times,Y_hat_whiten_unnorm,1./TR);
[ H_est, A_est1 ] = SuccessiveFeasble( Y_hat_whiten_unnorm,K );
A_est= B*A_est1 ;
time_Proposed = toc;

Y_normalized = Y_hat_whiten_unnorm;
S_est = pinv(A_est)*X; % simple zero-forcing unmixing
end

function [ H, A ,index] = SuccessiveFeasble( Y,K )
% The program is for successive feasibility implementation
% Input:
% Y: dataset Y
% K: the order of the model
% R: the local covariance matrices
% Us: The subspace projection matrix, which is optional
% Output:
% H: the estimation of the KR(A,A)
% A: the estimation of the mixing matrix
% Programmed by Fu Xiao
% xfu@ee.cuhk.edu.hk

[n,m]=size(Y); 
H_est=[]; 
index = []; 
P = eye(n);                         
for k=1:K
    [~, ind]=max(sum(abs(P*Y).^2).^(1/2));
    H_est = [H_est Y(:,ind)];
    A_temp = reshape(H_est(:,k),sqrt(n),sqrt(n));
    [U1 S1] = eig((A_temp+A_temp')/2);
    A_est(:,k) = U1(:,end);
    P = eye(n) - H_est*pinv(H_est);
    index = [index ind];
end




H = H_est;
A = A_est;






end

