function [Ax1, Ax2, Ax3,TotalTime,TotalTimeS] = accel(File1,File2,NumberOfFiles,PlotTime,YAxisRange);

% accel.m
% use as [Ax1, Ax2,Ax3,Time,TimeS] = accel(File1,File2,NumberOfFiles,PlotTime,YAxisRange);
% Inputs:   File1 : a string of character containing the name of the first
%                   file to be processed. No extension, in single quotation marks.
%           File2 : a string of character containing the name of the second
%                   file to be processed.  No extension, in single
%                   quotation marks. If you only want to process one file, enter 0.
%           NumberOfFiles : number of files to be processed (integer).
%           PlotTime: If PlotTime is set to 0, time in seconds from first
%                   measurement will be used for the plots. If it is set to 1, time
%                   in days from 01/01/1900 will be used (matlab time).
%           YAxisRange : The range between the upper and the lower Y axis
%                   limit on the plots.
%
% Outputs:  Ax1 : acceleration values on axis 1 in g's (array of integer)
%           Ax2 : acceleration values on axis 2 in g's(array of integer)
%           Ax3 : acceleration values on axis 3 in g's(array of integer)
%           TotalTime : time of each measurement in matlab time(array of integer)
%           TotalTimeS : time of each measurement in seconds, from time of
%           first measurement(array of integer)
%
%
% This function processes data from Andy Hamilton's accelerometer and plots
% it.
% Script was first created by Andy Hamilton, then modified into a function
% by Francois Cazenave, June 26, 2006
% To run this function, make sure all the files you want to process are in
% your 'current directory'. The files must be in .dat format, raw from the
% accelerometer. File names must not be modified and the time between each
% file must be exactly the same.
% Processing too many files at once may take a long time or go over the
% computer memory limits.
% After processing the raw files from the accelerometer, this function
% plots several diagrams to show acceleration on 3 axis. The time axis can
% be in matlab time or in seconds (see input variable description for details).
% The function also removes bad values (ie values that are out of a defined
% range. This range may be modified in the code below.
% The axis limit may also be change easily in the code.
%
% Example:
% [Ax1, Ax2, Ax3,Time,TimeS] = accel('44AB7180','44AB73D8',5,1,4);
% In the line above, 44AB7180 is the name on the first file and 44AB73D8 the
% name of the second file. 5 files will be processed, including 44AB7180.DAT
% and 44AB73D8.DAT
% X axis on plots will be time in matlab time.
% The Y axis range will be 4 g's.

% Create empty arrays in which the resulting values will be placed
Ax1=[];
Ax2=[];
Ax3=[];
TotalTime=[];

% Find the period length in second
% Files name from the accelerometer are in seconds, in hexadecimal.
% With the name of the first and the second file, we can find the period length
% (time of each cycle for the accelerometer).
Period=[];

if File2~=0; % If user wants to process one file only, next line is skipped 
    Period=hex2dec(File2)-hex2dec(File1);
end

DecFileName=hex2dec(File1);


CurrentFileName=File1;

for i=1:NumberOfFiles;

    % The next block of line was written by Andy Hamilton. What comes out
    % of it is the accceleration along 3 axis and the time of measurements
    % in seconds from 01/01/1970 for the current file
    fid = fopen([CurrentFileName '.DAT'],'r','ieee-be');
    % In case NumberOfFiles is larger than the number of file in the
    % current directory, this exits the loop and goes on with the function.
    if fid==-1;
        break
    end
    StartTime = fread(fid,1,'uint32');
    meas_freq = fread(fid,1,'float32');
    number_of_measurements = fread(fid,1,'uint32');
    calibration_bias1 = fread(fid,1,'float64');
    calibration_offset1 = fread(fid,1,'float64');
    calibration_bias2 = fread(fid,1,'float64');
    calibration_offset2 = fread(fid,1,'float64');
    calibration_bias3 = fread(fid,1,'float64');
    calibration_offset3 = fread(fid,1,'float64');
    Accel1 = fread(fid,number_of_measurements,'uint16');
    Accel2 = fread(fid,number_of_measurements,'uint16');
    Accel3 = fread(fid,number_of_measurements,'uint16');
    Accel1 = calibration_bias1*(Accel1/65535)+calibration_offset1;
    Accel2 = calibration_bias2*(Accel2/65535)+calibration_offset2;
    Accel3 = calibration_bias3*(Accel3/65535)+calibration_offset3;
    Time = 0:length(Accel1)-1;
    Time = Time/meas_freq;

    % These next 3 lines concatenate the acceleration values into 3 big
    % files which are to be returned by the function.
    Ax1=cat(1,Ax1,Accel1);
    Ax2=cat(1,Ax2,Accel2);
    Ax3=cat(1,Ax3,Accel3);

    % Time is given by the accelerometer in seconds from 01/01/1970. This next line
    % calculates the start time of the current file in matlab time.
    StartTime = datenum('01/01/1970') + StartTime/(3600*24);

    % Change the time to matlab time to make plots easier to read
    MatlabTime = Time'./(3600*24) + StartTime;

    % Concatenate the time of each file into one file to be returned by the
    % function.
    TotalTime=cat(1,TotalTime,MatlabTime);

    % Increment the file name, to go on to the next file.
    DecFileName=hex2dec(CurrentFileName)+Period;
    CurrentFileName=dec2hex(DecFileName);

end

% Make a file containing time in seconds from 1st measurement
TotalTimeS=(TotalTime-TotalTime(1))*3600*24;

% Remove bad values. This may be modified by the user depending on the
% files being processed.
Ax1(find(Ax1>10 & Ax1<-10 & Ax1==0))=NaN;
Ax2(find(Ax2>10 & Ax2<-10 & Ax2==0))=NaN;
Ax3(find(Ax1>10 & Ax3<-10 & Ax3==0))=NaN;

% Set the y axis limit for the plots.
% Find which axis in facing down or up and sets the axis limit accordingly.
MeanAx(1)=mean(Ax1);
MeanAx(2)=mean(Ax2);
MeanAx(3)=mean(Ax3);
MeanAx=round(MeanAx);

UpperYLim=YAxisRange/2;
LowerYLim=-UpperYLim;

% If PlotTime is set to 0, time in seconds from first measurement will be
% used for the plots. If it is set to 1, time in days from 01/01/1900 will be used (matlab time).

if PlotTime==0; % use time is seconds from first measurement for x axis
    % Plot the acceleration on 3 subplots
    % Plot acceleration along axis 1
    subplot(3,1,1)
    plot(TotalTimeS,Ax1);
    title(['Unfiltered acceleration along 3 axis from ' datestr(TotalTime(1)) ' to ' datestr(TotalTime(end))]);
    ylim([LowerYLim+MeanAx(1) UpperYLim+MeanAx(1)]);
    ylabel('axis 1 (g)');

    % Plot acceleration along axis 2
    subplot(3,1,2)
    plot(TotalTimeS,Ax2);
    ylim([LowerYLim+MeanAx(2) UpperYLim+MeanAx(2)]);
    ylabel('axis 2 (g)');

    % Plot acceleration along axis 3
    subplot(3,1,3)
    plot(TotalTimeS,Ax3);
    ylim([LowerYLim+MeanAx(3) UpperYLim+MeanAx(3)]);
    xlabel('time (s)')
    ylabel('axis 3 (g)');


elseif PlotTime==1; % use matlab time for x axis
    % Plot the acceleration on 3 subplots
    % Plot acceleration along axis 1
    subplot(3,1,1)
    plot(TotalTime,Ax1);
    title(['Unfiltered acceleration along 3 axis from ' datestr(TotalTime(1)) ' to ' datestr(TotalTime(end))]);
    datetick;
    xlim([TotalTime(1) TotalTime(end)]);
    ylim([LowerYLim+MeanAx(1) UpperYLim+MeanAx(1)]);
    ylabel('axis 1 (g)');

    % Plot acceleration along axis 2
    subplot(3,1,2)
    plot(TotalTime,Ax2);
    datetick;
    xlim([TotalTime(1) TotalTime(end)]);
    ylim([LowerYLim+MeanAx(2) UpperYLim+MeanAx(2)]);
    ylabel('axis 2 (g)');

    % Plot acceleration along axis 3
    subplot(3,1,3)
    plot(TotalTime,Ax3);
    datetick;
    xlim([TotalTime(1) TotalTime(end)]);
    ylim([LowerYLim+MeanAx(3) UpperYLim+MeanAx(3)]);
    xlabel('date');
    ylabel('axis 3 (g)');

end


