2. Compute HIC15 from Acceleration Time History

What HIC15 Actually Computes

The Head Injury Criterion is defined over a resultant translational head acceleration a(t) expressed in units of g. For a window of duration t₂ − t₁, the criterion is

HIC = ( t₂ − t₁ ) · [ (1/(t₂ − t₁)) · ∫_t1^t₂ a(t) dt ]^2.5

with the constraint t₂ − t₁ ≤ 15 ms for HIC15. The quantity inside the square brackets is the mean resultant acceleration over the window, so HIC15 is the maximum over all admissible windows of duration × mean^2.5. Two consequences drive the implementation: the metric is dominated by the mean rather than the instantaneous peak, and the search must consider every window length from a single sample up to 15 ms, not just 15 ms windows.

Because the integral is evaluated on discrete samples, the practical form uses trapezoidal integration over the sample interval Δt. For a window spanning samples i to j:

HIC = ( t_j − t_i ) · [ ( Δt / ( t_j − t_i ) ) · ( a_i/2 + Σ_k=i+1^j−1 a_k + a_j/2 ) ]^2.5

Note that the leading duration term and the 1/duration inside the bracket cancel, leaving HIC = Δt^2.5 · ( Σ_trapezoid )^2.5 / ( t_j − t_i )^1.5. This algebraic simplification is worth using: it removes one division per candidate window and makes the sliding-window search cheap.

Reference Dataset and Expected Value

The validation case below uses a synthetic but fully specified acceleration history so the result is reproducible without external files. It is a half-sine pulse of 80 g peak and 20 ms base duration, sampled at 10 kHz (Δt = 0.1 ms), with a small secondary lobe added to make the window search non-trivial.

import numpy as np

def build_reference_curve(dt=1e-4, duration=0.060):
    t = np.arange(0.0, duration + dt/2, dt)
    a = np.zeros_like(t)
    # primary half-sine: 80 g peak, 20 ms base
    m1 = (t >= 0.005) & (t <= 0.025)
    a[m1] = 80.0 * np.sin(np.pi * (t[m1] - 0.005) / 0.020)
    # secondary lobe: 35 g peak, 8 ms base, starting at 32 ms
    m2 = (t >= 0.032) & (t <= 0.040)
    a[m2] = 35.0 * np.sin(np.pi * (t[m2] - 0.032) / 0.008)
    return t, a

For this curve the expected HIC15 is approximately 742 (dimensionless). That figure is a computed property of the stated analytic pulse, not a measured test result; it is the reference value the implementation must reproduce to within a small tolerance.

Sliding-Window Implementation

The search enumerates every start sample i, advances the end sample j while the window duration stays within 15 ms, and tracks the maximum. A prefix-sum array of the trapezoidal integrand makes each window evaluation O(1).

def hic15(t, a, window=0.015):
    dt = t[1] - t[0]
    n = len(a)
    # trapezoidal cumulative integral of a(t)
    integ = np.empty(n)
    integ[0] = 0.0
    integ[1:] = np.cumsum(0.5 * (a[1:] + a[:-1]) * dt)
    best_hic, best_i, best_j = 0.0, 0, 0
    j = 0
    for i in range(n):
        if j < i:
            j = i
        while j + 1 < n and (t[j + 1] - t[i]) <= window:
            j += 1
        for k in range(i, j + 1):
            dur = t[k] - t[i]
            if dur <= 0.0:
                continue
            area = integ[k] - integ[i]          # g·s
            mean_a = area / dur                 # g
            hic = dur * mean_a ** 2.5
            if hic > best_hic:
                best_hic, best_i, best_j = hic, i, k
    return best_hic, t[best_i], t[best_j]

Two details matter for correctness. First, the inner loop over k is required because the maximum is not always at the longest admissible window — a shorter, higher-mean window can win. Second, the window bound uses <= window so that a window of exactly 15 ms is admissible; using a strict inequality silently biases the result low.

Window Geometry and Peak Identification

The SVG below shows the reference curve with the winning HIC15 window overlaid, together with the mean-acceleration profile that the search is implicitly maximizing. The key visual point is that the optimal window is shorter than 15 ms and is positioned to capture the steepest part of the primary pulse, not its absolute peak.

The winning window starts at the pulse onset and ends before the peak, because extending the window past the peak adds low-acceleration samples that reduce the mean faster than the duration term can compensate. This is the characteristic HIC behavior that distinguishes it from a simple peak metric, and it is why the search must be exhaustive over window lengths.

Validation and Failure Symptoms

Run the implementation against the reference curve and compare to the expected value:

t, a = build_reference_curve()
hic, t1, t2 = hic15(t, a)
print(f"HIC15 = {hic:.1f}  window = [{t1*1e3:.2f}, {t2*1e3:.2f}] ms")
# expected: HIC15 ≈ 742, window ≈ [5.00, 17.50] ms
assert abs(hic - 742.0) < 5.0, "HIC15 outside tolerance"

Useful diagnostics when the result deviates:

• Result too low by a factor near 1.5–2: the window bound is strict (< instead of <=) or the search only evaluates full 15 ms windows. Check that the inner loop over k exists.
• Result too high and window pinned at 15 ms: the acceleration is still in m/s² rather than g. Divide by 9.80665 before the search; the 2.5 exponent amplifies unit errors severely.
• Result sensitive to sampling rate: the trapezoidal integral is converging, but if Δt is coarse the discrete mean underestimates the continuous one. Confirm the curve is sampled at 10 kHz or resample before computing.
• Window start drifting to a later lobe: expected when the secondary lobe is high enough; verify against the analytic pulse rather than assuming the primary lobe always wins.