disp('   ');
disp(' ttl_cps.m  ver 1.1  May 1, 2014 ');
disp(' by Tom Irvine');
disp('  ');
disp(' The input file must have two columns:  time(sec) & TTL ');
disp('  ');
%
clear t;
clear ts;
clear ff;
clear y;
clear cps;
clear length;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
disp('   3=Excel file ');
file_choice = input('');
%
if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        THM=THM';
end
if(file_choice==2)
        FS = input(' Enter the matrix name:  ','s');
        THM=evalin('caller',FS);
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
t=THM(:,1);
y=THM(:,2);
%
n=length(t);
dt=(t(n)-t(1))/(n-1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
	vel(1)=( -y(3)+4.*y(2)-3.*y(1) )/(2.*dt);
	vel(2)=( -y(4)+4.*y(3)-3.*y(2) )/(2.*dt);
%
	ddt=12.*dt;
%
    vel(3:(n-2))=(-y(5:n)+8*y(4:(n-1))-8*y(2:(n-3))+y(1:(n-4)))/ddt;
%
    vel(n-1)=( +y(n-1)-y(n-3) )/(2.*dt);
    vel(n)  =( +y(n-1)-y(n-2) )/dt;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
threshold=0.5*max(abs(vel));
%
j=1;
%
pflag=0;
%
for i=2:n
    if( vel(i)>threshold && pflag==0)
        ts(j)=t(i);
        pflag=1;
        j=j+1;
    end    
%
    if( vel(i)<-threshold && pflag==1)
        pflag=0;
    end
%
end
%
num=j-1;
%
ff=zeros(num,1);
%
for i=1:(num-1)
    ff(i)=1/(ts(i+1)-ts(i));
end
ff(num)=ff(num-1);
%
cps=zeros(n,1);
%
jj=1;
for i=1:n
    if(t(i)<=ts(1))
        cps(i)=ff(1);
    end
    if(t(i)>=ts(num))
        cps(i)=ff(num);
    end
    if(t(i)>ts(1) && t(i)<ts(num))
        for j=jj:(num-1)
            if(t(i)>=ts(j) && t(i)<ts(j+1))
                cps(i)=ff(j);
                jj=j;
                break;
            end
        end
    end    
end
%
figure(1);
plot(t,y);
title('TTL');
xlabel('Time(sec)');
grid on;
%
figure(2);
plot(t,cps);
title('CPS')
ylabel('Hz');
xlabel('Time(sec)');
grid on;
%
disp(' Output:  cps_array ');
cps_array=[t cps];
%
min_cps=min(cps);
mean_cps=mean(cps);
max_cps=max(cps);
%
disp(' ');
out1=sprintf(' max  = %8.4g Hz ',max_cps);
out2=sprintf(' mean = %8.4g Hz ',mean_cps);
out3=sprintf(' min  = %8.4g Hz ',min_cps);
disp(out1);
disp(out2);
disp(out3);
disp(' ');
out1=sprintf(' max  = %8.4g RPM ',max_cps*60);
out2=sprintf(' mean = %8.4g RPM ',mean_cps*60);
out3=sprintf(' min  = %8.4g RPM ',min_cps*60);
disp(out1);
disp(out2);
disp(out3);