function [stepCount, strideTime, stepTime, swingTime, stanceTime, cadence, snLimbStanceTime, dbLimbStanceTime, ...
    strideLen, stepLen, strideReg, gaitSpeed, speed, pressAvg, pressPeak, copSpeed, cop] = processExercise(filename, actId)

	close all;
	
	% variable indexes
	IDX.LACCX = 2;
	IDX.LACCY = 3;
	IDX.LACCZ = 4;
	IDX.LGYRX = 5;
	IDX.LGYRY = 6;
	IDX.LGYRZ = 7;
	IDX.LMAGX = 8;
	IDX.LMAGY = 9;
	IDX.LMAGZ = 10;
	IDX.LP1   = 11;
	IDX.LP2   = 12;
	IDX.LP3   = 13;
	IDX.LP4   = 14;
	IDX.LP5   = 15;
	IDX.RACCX = 16;
	IDX.RACCY = 17;
	IDX.RACCZ = 18;
	IDX.RGYRX = 19;
	IDX.RGYRY = 20;
	IDX.RGYRZ = 21;
	IDX.RMAGX = 22;
	IDX.RMAGY = 23;
	IDX.RMAGZ = 24;
	IDX.RP1   = 25;
	IDX.RP2   = 26;
	IDX.RP3   = 27;
	IDX.RP4   = 28;
	IDX.RP5   = 29;
    
    I.ACCX = 1; I.ACCY = 2; I.ACCZ = 3;
    I.GYRX = 4; I.GYRY = 5; I.GYRZ = 6;
    I.MAGX = 7; I.MAGY = 8; I.MAGZ = 9;
    I.P1   = 10; I.P2  = 11; I.P3  = 12; I.P4 = 13; I.P5 = 14;
	
	COL.BLUE  = [0      0.4470 0.7410];
	COL.RED   = [0.8509 0.3294 0.1019];
    COL.GREEN = [0.4660 0.6740 0.1880];
    COL.RED2  = [0.6350 0.0780 0.1840];

	data = csvread(filename, 1, 0);
    
    ts = data(:,1) / 1000;
    
	wl = sum(data(:, IDX.LP1:IDX.LP5), 2);
	wr = sum(data(:, IDX.RP1:IDX.RP5), 2);
	
    [lhs, lto] = computeHS_TO_new (data(:,  2:15), COL.BLUE, 'L');
	[rhs, rto] = computeHS_TO_new (data(:, 16:29), COL.RED,  'R');
    
	stridesTime = diff(rhs);
    stridesReg = [];
    for h = 1:size(stridesTime)
        s = stridesTime(max(1, h-3):h);
        mu = sum(s) / size(s, 1);
        sigma = sqrt(sum(abs(stridesTime(h) - mu)) / size(s, 1));
        stridesReg = [stridesReg; sigma / mu];
    end
	
	[lst, rst] = computeStepTimes (lhs, rhs);
	stepsTime = [lst.time; rst.time];
    
    lswt = computeSwingTimes (lto, lhs);
    rswt = computeSwingTimes (rto, rhs);
    swingsTime = [lswt.time; rswt.time];
	
    lstt = computeStanceTimes (lto, lhs);
    rstt = computeStanceTimes (rto, rhs);
    stancesTime = [lstt; rstt];
    
    cadences = 60./stridesTime;
    
    % --- SPEED-RELATED STATS
    [lThetaInit, rThetaInit] = calibrate (2);  % use the first 2 seconds for calibration data
    [lStepLen, lStrideLen, lGaitSpd, lSpd] = computeDistance (data(:,  IDX.LACCX:IDX.LP5), lThetaInit, lswt, 'Left');
    [rStepLen, rStrideLen, rGaitSpd, rSpd] = computeDistance (data(:,  IDX.RACCX:IDX.RP5), rThetaInit, rswt, 'Right');
    stepsLen = [lStepLen; rStepLen];
    stridesLen = [lStrideLen; rStrideLen];
    figure; histogram (stridesLen, [0.2:0.2:2.0]);
    speeds = [lSpd; rSpd];
    gaitsSpeed = [lGaitSpd; rGaitSpd];
    % -------------------------
    
    ldbst = computeStanceTimes (lto, rhs);
    rdbst = computeStanceTimes (rto, lhs);
    dblst = [ldbst; rdbst];
    
    [pavg, peaks] = computePressAvgAndPeaks (wl, wr);
    
    lcop = computeCoP(data(:, IDX.LP1:IDX.LP5), 'l');
    rcop = computeCoP(data(:, IDX.RP1:IDX.RP5), 'r');
    
    lcopSpeed = computeCoPSpeed (ts, lcop);
    rcopSpeed = computeCoPSpeed (ts, rcop);
    copSpeed = [lcopSpeed(~isnan(lcopSpeed)); rcopSpeed(~isnan(rcopSpeed))]; % remove NaN values (no pressure)
    
    copaggr = computeCoPAggr (lcop, rcop);
    
    % output values assignments
	stepCount = size (stepsTime, 1); 
    strideTime = getAggrStat(stridesTime);
	stepTime = getAggrStat(stepsTime);
    swingTime = getAggrStat(swingsTime);
    stanceTime = getAggrStat(stancesTime);
    cadence = getAggrStat(cadences);
    snLimbStanceTime = swingTime;
    dbLimbStanceTime = getAggrStat(dblst);
    strideLen = getAggrStat(stridesLen);
    stepLen = getAggrStat(stepsLen);
    strideReg = getAggrStat(stridesReg);
    gaitSpeed = getAggrStat(gaitsSpeed);
    speed = getAggrStat(speeds);
    pressAvg = getAggrStat(pavg);
    pressPeak = getAggrStat(peaks);
    copSpeed = getAggrStat(copSpeed);
    cop = getAggrStat(copaggr);
    
    res = [strideTime; stepTime; cadence; strideLen; stepLen; gaitSpeed; speed];
    confInt = res(:, 3) - res(:,2);
    res = [res(:, 1) confInt];
    res = reshape(res',1,[]);
    res = [stepCount res];
    
    disp (res);
    
    % ---- support functions ----
    
    function [lThetaInit, rThetaInit] = calibrate (startInterval)
        len = startInterval / 0.02;   % samples
        
        % check if the first <len> samples have low variance i.e. can be used for calibration
        gyrnorm = data(1:len, IDX.LGYRX:IDX.LGYRZ) ./ (max(data(:, IDX.LGYRX:IDX.LGYRZ)) - min(data(:, IDX.LGYRX:IDX.LGYRZ)));
        gyrmag = sqrt(sum((gyrnorm).^2, 2));
        gyrvar = var(gyrmag);
        
        if gyrvar < 0.0001
            % remove the initial offset of the gyro samples
            data(:, IDX.LGYRX:IDX.LGYRZ) = data(:, IDX.LGYRX:IDX.LGYRZ) - mean(data(1:len, IDX.LGYRX:IDX.LGYRZ));
            data(:, IDX.RGYRX:IDX.RGYRZ) = data(:, IDX.RGYRX:IDX.RGYRZ) - mean(data(1:len, IDX.RGYRX:IDX.RGYRZ));

            % rimodula acc x/y/z su modulo gravita' (a riposo)
            lgOff = mean(data(1:len,IDX.LACCX:IDX.LACCZ) - (data(1:len,IDX.LACCX:IDX.LACCZ) ./ sqrt(sum((data(1:len,IDX.LACCX:IDX.LACCZ)/1000).^2, 2))));
            rgOff = mean(data(1:len,IDX.RACCX:IDX.RACCZ) - (data(1:len,IDX.RACCX:IDX.RACCZ) ./ sqrt(sum((data(1:len,IDX.RACCX:IDX.RACCZ)/1000).^2, 2))));
            %data(:, IDX.LACCX:IDX.LACCZ) = data(:, IDX.LACCX:IDX.LACCZ) - lgOff;
            %data(:, IDX.RACCX:IDX.RACCZ) = data(:, IDX.RACCX:IDX.RACCZ) - rgOff;

            % compute initial resting angles of the shoes
            lThetaInit = mean(atan(-data(1:len,IDX.LACCY) ./ data(1:len,IDX.LACCZ)));
            rThetaInit = mean(atan(-data(1:len,IDX.RACCY) ./ data(1:len,IDX.RACCZ)));
        else
            disp ('[WARN] Cannot calibrate using the starting sequence');
            lThetaInit = 0;
            rThetaInit = 0;
        end
    end
    
    function [stepLen, strideLen, gaitSp, speed] = computeDistance (d, thetaInit, swings, side)

        startIdx = swings.startIdx;
        endIdx = swings.endIdx;
        
        figure;
        
        % filtering gyrX
        filtCutOff = 0.001;
        [b, a] = butter(1, (2*filtCutOff) / 50, 'high');
        d(:, I.GYRX) = filtfilt(b, a, d(:, I.GYRX));
        
        filtCutOff = 5;
        [b, a] = butter(1, (2*filtCutOff) / 50, 'low');
        d(:, I.GYRX) = filtfilt(b, a, d(:, I.GYRX));
        
        theta = cumtrapz(ts, d(:, I.GYRX) ./ 1000.0) + thetaInit;   % mrad/s to rad/s
        plot(ts, rad2deg(theta), 'DisplayName', 'theta(deg)');
        hold on;
        
        hpTh = plot(ts(1), 0, 'DisplayName', 'Corrected step theta(deg)');
        hpStart = plot(ts(1), 0, '--', 'Color', COL.GREEN, 'DisplayName', 'Step Start');
        hpEnd = plot(ts(1), 0, '--', 'Color', COL.RED2, 'DisplayName', 'Step End');
        hpSpd = plot(ts(1), 0, 'DisplayName', 'Step Speed');
        hpAcc = plot(ts(1), 0, 'DisplayName', 'Step Accel');
        hpDst = plot(ts(1), 0, 'DisplayName', 'Step Dist');
        
        magAccGyr = sqrt(sum(d(:,I.GYRX:I.GYRZ).^2, 2)) .* sqrt(sum(data(:,I.ACCX:I.ACCZ).^2, 2));
        
        stepLen = [];
        strideLen = [];
        speed = [];
        gaitSp = [];
        
        for i = 1:min(size(startIdx, 1), size(endIdx, 1))

            % zero-movement instant finder method
            backwardLookingWindow = max(startIdx(i) - 25, 1);
            [a, b] = min(magAccGyr(backwardLookingWindow:startIdx(i)));
            zuptBackIdx = backwardLookingWindow - 1 + b;
            startTheta = atan(-d(zuptBackIdx,I.ACCY) ./ d(zuptBackIdx,I.ACCZ));

            % find forward zero-movement instant
            forwardLookingWindow = min(endIdx(i) + 25, size(ts, 1));
            [a, b] = min(magAccGyr(endIdx(i):forwardLookingWindow));
            zuptFwdIdx = endIdx(i) - 1 + b;
            endTheta = atan(-d(zuptFwdIdx,I.ACCY) ./ d(zuptFwdIdx,I.ACCZ));
            
            % select and integrate stride motion window
            idx = zuptBackIdx:zuptFwdIdx;
            th = cumtrapz(ts(idx), d(idx, I.GYRX) ./ 1000.0) + startTheta;
            
            % correct theta drift (align start and end with zero)
            drift = (th(end) - th(1) - (endTheta - startTheta)) / size(idx, 2);
            th = th - [0; ((1:(size(idx,2)-1))' * drift)];
            
            set(hpTh,'XData',[get(hpTh,'XData') NaN ts(idx)'], 'YData', [get(hpTh,'YData') NaN rad2deg(th)']);
            set(hpStart,'XData', [get(hpStart,'XData') NaN [ts(startIdx(i))' ts(startIdx(i))']], ...
                        'YData', [get(hpStart,'YData') NaN [rad2deg(min(theta(startIdx(i))-0.1, -0.87)) rad2deg(max(theta(startIdx(i))+0.1, 0.87))]]);
            set(hpEnd,'XData', [get(hpEnd,'XData') NaN [ts(endIdx(i))' ts(endIdx(i))']], ...
                      'YData', [get(hpEnd,'YData') NaN [rad2deg(min(theta(endIdx(i))-0.1, -0.87)) rad2deg(max(theta(endIdx(i))+0.1, 0.87))]]);
            
            fwacc = -1 * (((-d(idx, I.ACCY) * 0.00981) .* cos(th)) - ((d(idx, I.ACCZ) * 0.00981) .* sin(th)));
            
            fwspd = cumtrapz(ts(idx), fwacc);
            
            % correct speed drift (align start and end with zero)
            drift = (fwspd(end) - fwspd(1)) / size(idx, 2);
            fwspd = fwspd - [0; ((1:(size(idx,2) - 1))' * drift)];
            
            % compute forward distance
            fwdst = cumtrapz(ts(idx), fwspd);
            
            strideLen = [strideLen; fwdst(end)];
            
            [~, stepIdxStart] = min(sqrt(sum(d(idx,I.ACCX:I.ACCZ).^2, 2)));
            stepLen = [stepLen; fwdst(end) - fwdst(stepIdxStart)];

            strideSpd = strideLen / (ts(idx(end)) - ts(idx(1)));
            speed = [speed; strideSpd];
            
            if strideSpd < 0.2
                gaitSp = [gaitSp; 1];
            elseif strideSpd < 0.6
                gaitSp = [gaitSp; 2];
            else
                gaitSp = [gaitSp; 3];
            end
            
            set(hpAcc,'XData',[get(hpAcc,'XData') NaN ts(idx)'], 'YData', [get(hpAcc,'YData') NaN fwacc']);
            set(hpSpd,'XData',[get(hpSpd,'XData') NaN ts(idx)'], 'YData', [get(hpSpd,'YData') NaN fwspd']);
            set(hpDst,'XData',[get(hpDst,'XData') NaN ts(idx)'], 'YData', [get(hpDst,'YData') NaN fwdst']);
        end
        
        title([side sprintf('Foot [StrideLen Avg: %.2f Median: %.2f Min: %.2f Max: %.2f],    AvgStepLen: %.2fm', ...
            mean(strideLen), median(strideLen), min(strideLen), max(strideLen), mean(stepLen))]);
        
    end
    
    function [hs, to] = computeHS_TO_new (d, color, side)
		
        dmin = min(d);
        dmax = max(d);
        
        dn = d ./ (dmax - dmin);
        dn(isnan(dn)) = 0;       % NaN results when min and max equals 0
        
        wt = sum(dn(:, 10:14), 2);
        wt = wt ./ max(wt);
        
        wtoe = sum(dn(:, IDX.RP3:IDX.RP5), 2);
        wheel = sum(dn(:, IDX.RP1:IDX.RP2), 2);
        
        [pp, pts, ww, pr] = findpeaks (-dn(:, I.GYRX), ts, 'MinPeakDistance', 0.2, 'MinPeakProminence', 0.05, 'MinPeakHeight', 0.05);
        
        mmin = movmin(-dn(:, I.GYRX), [10 0]) < -0.05;
        hs = ts(mmin & ismember(ts, pts));
        to = ts(~mmin & ismember(ts, pts));
        
        figure;
        plot(ts, wt, 'DisplayName', [side 'Weight']); hold on;
        plot(ts, -dn(:, I.GYRX), 'DisplayName', ['-' side 'GyrX']);
        plot(hs, wt(ismember(ts, hs)), 'x', 'MarkerSize', 10, 'MarkerEdgeColor', color, 'DisplayName', [side 'HeelStrike']);
        plot(to, wt(ismember(ts, to)), 'o', 'MarkerSize', 10, 'MarkerEdgeColor', color, 'DisplayName', [side 'ToeOff']);
    end
    
    % compute Heel Strike and Toe Off events
	function [hs, to] = computeHS_TO (tsp, w, color, side)
        
		% heel strike
		dw = diff(w + 1)./(w(1:end-1) + 1);
		[pks, tsHS, ww, pp] = findpeaks ([0; dw], tsp, 'MinPeakDistance', 0.3, 'MinPeakProminence', 10, 'MinPeakHeight', 16);
		plot(tsHS, w(ismember(tsp, tsHS)), 'x', 'MarkerSize', 10, 'MarkerEdgeColor', color, 'DisplayName', [side 'HeelStrike']);
		hs = tsHS;
		
		% toe off
		fw = flip (w);
		dw = diff(fw + 1)./(fw(1:end-1) + 1);
		dw = [0; dw].*fw;
		[pks, idx, ww, pp] = findpeaks (dw, 'MinPeakDistance', 15, 'MinPeakProminence', 50, 'MinPeakHeight', 10000);
		idx = size(tsp, 1) - idx + 2;
		plot(tsp(idx), w(idx), 'o', 'MarkerSize', 10, 'MarkerEdgeColor', color, 'DisplayName', [side 'ToeOff']);
		to = flip(tsp (idx));
	end

	% compute left and right step times, intervals [RHS <---> LHS] and [LHS <---> RHS]
	function [lst, rst] = computeStepTimes (lhs, rhs)
		hs       = [lhs; rhs];
		hsisleft = [ones(size(lhs)); zeros(size(rhs))];
		[fs, so] = sort (hs);
		ss = hsisleft(so);
		
		i = 2;
		lst.startIdx = [];
        lst.endIdx = [];
        lst.time = [];
		rst.startIdx = [];
        rst.endIdx = [];
        rst.time = [];
		while i <= size (ss, 1)
			if ss(i) ~= ss(i-1)
				if (ss(i) == 1) % is left step
					lst.time     = [lst.time;     fs(i) - fs(i-1)];
                    lst.startIdx = [lst.startIdx; find(ts == fs(i-1))];
                    lst.endIdx   = [lst.endIdx;   find(ts == fs(i))];
				elseif (ss(i) == 0) % is right step
					rst.time     = [rst.time;     fs(i) - fs(i-1)];
                    rst.startIdx = [rst.startIdx; find(ts == fs(i-1))];
                    rst.endIdx   = [rst.endIdx;   find(ts == fs(i))];
				end
			end
			i = i + 1;
		end
    end

    % compute swing times, intervals [TO <----> HS]
	function swt = computeSwingTimes (to, hs)
		seq       = [to; hs];
		seqisTO   = [ones(size(to)); zeros(size(hs))];
		[fs, so]  = sort (seq);
		ss = seqisTO(so);
        i = 2;
        
        swt.startIdx = [];
        swt.endIdx = [];
        swt.time = [];
        while i <= size (ss, 1)
            if (ss(i) == 0 && ss(i-1) == 1)
                swt.time     = [swt.time;     fs(i) - fs(i-1)];
                swt.startIdx = [swt.startIdx; find(ts == fs(i-1))];
                swt.endIdx   = [swt.endIdx;   find(ts == fs(i))];
            end
            i = i + 1;
        end
    end

    % compute stance times, intervals [HS <----> TO]
	function stt = computeStanceTimes (to, hs)
		seq       = [to; hs];
		seqisTO   = [ones(size(to)); zeros(size(hs))];
		[fs, so]  = sort (seq);
		ss = seqisTO(so);
        i = 2;
        stt = [];
        while i <= size (ss, 1)
            if (ss(i) == 1 && ss(i-1) == 0)
                stt = [stt; fs(i) - fs(i-1)];
            end
            i = i + 1;
        end
    end

    % compute pressure average during pressure
	function [pavgs, peaks] = computePressAvgAndPeaks (wl, wr)
        
        PRESS_THRESH  = 100;     % total pressure threshold
        MIN_PRESS_DUR = 1;       % minimum pressure interval duration
        
        wt = wl + wr;
        wti = wt > PRESS_THRESH; % thresholding
        
        [start, len, count] = extractSequences (wti);
        
        start = start (len > MIN_PRESS_DUR / 0.02);
        len = len(len > MIN_PRESS_DUR / 0.02);    % select only intervals of at least 1 sec
        
        pavgs = [];
        peaks = [];
        for k = 1:size(start,2)
            pavgs = [pavgs; mean(wt(start(k):start(k)+len(k)-1))];
            peaks = [peaks; max(wt(start(k):start(k)+len(k)-1))];
        end
    end

    % compute center of pressure
	function [cop] = computeCoP (p, side)
        SPOS = [ 0.78 -9.42;   % expressed in cm
                -3.11 -9.25;
                -2.86  9.36;
                 0.11  7.07;
                 3.30  9.27];
        SPOS = SPOS / 100.0;   % converts to meters
        
        if side == 'r'
            SPOS(:,1) = -SPOS(:,1);
        end
        
        copx = [];
        copy = [];
        for i = 1:size(p,1)
            copx = [copx; dot(p(i, 1:5), SPOS(:, 1)) / sum(p(i, 1:5))];
            copy = [copy; dot(p(i, 1:5), SPOS(:, 2)) / sum(p(i, 1:5))];
        end
        cop = [copx copy];
        
    end

    % compute center of pressure speed
	function copSpeed = computeCoPSpeed (ts, cop)
        
        [start, len, count] = extractSequences (~any(isnan(cop), 2));
        
        copSpeed = NaN([size(ts, 1) 1]);
        for i = 1:count
            if (len(i) >= 2)
                for j = start(i)+1:start(i)+len(i)-1 
                    tdiff = ts(j) - ts(j-1);
                    sdiff = sqrt((cop(j, 1) - cop(j-1, 1))^2 + (cop(j, 2) - cop(j-1, 2))^2);
                    copSpeed(j) = sdiff/tdiff; % expressed as cm/s
                end
            end
        end
    end

    % compute cop aggregated stats
	function copaggr = computeCoPAggr (lcop, rcop)
        % calcolo della radice quadrata di (Xd-Xs)^2+(Yd-Ys)^2, assumendo che a piede alzato
        % si assegni un valore fisso nominale di X ed Y fuori dai valori assumibili
        % a piede a terra (es. Y=20cm X=10cm).
        
        % put 10cm where X is NaN
        lcop(isnan(lcop(:,1)), 1) = 0.10;
        rcop(isnan(rcop(:,1)), 1) = 0.10;
        
        % put 20cm where Y is NaN
        lcop(isnan(lcop(:,2)), 2) = 0.20;
        rcop(isnan(rcop(:,2)), 2) = 0.20;
        
        copaggr = sqrt((rcop(:,1) - lcop(:,1)).^2 + (rcop(:,2) - lcop(:,2)).^2);
    end

    function [start, len, count] = extractSequences (idx)
        n = length(idx);
        start = [];            % where each 1-string starts
        len = [];              % length of each 1-string
        count = 0;             % count of 1-strings
        inOnes = false;
        for k = 1:n
            if idx(k) == 0             % not in 1-string
                inOnes = false;
            elseif ~inOnes             % start of new 1-string
                inOnes = true;
                count = count + 1;
                start = [start k];
                len = [len 1];
            else                       % still in 1-string
                len(count) = len(count) + 1;
            end
        end
    end
    
    % compute mean and confidence intervals of an array of data
    function aggr = getAggrStat (v)
        m = mean(v);
        
        % z statistic (we know the stdev of the population)
        %ci = 1.96 * (std(v) / sqrt(size(v, 1)));
        
        % t statistic (we know the stdev of the sample)
        ci = 2.1 * (std(v) / sqrt(size(v, 1)));
        
        aggr = [m max(0, m-ci) m+ci];
    end

    function plotData (d)
        
        t = d(:, 1);
        dmin = min(d);
        dmax = max(d);
        
        dn = d ./ (dmax - dmin);
        
        wsx = sum(dn(:, IDX.LP1:IDX.LP5), 2);
        wsx = wsx ./ max (wsx);
        
        wdx = sum(dn(:, IDX.RP1:IDX.RP5), 2);
        wdxtoe = sum(dn(:, IDX.RP3:IDX.RP5), 2);
        wdxheel = sum(dn(:, IDX.RP1:IDX.RP2), 2);
        wdx = wdx ./ max (wdx);
        
        [pp, pts, ww, pr] = findpeaks (-dn(:, IDX.RGYRX), ts, 'MinPeakDistance', 0.2, 'MinPeakProminence', 0.1, 'MinPeakHeight', 0.1);

        figure;
        plot (ts, wdx); hold on;
        %plot ([pts pts]', [-ones(1, size(pts,1),1); ones(1, size(pts,1),1)], 'x', 'MarkerSize', 10, 'MarkerEdgeColor', 'r', 'DisplayName', ['HeelStrike']);
        %plot (ts, wdxtoe); hold on;
        %plot (ts, wdx); hold on;
        %plot (ts, dn(:, IDX.RGYRX));
        %plot (ts, dn(:, IDX.RACCZ));
        
    end
end

