How to rid HIC calculation of nested for loops

38 次查看(过去 30 天)
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?

采纳的回答

Ma Ma
Ma Ma 2020-12-10
编辑:Ma Ma 2020-12-11
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

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by