disp(' ');  
disp(' differ.m, ver 1.3,  March 14, 2013');
disp('  by Tom Irvine ');
disp('  Email:  tom@vibrationdata.com ');
disp(' ');
disp(' This program differentiates a time history.');
%
clear THM;
clear vel;
clear t;
clear y;
clear v;
clear vv;
clear c;
clear a;
clear diff_data;
%
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)
        THM = input(' Enter the matrix name:  ');
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
t=THM(:,1);
y=THM(:,2);
%
tmx=max(t);
tmi=min(t);
n = max(size(y));
a=zeros(n,1);
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
disp(' ')
disp(' Time Step ');
clear difft;
difft=diff(t);
dtmin=min(difft);
dtmax=max(difft);
%
out4 = sprintf(' dtmin  = %8.4g sec  ',dtmin);
out5 = sprintf(' dt     = %8.4g sec  ',dt);
out6 = sprintf(' dtmax  = %8.4g sec  ',dtmax);
disp(out4)
disp(out5)
disp(out6)
%
disp(' ')
disp(' Sample Rate ');
out4 = sprintf(' srmin  = %8.4g samples/sec  ',1/dtmax);
out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
out6 = sprintf(' srmax  = %8.4g samples/sec  \n',1/dtmin);
disp(out4)
disp(out5)
disp(out6)
%
ncontinue=1;
if(((dtmax-dtmin)/dt)>0.01)
    disp(' ')
    disp(' Warning:  time step is not constant.  Continue calculation? 1=yes 2=no ')
    ncontinue=input(' ');
end
if(ncontinue==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;
%  	
%%	for(i=3:(n-2))   
%%        vel(i)=( -y(i+2) +8.*y(i+1) -8.*y(i-1) +y(i-2) ) / ddt;
%%  end
%
    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;
end    
%
disp(' Select input data type: ');
disp('   1=displacement  ');   
disp('   2=velocity      ');
disp('   3=acceleration  ');
disp('   4=other         ');
idt=input(' ');
%
figure(1);
plot(t,y);
if(idt==1)
   title('Displacement');   
end
if(idt==2)
   title('Velocity');    
end
if(idt==3)
   title('Acceleration'); 
end
if(idt==4)
   title('Original'); 
end
%
xlabel('Time(sec)');
grid on;
figure(2);
plot(t,vel);
if(idt==1)
   title('Velocity');   
end
if(idt==2)
   title('Acceleration');    
end
if(idt==3)
   title('Jerk'); 
end
if(idt==4)
   title('Differentiated'); 
end
xlabel('Time(sec)');
grid on;
diff_data=[t, vel'];
disp(' The matrix name is:  diff_data ');