count_final=lc_fast(signal,levels)

function count_final=lc_fast(signal,levels)
%
%
%LC_FAST calculates how many times a given input signal crosses a set of levels.
% Input variable "signal" can be also a mxn matrix and in that case LC_FAST calculation is done for each column.
% Input variable "levels" must be a column vector qx1.
% Output variable "count_final" is a qxn array where count(i,j) is the number of crossings between signal(:,j) and levels(i).
%
%
% Author: Roberto Madia
% Version 1.0 made on Nov 2009
%
%
%
% Comments:
% This algorithm is generally faster than traditional ones and it has been specifically developed to handle big dimension datasets.
% Example:
% data=rand(200000,12); 12 random signals made by 200,000 samples each.
% levels=rand(1000,1); 1000 random levels to be checked.
% count_final=lc_fast(data,levels); % computation done few seconds!
%
%
% Technical notes: 
% This routines calculates the number of crossings considering also the cases where the signal reaches a given level without going
% above it (just level touching).
%
%


cells=length(levels);
[n,speeds]=size(signal);
count=zeros(cells+1,speeds);
EPS=1e-13;
[levels2,ind_lev]=sort(levels);
signal=[ones(1,speeds)*EPS; signal; ones(1,speeds)*EPS];
% calculation of exact crossings and dataset preparation for next histogram
for i=1:cells
    ind=find(signal==levels2(i));
    if not(isempty(ind))
        signal(ind(signal(ind-1)>=levels2(i)))=levels2(i)+EPS;
        signal(ind(signal(ind-1)<=levels2(i)))=levels2(i)-EPS;
        ind11=ind(signal(ind-1)>levels2(i)&signal(ind+1)<levels2(i));
        ind22=ind(signal(ind-1)<levels2(i)&signal(ind+1)>levels2(i));
        ind_temp=[ind11;ind22];
        [r,c]=ind2sub([n+2 speeds],ind);
        for j=1:length(c)
            count(i,c(j))=count(i,c(j))+1;
        end
        if not(isempty(ind_temp))
            [r2,c2]=ind2sub([n+2 speeds],ind_temp);
            for j=1:length(c2)
                count(i,c2(j))=count(i,c2(j))-1;
            end
        end
    end
end


[cc,bin]=histc(signal(2:end-1,:),[-inf; levels2; inf]);
% calculation of crossings
for i=1:n-1
    for j=1:speeds
        increment=bin(i+1,j)-bin(i,j);
        if sign(increment)==-1
            count(bin(i+1,j):bin(i,j)-1,j)=count(bin(i+1,j):bin(i,j)-1,j)+1;
        elseif sign(increment)==1
            count(bin(i,j):bin(i+1,j)-1,j)=count(bin(i,j):bin(i+1,j)-1,j)+1;
        end
    end
end
count_final=count(ind_lev,:);