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?
0 个评论
采纳的回答
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 个评论
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!