function [s,B,SpikeTrain] = FastICA(EMG,M,Fs)
% EMG (N * M): M is the length of signal. N is the number of channels
% M: number of total iteration
% Fs: sample frequency
% s: separated source after decomposition
% B: separation vector
% SpikeTrain: motor unit spike trains after decomposition
% Default Settings
Tolx = 10^-4;
[NumCh,N] = size(EMG);
s = zeros(N,M);
stemp = zeros(N,M);
B = zeros(NumCh,1);
SpikeTrain = zeros(N,M);

% W = eye(NumCh);
% [b,a] = butter(4,100/(Fs/2),'low');
% Gx = x^3/3; gx = x^2;; gx' = 2x;
for i = 1:M
    w = [];
    %     w(:,1) = zeros(NumCh,1);
    %     w(:,2) = W(:,i);
    w(:,1) = randn(NumCh,1);
    w(:,2) = randn(NumCh,1);
    for n = 2:100
        if abs(w(:,n)'*w(:,n-1)-1)>Tolx
            A = mean(2*w(:,n)'*EMG);
            w(:,n+1) = EMG*(((w(:,n)'*EMG)').^2)-A*w(:,n);
            w(:,n+1) = w(:,n+1) - B*B'*w(:,n+1);
            w(:,n+1) = w(:,n+1)/norm(w(:,n+1));
        else
            break;
        end
    end
    CoV(1) = 1;
    CoV(2) = 0.99;
    %     for m = n:100
    %         if abs(CoV(m-n+2)-CoV(m-n+1))>Tolx
%     stemp(:,i) = w(:,n)'*EMG;
%     s(:,i) = filtfilt(b,a,stemp(:,i));
    s(:,i) = w(:,n)'*EMG;
    [pks,loc] = findpeaks(s(:,i).^2);
    
    
    [idx,C] = kmeansplus(pks',2);
%     [idx,~] = myCluster2(pks);
    
    
    if sum(idx==1)<=sum(idx==2)
        SpikeLoc = loc(idx==1);
    else
        SpikeLoc = loc(idx==2);
    end
    %             SpikeInterval = diff(SpikeLoc);
    %             CoV(m-n+3) = std(SpikeInterval)/mean(SpikeInterval);
    %             w(:,m+1) = mean(EMG(:,SpikeLoc),2);
    %         else
    %             break;
    %         end
    %     end
    SpikeTrain(SpikeLoc,i) = 1;
    B(:,i) = w(:,end);
end