function[Gf] = EngrLogger()
    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 = 'H:\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]
    
    
    %% Process IMU data - output floats
        function[data_int16, data_double] = bytes2float(fileName)
        
        % SD Card directory
        filePath = 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
    
    %% Remove offsets
    
    % reference: http://cache.freescale.com/files/sensors/doc/app_note/AN4399.pdf
    g = 981; % [mm/sec^2]
    
    % Turn calibratio on/off
    %calibrateSensor = 1;
    
    % Run offsetCalcAccel('fileName.txt') to obtain Gf values
    % if calibrateSensor == 1
    %     [Gfx, Gfy, Gfz] = offsetCalcAccel(Gf)
    % end
    
    % Calibration values obtained from calibration in six positions
    Gfx = [948.6, -989.1];
    Gfy = [1129.0, -792.2];
    Gfz = [984.7, -981.4];
    
    % Gain value
    Wx = 2 / (Gfx(1)/1000 - Gfx(2)/1000);
    Wy = 2 / (Gfy(1)/1000 - Gfy(2)/1000);
    Wz = 2 / (Gfz(1)/1000 - Gfz(2)/1000);
    
    % Linear offset value
    Vx = -(Gfx(1) + Gfx(2))/(Gfx(1) - Gfx(2))*g;
    Vy = -(Gfy(1) + Gfy(2))/(Gfx(1) - Gfy(2))*g;
    Vz = -(Gfz(1) + Gfz(2))/(Gfx(1) - Gfz(2))*g;
    
    Gx = zeros(length(accelX_double),1);
    Gy = zeros(length(accelX_double),1);
    Gz = zeros(length(accelX_double),1);
    
    % Apply gain and offsets
    Gx(:,1) = Wx.*Gf(:,1) + Vx;
    Gy(:,1) = Wy.*Gf(:,2) + Vy;
    Gz(:,1) = Wz.*Gf(:,3) + Vz;
    
    % Plot the corrected data
    figure;
    hold on
    plot(Gx, 'r');
    plot(Gy, 'g');
    plot(Gz, 'b');
    
    % plot(Gf(:,1), 'Color', r2);
    % plot(Gf(:,2), 'Color', g2);
    % plot(Gf(:,3), 'Color', b2);
    
    xlabel('Samples');
    ylabel('Acceleration (m/sec^2)');
    title('Accelerometer Calibration');
    legend('x-axis', 'y-axis', 'z-axis', 'Location', 'best');
    grid
    hold off
    
    
    %% Filter the data
    
    % De - noise filter
    mu = 0.005;
    Dx = denoiseFilter(Gx, mu);
    Dy = denoiseFilter(Gy, mu);
    Dz = denoiseFilter(Gz, mu);
    
    % Moving average filter
    MAx = movingAverageFilter(Gx, 100);
    MAy = movingAverageFilter(Gy, 100);
    MAz = movingAverageFilter(Gz, 100);
    
    figure;
    hold on
    plot(DNx, 'r');
    plot(DNy, 'g');
    plot(DNz, 'b');
    plot(MVAx, 'Color', r2);
    plot(MVAy, 'Color', g2);
    plot(MVAz, 'Color', b2);
    hold off
    title('Filtered');
    xlabel('Samples');
    ylabel('Acceleration (mm/sec^2)');
    legend('x axis - denoise', 'y axis - denoise', 'z axis - denoise', 'x axis - MVA', 'y axis - MVA', 'z axis - MVA');
    grid on           
     
         
    %% Velocity Estimation
    
    %reference: http://www.mathworks.com/help/matlab/ref/cumtrapz.html
    
    % Numerically integrate acceleration
    velx = [cumtrapz(DNx).*(1/fs), cumtrapz(MVAx).*(1/fs)]; % [mm/sec]
    vely = [cumtrapz(DNy).*(1/fs), cumtrapz(MVAy).*(1/fs)]; % [mm/sec]
    velz = [cumtrapz(DNz).*(1/fs), cumtrapz(MVAz).*(1/fs)]; % [mm/sec]

    % Plot velocity profile and velocity estimates
    figure(3);
    hold on
    plot(velx(:,1),'r');
    plot(vely(:,1),'g');
    plot(velz(:,1),'b');
    
    plot(velx(:,2),'Color', r2);
    plot(vely(:,2),'Color', g2);
    plot(velz(:,2),'Color', b2);
    
    title('Velocity');
    xlabel('Samples');
    ylabel('Velocity (mm/sec)');
    legend('vel estiamte - denoise', 'vel estimate - MVA', 'Location', 'best');
    
    
    %% Plot the data
    
    
    %% Calibration Sequence
    function[ ] = offsetCalcAccel(calibrationData)
        
        % Description: calcualte the calibration constands for accel
        
        % Plot the calibration data
        samples = linspace(0, size(calibrationData, 1), size(calibrationData, 1))';
        
        figure(1)
        plot(samples, calibrationData(:,1), 'r-', samples, calibrationData(:,2), 'g-', samples, calibrationData(:,3), 'b-' );
        xlabel('Samples');
        ylabel('Acceleration (m/sec^2)');
        title('Accelerometer Calibration');
        legend('x-axis', 'y-axis', 'z-axis', 'Location', 'best');
        
        % Identify the data ranges for each calbration position - manually input
        index = [
            50, 1800;
            2345, 3750;
            4260, 5675;
            6228, 7575;
            8061, 9651;
            10030, 11870;
            ];
        
        % Calcualte averages
        data_avg = zeros(6,3);
        
        for i = 1:6
            
            for j = 1:3
                
                data_avg(i,j) = mean(calibrationData(index(i,1):index(i,2),j));
            end
        end
        
        data_avg
        
        % Calcualte offset values
        for i=1:3
            temp  = 0;
            for j = 1:6
                if data_avg(j,i) < 900 && data_avg(j,i) > -900
                    temp = temp + data_avg(j,i);
                end
            end
            offset(1,i) = temp/4;
        end
        
        offset
    end
    
     %% Denoise signal function
    function [ fs ] = denoiseFilter( g, mu )
        
        % Apply the de-noise filter
        fs = denoisetv(g,mu);
        
        function [f] = denoisetv(g,mu)
            I = length(g);
            u = zeros(I,1);
            y = zeros(I,1);
            rho = 10;
            
            eigD = abs(fftn([-1;1],[I 1])).^2;
            
            for k=1:100
                f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
                v = D(f)+(1/rho)*y;
                u = max(abs(v)-1/rho,0).*sign(v);
                y = y - rho*(u-D(f));
            end
        end
        
        function y = D(x)
            y = [diff(x);x(1)-x(end)];
        end
        
        function y = Dt(x)
            y = [x(end)-x(1);-diff(x)];
        end
        
        
    end
    
    %% Apply moving average filter
    function [ mva ] = movingAverageFilter(inputSignal, movingAverageSize)

        % Define moving average filter
        a = 1;
        b = ones(1,movingAverageSize)./movingAverageSize;
        
        % Frequency response of filter
%         [H,F] = freqz(b,a,1024,fs);
%         figure;
%         plot(F,H);
        
        % Apply moving average filter
        mva = filter(b,a,inputSignal);
    end
    
end % parent function

