function Ke=electrodeKernel(K,startTail)
% Ke=electrodeKernel(K,startTail)
%
% Extracts the electrode kernel Ke from the raw kernel K
% by removing the membrane kernel, estimated from the
% indexes >= startTail of the raw kernel.
%
% K-Ke is the membrane kernel filtered by the electrode
%
% R. Brette, Sep 2005

    function [err,Kel]=removeKm(RawK,x,tau,dt,tail)
    % Solves Ke = RawK - Km * Ke/Re
    % for an exponential Km
    % Returns the error on the tail
    
        alpha=x*dt/tau;
        lambda=exp(-dt/tau);
        Y=0*RawK;
        Y(1)=alpha/(alpha+1)*RawK(1); % K(1) == 0 ?
        for i=2:length(RawK)
            Y(i)=(alpha*RawK(i)+lambda*Y(i-1))/(1+alpha);
        end
        Kel=RawK-Y;
        err=Kel(tail)'*Kel(tail);
    end

    dt=0.1;

    % Exponential fit of the tail
    mymodel=fittype('a*exp(b*x)');
    tail=startTail:length(K);
    t=((0:(length(K)-1))*dt)';
    rfit=fit(t(tail),K(tail),mymodel,'start',[1,-0.2]);
    % Eliminate the offset (??)
    %K=K-rfit.c;

    % Membrane time constant and resistance
    tau = -1/rfit.b
    R = rfit.a*tau/dt - sum(rfit.a*exp(rfit.b*t(1:2)));
    % NB: the first two steps are ignored

    % Electrode resistance
    Re = sum(K(3:startTail))-sum(rfit.a*exp(rfit.b*t(3:startTail)));

    % Replace tail by fit
    K(tail)=feval(rfit,t(tail));

    % Clean the beginning (optional)
    K(1:2)=0;

    Km=feval(rfit,t);
%    Km(tail)=K(tail);

    % Optimize the membrane kernel
    z=fminbnd(@(x)removeKm(K,x,tau,dt,tail),0.5*R/Re,2*R/Re);
    [err,Ke]=removeKm(K,z,tau,dt,tail);

    % Clean the electrode kernel (remove negative numbers)
    %Ke=Ke.*(Ke>0);

    % Display electrode resistance
    Re=sum(Ke)
    R=z*Re
    
end
