How to rid HIC calculation of nested for loops

HIC (Head Injury Criteria) is a biomechanics calculation performed on an acceleration time history (https://en.wikipedia.org/wiki/Head_injury_criterion).

Since my FORTRAN days, I've always used nested for loops such as:

function [HICvalue, iT1, iT2, avgAcc] = HIC(RAcc, DeltaT, Win)
    % RAcc is a resultant in "g"s (always +)
    % DeltaT is uniform number of seconds between every point in RAcc
    % Win is a time (in mS) window limit which the calculation is not to exceed
    
    winpts = round((Win/1000)./DeltaT);
    pts = length(RAcc);
    area = cumtrapz(RAcc(:,1)).*DeltaT;
    HICvalue = -1;
    for i=1:winpts
        for j = 1:pts-i
            Temp = (( (1.0/(i.*DeltaT)) * (area(j+i) - area(j)) )^2.5)* i * DeltaT;
            if(Temp > HICvalue)
                HICvalue = Temp;
                iT1 = j;
                iT2 = j+i;
                avgAcc = (1.0/(i.*DeltaT)) * (area(j+i) - area(j));
            end
        end
    end

end

Now that sample frequencies are in the millions, making the acceleration arrays huge, this method has lost its shine [and takes too long]!

Might anyone please help an old dog simplify this [using matrices...], getting rid of the point-wise nested looping?

Answer 1

[Ma Ma]

Thanks for sharing the HIC function.

Below is a vectorized version, The biggest speed-up comes from replacing the power-function with multiplications.

On my machine, the vectorized function fed with 5000 sample points and Win = 15 reduces the computational time by 65%, compared the nested loop function.

function [HICvalue, iT1, iT2,avgAcc] = HIC_vectorized(RAcc, DeltaT, Win)
    % vectorized version of https://www.mathworks.com/matlabcentral/answers/528063-how-to-rid-hic-calculation-of-nested-for-loops
    % Originally by Michael Schlick
    % RAcc is a resultant in "g"s (always +)
    % DeltaT is uniform number of seconds between every point in RAcc
    % Win is a time (in mS) window limit which the calculation is not to exceed
    
    winpts = round((Win/1000)./DeltaT);
    pts = length(RAcc);
    area = cumtrapz(RAcc(:,1)).*DeltaT;
    
    list = pts - (1:winpts);
    num_combinations = sum(list);
    
    I = zeros(num_combinations,1);
    J = zeros(num_combinations,1);
    F = zeros(num_combinations,1);
    r_start = 1;
    for i = 1:winpts
        r_end = r_start-1+pts-i;
        I(r_start:r_end) = i;
        J(r_start:r_end) = 1:(pts-i);
        F(r_start:r_end) = 1/(i*DeltaT);
        r_start = r_end+1;
    end
    
    avgAcc_vec = F.*(area(J+I)-area(J));
    Power2_5 = avgAcc_vec.*avgAcc_vec.*sqrt(avgAcc_vec);
    Temp = Power2_5.*I*DeltaT;
    [HICvalue, index] = max(Temp);
    iT1 = J(index);
    iT2 = J(index) + I(index);
    avgAcc = avgAcc_vec(index);
    
end