// Active Electrode Compensation
//
// R. Brette (brette@di.ens.fr)
// Apr 2006
//

function [vi,ii]=correlations(v,i,ksize)
// [vi,ii]=correlations(v,i,ksize)
//
// Calculates the correlation vector <v(n)i(n-k)>
// and the autocorrelation vector <i(n)i(n-k)>
// from the recording v and the current i.
// ksize the size of the full kernel (typically 150-200 points)

  vi=zeros(ksize,1);
  ii=vi;
  v=v-mean(v);
  i=i-mean(i);
  for k=1:ksize
    vi(k)=sum(v(k:$).*i(1:$-k+1));
    ii(k)=sum(i(k:$).*i(1:$-k+1));
  end
endfunction

function K=rawKernel(vi,ii)
// K=rawKernel(VI,II)
//
// Calculates the raw kernel from the autocorrelation of the input current
// II and correlation between current and recorded potential VI
// i.e. finds the best linear filter which transforms I into V
//
// R. Brette, Sep 2005

  K=toeplitz(ii)\vi;
endfunction

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);
    endfunction

    dt=0.1;

    // Exponential fit of the tail
    deff('e=mymodel(p,z)','e=p(1)*exp(p(2)*z(1))-z(2)');
    tail=startTail:length(K);
    t=((0:(length(K)-1))*dt)';
    p=datafit(mymodel,[t(tail);K(tail)],[1;-0.2]);

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

    // Electrode resistance
    Re = sum(K(3:startTail))-sum(p(1)*exp(p(2)*t(3:startTail)));

    // Replace tail by fit
    deff('y=fitexp(x)','y=p(1)*exp(p(2)*x)');
    K(tail)=feval(t(tail),fitexp);

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

    Km=feval(t,fitexp);
//    Km(tail)=K(tail);

    // Optimize the membrane kernel
    deff('e=remK(x)','e=removeKm(K,x,tau,dt,tail)');
    z=fminbnd(remK,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;
    printf("Re = %f\tR=%f",Re,R);
    
endfunction