function[Gf] = TableTiltError()
    close all
    clear fig
    
    % Description:
    %
    % Step 1: import byte data from text file
    % Step 2: convert bytes to floats
    % Step 3: apply correction offsets and scaling factors
    % Step 4: filter the data
    % Step 5: numerically integrate to obtain velocity estimate
    % Step 6: compute error based on desired velocity profile
    
    %--------------------------------------------------------------------------
    %--------------------------------------------------------------------------
    
    % Custom colors for plotting
    r2 = [0.6902 0.0275 0.2706];
    g2 = [59/255 69/255 150/255];
    b2 = [7/255 176/255 159/255];

    %Name of file to import
    fileName = 'C:\Users\nraymond\Desktop\sd temp\01-Jun-2011-22-50-39';
        
    [accelX_int16, accelX_double] = bytes2float(strcat(fileName,'_accelX.txt'));
    [accelY_int16, accelY_double] = bytes2float(strcat(fileName,'_accelY.txt'));
    [accelZ_int16, accelZ_double] = bytes2float(strcat(fileName,'_accelZ.txt'));
    
    % Raw sensor output
    Gf = [accelX_int16, accelY_int16, accelZ_int16];
    
    % Sample rate 
    fs = 100; % [Hz]
    
    
    %% Import IMU data - output floats
        function[data_int16, data_double] = bytes2float(fileName)
        
        % SD Card directory
        filePath = fileName; %strcat('H:\', fileName);
        fileID = fopen(filePath);
        
        % Directory on computer - DOWNLOADS
        %fileID = fopen('C:\Users\nraymond\Downloads\01-Jun-2011-01-00-30_accelZ.txt');
        
        % Read the file, assuming signed 8 bit integers
        A = fread(fileID);
        
        % Close the file
        fclose('all');
        
        numSamples = length(A)/2;
        
        %Data_byte = zeros(numSamples,2);
        data_double = zeros(numSamples,1);
        
        index = 1;
        for i=1:2:numSamples*2
            MyByteArray = [A(i) A(i+1)];
            MyByteArray = uint8(MyByteArray);
            MySingle = typecast(MyByteArray,'int16');
            accelDouble = double(MySingle)/100.0; % [m/sec^2]
            
            data_int16(index,1) = MySingle;
            data_double(index,1) = accelDouble;
            
            index = index + 1;
        end
        
        end
    
        %% Identify the data range of samples for processing data

        readRawData = [
                        1, 947;
                        1279, 3179;
                        3793, 5644;
                        6020, 7874;
                        8296, 10090;
                        10530, 12360;
                        12750, 14620;
                        15050, 16850;
                        17240, 19120;
                        19490, 21370;
                        21750, 23600;
                        ];
        
%% Calcualte averages
    calibrationData = Gf;
    
    data_avg = zeros(length(readRawData),3);
    data_rms = zeros(length(readRawData),3);
    data_std = zeros(length(readRawData),3);
    data_CI = zeros(length(readRawData),3);
    distance = zeros(length(readRawData),1);
    
    alpha = 0.05;

    for p = 1:length(readRawData)

        for q = 1:3
            x = calibrationData(readRawData(p,1):readRawData(p,2),q);
            data_avg(p,q) = mean(double(x));
            data_rms(p,q) = rms(double(x));
            data_std(p,q) = std(double(x));
            data_CI(p,q) = confidence_T(alpha,  double(x), 2); % use 2-tail t-test
            
            % calcualte pitch and roll angles from accel measurements
            % NEED TO REMOVE OFFSETS FIRST
%             if q==2
%                 angle_x(p,1) = 90 - atan2d(9807, data_avg(p,q)); % [deg]
%             end
%             
%             if q ==3
%                 angle_y(p,1) = 90 - atan2d(9807, data_avg(p,q)); % [deg]
%             end
            
        end
        
        distance(p,1) = p*50; % [mm]
        
    end
    
 
    
    %% Plot the corrected data
    figure(1);
    hold on
    plot(Gf(:,1), 'r');
    plot(Gf(:,2), 'g');
    plot(Gf(:,3), 'b');
   
    xlabel('Samples');
    ylabel('Acceleration (m/sec^2)');
    title('Linear Translation Stage - Flatness Test');
    legend('x-axis', 'y-axis', 'z-axis', 'Location', 'best');
    grid
    hold off
    
    %% Plot the average values
    figure(2);
    hold on
    plot(distance, data_avg(:,1), 'ro-');
    plot(distance, data_avg(:,2), 'go-');
    plot(distance, data_avg(:,3), 'bo-');
   
    xlabel('Linear Table Distance (mm)');
    ylabel('Acceleration (mm/sec^2)');
    title('Average Sensor Output Values (int16)');
    legend('x-axis', 'y-axis', 'z-axis', 'Location', 'best');
    grid
    hold off
    
    %% Plot the RMS values
     figure(3);
    hold on
    plot(distance, data_rms(:,1), 'ro-');
    plot(distance, data_rms(:,2), 'go-');
    plot(distance, data_rms(:,3), 'bo-');
   
    xlabel('Linear Table Distance (mm)');
    ylabel('Acceleration (mm/sec^2)');
    title('RMS of Sensor Ouput Values');
    legend('x-axis', 'y-axis', 'z-axis', 'Location', 'best');
    grid
    hold off
    
    %% Plot the STD values
    figure(4);
    hold on
    plot(distance, data_std(:,1), 'ro-');
    plot(distance, data_std(:,2), 'go-');
    plot(distance, data_std(:,3), 'bo-');
    
    xlabel('Linear Table Distance (mm)');
    ylabel('Acceleration (mm/sec^2)');
    title('Standard Deviation of Sensor Output Values');
    legend('x-axis', 'y-axis', 'z-axis', 'Location', 'best');
    grid
    hold off
    
    
    
    %% Plot the CI values
    figure(5);
    hold on
    errorbar(distance, data_avg(:,1), data_CI(:,1), 'ro-');
    errorbar(distance, data_avg(:,2), data_CI(:,2), 'go-');
    errorbar(distance, data_avg(:,3), data_CI(:,3), 'bo-');
    
    xlabel('Linear Table Distance (mm)');
    ylabel('Acceleration (mm/sec^2)');
    title('95% CI Values');
    legend('x-axis', 'y-axis', 'z-axis', 'Location', 'best');
    grid
    hold off
    %% Function Definitions
    
    function[CI] = confidence_T(alpha, input_data, tail)
        
        % reference: t_inv - https://www.mathworks.com/matlabcentral/answers/260564-how-to-implement-tinv-manually-without-statistics-and-machine-learning-toolbox
        % Inputs:
        %
        % alpha
        % standard_dev: standard deviation of sample
        % size: number of data points in sample
        % tail:
        %    if tail = 1; 1-tailed t-distribution
        %    if tail = 2; 2-tailed t-distribution
        
        % Local functions
        tdist2T = @(t,v) (1-betainc(v/(v+t^2),v/2,0.5));        % 2-tailed t-distribution
        tdist1T = @(t,v) 1-(1-tdist2T(t,v))/2;                  % 1-tailed t-distribution
        
        % Calculations  
        v = length(input_data)-1;                               % degrees of freedom
        SEM = std(input_data)/sqrt(length(input_data));         % Standard Error
        
        if tail == 1
            t_inv = @(alpha,v) fzero(@(tval) (max(alpha,(1-alpha)) - tdist1T(tval,v)), 5);  % T-Statistic Given Probability ‘alpha’ & Degrees-Of-Freedom ‘v’
            ts = t_inv(alpha,v);                                % T-Score
        end
        
        if tail == 2
            t_inv = @(alpha,v) fzero(@(tval) (max(alpha,(1-alpha)) - tdist2T(tval,v)), 5);  % T-Statistic Given Probability ‘alpha’ & Degrees-Of-Freedom ‘v’
            ts = t_inv(alpha,v);                                % T-Score
        end
        
        CI = ts*SEM;                         % Confidence Intervals
    end
    
end