function [A_est,iter] = PAPA_mod(x,K,L,err_bound,iter_max)
%======================time line of coding =======================
% Modefied by Fu Xiao on Nov.17,2012
% The Chinese University of Hong Kong.
% xfu@ee.cuhk.edu.hk

% KaKit, 22Jan 2012.: Slight update

% Modified by Ken on Aug. 14, 2010

% Programmer: Yi-Lin Chiour
% National Tsing Hua University, Hsinchu, Taiwan 
% Date: June. 08, 2010

%======================         Usage    ==================================

% [A_est,iter,time_PAPA] = PAPA(x,K,L,err_bound,iter_max)

%====================== Parameters ========================================
% Outputs: 
% A_est     : the estimated mixing matrix
% iter      : averge iteration for estimating one column 
% time_PAPA : the running time of PAPA
%---------------------------------------------------
% Inputs:
% x          : the data.
% K          : number of sources
% L          : frame length
% err_bound : the tolerance between the currently estimated channel and the previously estimeated channel.
% iter_max  : the maximal number of iteration. 





%------------------------ local covariance compute --------------------
[N,T]=size(x);
if N<K
    error('PAPA is for overdetermined case')
end


M = floor(T/L); % Number of frames
overlap = 1;  % means we want to overlap time blocks by 50%
if overlap == 0;
    M_overlap = M;
elseif overlap == 1;
    M_overlap = 2*M-1; % number of frames after overlapping
end
Rx = zeros(N, N, M_overlap);
for mmm = 1:M-1
    R(:,:,mmm) = x(:, (mmm-1)*L+1:mmm*L)*x(:, (mmm-1)*L+1:mmm*L)'/L;
    R(:,:,M+mmm) = x(:, (mmm-1)*L+1+(L/2):mmm*L+(L/2))*x(:, (mmm-1)*L+1+(L/2):mmm*L+(L/2))'/L;
end
R(:,:,M) = x(:,(M-1)*L+1:M*L)*x(:,(M-1)*L+1:M*L)'/L;


% %--------------------- Noise Removal (you may comment it)--------
sigma_est = zeros(M_overlap,1);
for mmm = 1:M_overlap
    [eigvec, eigval] = eig(R(:,:,mmm));
    if min(diag(eigval)) < sigma_est
        sigma_est(mmm) = min(diag(eigval));
    end
end

for mmm = 1:M_overlap
    R_noisefree(:,:,mmm) = R(:,:,mmm) - min(sigma_est)*eye(N);
end

% %------------------------ Pre-whitening stage (if overdetermined)------
R = R_noisefree;

[N,N,M]=size(R);
z_temp = R(:);
R_reshape = reshape(z_temp,N^2,M);

R_bar_vec = reshape(R,N^2,M);
R_bar_vec = (1/M)*sum(R_bar_vec,2);
R_bar = reshape(R_bar_vec,N,N);

[U S] = svd(0.5*(R_bar+R_bar'));
B = U(:,1:K)*sqrt(S(1:K,1:K));

PinvB=pinv(B);
R_vec = reshape(R,N^2,M);
R_whiten = kron(conj(PinvB),PinvB)*R_vec;

B_pre = B;

tic
%------------------------- Khatri-Rao subspace extraction ------

z = R_whiten(:);
Y = reshape(z,K^2,M);

iter = 0;
C = Y*Y';
[U, S] =eig(C+C');

Us = U(:, end-K+1:end);

for k = 1:K
    %---------------------------Random Initialization ---------------------
    h            = Us*randn(K,1);   
    
    %----------------------- Alternating Projection ----------------------
    [g, g_m, iter_t]   = ALTPRJ(h,iter_max,err_bound,Us,K);  
    iter         = iter+ iter_t;
    G(:,k)       = g;

    
    if k < K
        %--------- Orthogonal Projection: (eye(K^2)-g_m*g_m')*Us------------
        UU = g_m'*Us;
        gUs = g_m*UU;
        Us =  Us - gUs;
%        Us = (eye(K^2)-g_m*g_m')*Us;
       %-----------------(Fu Xiao Nov.13,2012)-----------------------------
    end
    
end

time_PAPA = toc;
iter= iter/K; % to make the comparison fair 

A_est = B_pre*G;


function [a,a_m,iter] = ALTPRJ(h,iter_max,err_bound,Us,K)
%=====================================================================
% Programmer: Yi-Lin Chiour
% National Tsing Hua University, Hsinchu, Taiwan 
% Date: June. 08, 2010

% Modified by Fu Xiao, Nov.13,2012.
% The Chinese University of Hong Kong.
% xfu@ee.cuhk.edu.hk


%======================================================================
% [h] = ALTPRJ(h,iter_max,err_bound,Us,K)
%======================================================================
% Outputs: 
% a         : one column vector of the estimated mixing matrix.
%---------------------------------------------------
% Inputs:

% h         : a vector belongs to the range space of the basis matrix Us.
% iter_max  : the maximal number of iterations. 
% err_bound : for stopping criterion
% Us:         The Khatri-Rao subspace
% K:          number of sources

obj_val_old= inf;



for i=1:iter_max
      C = reshape(h,K,K); 
      C = (C+C')/2; % 
      [eigvec,eigval]   = eig(C);
      [eigval order_eig]= sort(diag(abs(eigval)),'descend');
      eigvec = eigvec(:,order_eig);

      a                 = eigvec(:,1);
      a_m               = kron(conj(a),a);
      alpha             = eigval(1)./abs(eigval(1));

      h                 = alpha* Us*Us'*a_m;
      %---------------------------------------------
      obj_val_cur       = sum(abs(h- alpha*a_m).^2);   
      %----------------Stopping Criterion -------------------------------
      if i > 1,
          if abs(obj_val_old- obj_val_cur)./abs(obj_val_old)  < err_bound,%
            iter= i; return; end;
      end;
      obj_val_old  = obj_val_cur;
end
iter= iter_max;