function [Mout,Fout,corrcoef]=FECG_extract_multichan(Fin,fs,peaks,L)

[nleads,nsamples]=size(Fin);
phase = PhaseCalculation(peaks);     % phase calculation

JJ = find(peaks);
fm = fs./diff(JJ);          % heart-rate
w = mean(2*pi*fm);          % average heart-rate in rads.
wsd = std(2*pi*fm,1);       % heart-rate standard deviation in rads.


%%
% Pi-ICA analysis
[indep_comp,W,A] = PiCA([Fin -Fin],[peaks peaks]);

indep_comp = indep_comp(:,1:nsamples);
num_comp=size(indep_comp,1);


%%
% MECG Extraction

bins = 250; % number of phase bins
ECGsd= zeros(L,bins);
nk= zeros(1,L); % a vector indicating number of Gaussian kernels for each independent component
Q= (wsd)^2; % the process noise covariance matrix initialized with its first component. Note: size of Q changes in the for loop.
Wmean= w; % the process noise mean vector initialized with its first component. Note: size of Wmean changes in the for loop.

for comp_ind=1:L
    [ECGmean,ECGsd(comp_ind,:),meanphase] = MeanECGExtraction(indep_comp(comp_ind,:),phase,bins,1); % mean ECG extraction
    ECGsmooth=smooth(ECGmean,6,'moving')';
     
     test1=diff(ECGsmooth);
     test1=sign(test1);
     
     test2=filter([1 1],1,test1);
     gauss_ind=find(test2==0);
     lowchange=abs(diff(ECGsmooth(gauss_ind)))<1e-3; % ECGmean replaced by ECGsmooth 
     gauss_ind=setdiff(gauss_ind,gauss_ind(find(lowchange==1)));
     
    
    theta = meanphase(gauss_ind);
    alpha = 1.2*ECGsmooth(gauss_ind);% ECGmean replaced by ECGsmooth
    b = .04*ones(size(alpha));
    
    InitParams1= [alpha b theta];

    options = optimset('TolX',1e-4,'TolFun',1e-4,'MaxIter',100);
    OptimumParams = nlinfit(meanphase,ECGmean,@ECGModel,InitParams1,options);
    
%         Model = ECGModelError(OptimumParams,ECGmean,meanphase,1);
% 
%     figure;
%     plot(Model,'k','LineWidth',2);
%     hold on
%     plot(ECGmean,'r');
%     plot(ECGsmooth);
%     plot(gauss_ind,ECGsmooth(gauss_ind),'*');
    
    
 
% KF modeling and Parameter setting
    N = length(OptimumParams)/3;% number of Gaussian kernels
    nk(comp_ind)= N;
    subQ= [(.1*OptimumParams(1:N)).^2; (.05*ones(1,N)).^2; (.05*ones(1,N)).^2];
    Q=[Q zeros(size(Q,1),N*3); zeros(N*3, size(Q,2)) diag(subQ(:))]; 
    kp= reshape(OptimumParams,[],3)'; % Kernel parameters reshaped into a matrix. 1st row: alpha, 2nd row: b, 3rd row: theta
    kp= [1 0 0;0 0 1; 0 1 0]*kp; % Kernel parameters reshaped into a matrix. 1st row: alpha, 2nd row: theta, 3rd row: b
    Wmean= [Wmean; kp(:)];    
end

Y = [phase ; indep_comp(1:L,:)]; % Observation matrix
X0 = [-pi zeros(1,L)]'; % mean of the initial state
P0 = diag([(2*pi)^2; (10*max(abs(indep_comp(1:L,:)),[],2)).^2]);
Q=[Q zeros(size(Q,1),L); zeros(L, size(Q,2)) diag((.05*mean(ECGsd(:,1:round(bins/10)),2)).^2)];
R = [(w/fs).^2/12 zeros(1,L) ;zeros(L,1) diag((mean(ECGsd(:,1:round(bins/10)),2)).^2)];
Wmean= [Wmean; zeros(L,1)];
Vmean= zeros(L+1,1);

% Reserve space for estimates.
Xminus = zeros(size(X0,1),size(Y,2));
Xekf = zeros(size(X0,1),size(Y,2));
Pminus = zeros(size(X0,1),size(X0,1),size(Y,2));
Pekf = zeros(size(X0,1),size(X0,1),size(Y,2));

X= X0;
P= P0;

params={Wmean, nk, fs};

% Estimate with EKF
for k=1:size(Y,2)

   
   [X,P] = ekf_predict1(X,P,ekf_df_dx(X,params),Q,ekf_cost(X,params), ekf_df_dw(X,params));
   
    if(abs(X(1)-Y(1,k)) > pi)
       X(1) = Y(1,k);
    end 
   
   Xminus(:,k)   = X;
   Pminus(:,:,k) = P;
   
   [X,P] = ekf_update1(X,P,Y(:,k),eye(L+1),R);
   Xekf(:,k)   = X;
   Pekf(:,:,k) = P;
     
end



[Xeks,Peks] = ekf_smooth1(Xekf,Pekf,Xminus,Pminus,@ekf_df_dx,params);

indep_comp_cleaned= indep_comp(1:L,:)-Xeks(2:L+1,:);
MECG_comp_cleaned= Xeks(2:L+1,:);



%%
% Inverse Periodic-ICA

FICA=[indep_comp_cleaned;indep_comp(L+1:4,:)];
MICA=[MECG_comp_cleaned;zeros(num_comp-L,nsamples)];
Fout=A*FICA; 
Mout=A*MICA;

corrcoef = zeros(4,1);
for i=1:4
    xx=dualSample(Fout(i,:),peaks,phase);
    corrcoef(i)=mean(xx.*Fout(i,1:length(xx)))/sqrt(mean(xx.^2)*mean(Fout(i,1:length(xx)).^2));
end

end