jmrp.io / PyOctaveBand
[READ-ONLY] Mirror of https://github.com/jmrplens/PyOctaveBand. [Python3] Octave-Band and Fractional Octave-Band filter. For signal in time domain. jmrplens.github.io/PyOctaveBand/
acoustics audio filter frequency frequency-analysis frequency-domain octave python3 signal time-domain
0
Fork 0

Configure Feed

Select the types of activity you want to include in your feed.

Human vibration exposure (ISO 8041-1, ISO 2631, ISO 5349, Directive 2002/44/EC) (#94)

New `human_vibration` module implementing the ISO 8041-1:2017 master
weighting cascade (all nine weightings Wb/Wc/Wd/We/Wf/Wh/Wj/Wk/Wm),
ISO 2631-1/-2/-4 whole-body metrics (running RMS, MTVV, VDV, MSDV, crest
factor, vibration total value, energy-equivalent acceleration), the
ISO 5349-1/-2 hand-arm A(8) arithmetic and vibration-white-finger latency,
and the Directive 2002/44/EC EAV/ELV assessment. Every result object exposes
.plot(); validated to <0.05% against the ISO 8041-1 Annex B tables and the
ISO 5349-2 Annex E worked examples. Includes docs page, figures, setup
diagram, and 7 conformance checks (100%).

authored by

José M. Requena Plens and committed by
GitHub (Jul 9, 2026, 2:43 AM +0200) 603fc8d4 a6d2960b

+2120 -7
.github/images/daily_vibration_exposure.png

This is a binary file and will not be displayed.

.github/images/daily_vibration_exposure_dark.png

This is a binary file and will not be displayed.

.github/images/daily_vibration_exposure_es.png

This is a binary file and will not be displayed.

.github/images/daily_vibration_exposure_es_dark.png

This is a binary file and will not be displayed.

+1
.github/images/diagram_human_vibration.svg
··· 1 + <svg xmlns="http://www.w3.org/2000/svg" width="900" height="580" viewBox="0 0 900 580"><rect width="900" height="580" fill="#ffffff"/><text x="450.0" y="30" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="26" font-weight="600" fill="#1a1a1a" text-anchor="middle">Whole-body vibration measurement chain (ISO 2631-1 / ISO 8041-1)</text><line x1="40" y1="510.0" x2="350" y2="510.0" stroke="#1a1a1a" stroke-width="2.2" stroke-linecap="round"/><line x1="40" y1="510.0" x2="32" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="64" y1="510.0" x2="56" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="88" y1="510.0" x2="80" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="112" y1="510.0" x2="104" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="136" y1="510.0" x2="128" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="160" y1="510.0" x2="152" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="184" y1="510.0" x2="176" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="208" y1="510.0" x2="200" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="232" y1="510.0" x2="224" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="256" y1="510.0" x2="248" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="280" y1="510.0" x2="272" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="304" y1="510.0" x2="296" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="328" y1="510.0" x2="320" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><rect x="118" y="424" width="132" height="18" rx="4" fill="#f0f2f5" stroke="#1a1a1a" stroke-width="2"/><rect x="118" y="336" width="16" height="90" rx="3" fill="#f0f2f5" stroke="#1a1a1a" stroke-width="2"/><line x1="184" y1="442" x2="184" y2="510.0" stroke="#1a1a1a" stroke-width="2.4" stroke-linecap="round"/><line x1="184" y1="506.0" x2="184.0" y2="461.0" stroke="#d62728" stroke-width="2.4" stroke-linecap="round"/><path d="M 184.0 452.0 L 187.6 461.0 L 180.4 461.0 Z" fill="#d62728" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184" y="498.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="17" fill="#d62728" text-anchor="middle" font-style="italic">vibration input</text><circle cx="178" cy="351.6" r="17.6" fill="#666666" stroke="none" stroke-width="1.5"/><line x1="178" y1="369.2" x2="178" y2="430.8" stroke="#666666" stroke-width="3" stroke-linecap="round"/><line x1="178" y1="430.8" x2="230.8" y2="430.8" stroke="#666666" stroke-width="2.4" stroke-linecap="round"/><line x1="230.8" y1="430.8" x2="230.8" y2="510.0" stroke="#666666" stroke-width="2.4" stroke-linecap="round"/><line x1="178" y1="386.8" x2="216.72" y2="413.2" stroke="#666666" stroke-width="2.4" stroke-linecap="round"/><rect x="167.0" y="412.0" width="18" height="16" rx="2" fill="#d62728" stroke="#1a1a1a" stroke-width="1.5"/><line x1="176.0" y1="412.0" x2="176.0" y2="371.0" stroke="#2ca02c" stroke-width="2.0" stroke-linecap="round"/><path d="M 176.0 362.0 L 179.6 371.0 L 172.4 371.0 Z" fill="#2ca02c" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184.0" y="366.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#2ca02c" text-anchor="start" font-weight="600">z</text><line x1="185.0" y1="420.0" x2="229.0" y2="420.0" stroke="#2ca02c" stroke-width="2.0" stroke-linecap="round"/><path d="M 238.0 420.0 L 229.0 423.6 L 229.0 416.4 Z" fill="#2ca02c" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="242.0" y="425.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#2ca02c" text-anchor="start" font-weight="600">x</text><line x1="169.0" y1="426.0" x2="139.17665747042085" y2="448.5690159683302" stroke="#2ca02c" stroke-width="2.0" stroke-linecap="round"/><path d="M 132.0 454.0 L 137.0 445.7 L 141.3 451.4 Z" fill="#2ca02c" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="124.0" y="464.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#2ca02c" text-anchor="end" font-weight="600">y</text><text x="150" y="544.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#1a1a1a" text-anchor="middle">Seat/body interface</text><rect x="490.0" y="96.0" width="320.0" height="72.0" rx="12" fill="#f0f2f5" stroke="#1f77b4" stroke-width="2"/><text x="650.0" y="127.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#1a1a1a" text-anchor="middle" font-weight="600">Triaxial accelerometer</text><text x="650.0" y="152.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#666666" text-anchor="middle">a_x , a_y , a_z (m/s²)</text><rect x="490.0" y="206.0" width="320.0" height="72.0" rx="12" fill="#f0f2f5" stroke="#1f77b4" stroke-width="2"/><text x="650.0" y="237.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#1a1a1a" text-anchor="middle" font-weight="600">Band limiting + Wk / Wd</text><text x="650.0" y="262.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#666666" text-anchor="middle">weighting (ISO 8041-1)</text><rect x="490.0" y="316.0" width="320.0" height="72.0" rx="12" fill="#f0f2f5" stroke="#1f77b4" stroke-width="2"/><text x="650.0" y="347.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#1a1a1a" text-anchor="middle" font-weight="600">Weighted r.m.s. a_w &amp; VDV</text><text x="650.0" y="372.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#666666" text-anchor="middle">(ISO 2631-1)</text><line x1="650.0" y1="168" x2="650.0" y2="197.0" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 206.0 L 646.4 197.0 L 653.6 197.0 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="278" x2="650.0" y2="307.0" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 316.0 L 646.4 307.0 L 653.6 307.0 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="252" y1="420.0" x2="478.3540341591462" y2="139.00878518174954" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 484.0 132.0 L 481.2 141.3 L 475.6 136.8 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="388" x2="650.0" y2="415.0" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 424.0 L 646.4 415.0 L 653.6 415.0 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><rect x="400" y="424" width="470" height="78" rx="12" fill="none" stroke="#d62728" stroke-width="2" stroke-dasharray="6,5"/><text x="635" y="452" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="20" fill="#1a1a1a" text-anchor="middle" font-weight="600">a_v = √(Σ k_j² a_wj²) → A(8) = a_v·√(T/T₀)</text><text x="635" y="480" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#d62728" text-anchor="middle">assessed vs EAV / ELV (Directive 2002/44/EC)</text></svg>
+1
.github/images/diagram_human_vibration_dark.svg
··· 1 + <svg xmlns="http://www.w3.org/2000/svg" width="900" height="580" viewBox="0 0 900 580"><rect width="900" height="580" fill="#0d1117"/><text x="450.0" y="30" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="26" font-weight="600" fill="#e6e6e6" text-anchor="middle">Whole-body vibration measurement chain (ISO 2631-1 / ISO 8041-1)</text><line x1="40" y1="510.0" x2="350" y2="510.0" stroke="#e6e6e6" stroke-width="2.2" stroke-linecap="round"/><line x1="40" y1="510.0" x2="32" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="64" y1="510.0" x2="56" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="88" y1="510.0" x2="80" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="112" y1="510.0" x2="104" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="136" y1="510.0" x2="128" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="160" y1="510.0" x2="152" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="184" y1="510.0" x2="176" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="208" y1="510.0" x2="200" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="232" y1="510.0" x2="224" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="256" y1="510.0" x2="248" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="280" y1="510.0" x2="272" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="304" y1="510.0" x2="296" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="328" y1="510.0" x2="320" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><rect x="118" y="424" width="132" height="18" rx="4" fill="#1c2128" stroke="#e6e6e6" stroke-width="2"/><rect x="118" y="336" width="16" height="90" rx="3" fill="#1c2128" stroke="#e6e6e6" stroke-width="2"/><line x1="184" y1="442" x2="184" y2="510.0" stroke="#e6e6e6" stroke-width="2.4" stroke-linecap="round"/><line x1="184" y1="506.0" x2="184.0" y2="461.0" stroke="#e46a6a" stroke-width="2.4" stroke-linecap="round"/><path d="M 184.0 452.0 L 187.6 461.0 L 180.4 461.0 Z" fill="#e46a6a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184" y="498.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="17" fill="#e46a6a" text-anchor="middle" font-style="italic">vibration input</text><circle cx="178" cy="351.6" r="17.6" fill="#9a9a9a" stroke="none" stroke-width="1.5"/><line x1="178" y1="369.2" x2="178" y2="430.8" stroke="#9a9a9a" stroke-width="3" stroke-linecap="round"/><line x1="178" y1="430.8" x2="230.8" y2="430.8" stroke="#9a9a9a" stroke-width="2.4" stroke-linecap="round"/><line x1="230.8" y1="430.8" x2="230.8" y2="510.0" stroke="#9a9a9a" stroke-width="2.4" stroke-linecap="round"/><line x1="178" y1="386.8" x2="216.72" y2="413.2" stroke="#9a9a9a" stroke-width="2.4" stroke-linecap="round"/><rect x="167.0" y="412.0" width="18" height="16" rx="2" fill="#e46a6a" stroke="#e6e6e6" stroke-width="1.5"/><line x1="176.0" y1="412.0" x2="176.0" y2="371.0" stroke="#5abf5a" stroke-width="2.0" stroke-linecap="round"/><path d="M 176.0 362.0 L 179.6 371.0 L 172.4 371.0 Z" fill="#5abf5a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184.0" y="366.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#5abf5a" text-anchor="start" font-weight="600">z</text><line x1="185.0" y1="420.0" x2="229.0" y2="420.0" stroke="#5abf5a" stroke-width="2.0" stroke-linecap="round"/><path d="M 238.0 420.0 L 229.0 423.6 L 229.0 416.4 Z" fill="#5abf5a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="242.0" y="425.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#5abf5a" text-anchor="start" font-weight="600">x</text><line x1="169.0" y1="426.0" x2="139.17665747042085" y2="448.5690159683302" stroke="#5abf5a" stroke-width="2.0" stroke-linecap="round"/><path d="M 132.0 454.0 L 137.0 445.7 L 141.3 451.4 Z" fill="#5abf5a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="124.0" y="464.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#5abf5a" text-anchor="end" font-weight="600">y</text><text x="150" y="544.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#e6e6e6" text-anchor="middle">Seat/body interface</text><rect x="490.0" y="96.0" width="320.0" height="72.0" rx="12" fill="#1c2128" stroke="#4da3d8" stroke-width="2"/><text x="650.0" y="127.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#e6e6e6" text-anchor="middle" font-weight="600">Triaxial accelerometer</text><text x="650.0" y="152.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#9a9a9a" text-anchor="middle">a_x , a_y , a_z (m/s²)</text><rect x="490.0" y="206.0" width="320.0" height="72.0" rx="12" fill="#1c2128" stroke="#4da3d8" stroke-width="2"/><text x="650.0" y="237.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#e6e6e6" text-anchor="middle" font-weight="600">Band limiting + Wk / Wd</text><text x="650.0" y="262.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#9a9a9a" text-anchor="middle">weighting (ISO 8041-1)</text><rect x="490.0" y="316.0" width="320.0" height="72.0" rx="12" fill="#1c2128" stroke="#4da3d8" stroke-width="2"/><text x="650.0" y="347.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#e6e6e6" text-anchor="middle" font-weight="600">Weighted r.m.s. a_w &amp; VDV</text><text x="650.0" y="372.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#9a9a9a" text-anchor="middle">(ISO 2631-1)</text><line x1="650.0" y1="168" x2="650.0" y2="197.0" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 206.0 L 646.4 197.0 L 653.6 197.0 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="278" x2="650.0" y2="307.0" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 316.0 L 646.4 307.0 L 653.6 307.0 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="252" y1="420.0" x2="478.3540341591462" y2="139.00878518174954" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 484.0 132.0 L 481.2 141.3 L 475.6 136.8 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="388" x2="650.0" y2="415.0" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 424.0 L 646.4 415.0 L 653.6 415.0 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><rect x="400" y="424" width="470" height="78" rx="12" fill="none" stroke="#e46a6a" stroke-width="2" stroke-dasharray="6,5"/><text x="635" y="452" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="20" fill="#e6e6e6" text-anchor="middle" font-weight="600">a_v = √(Σ k_j² a_wj²) → A(8) = a_v·√(T/T₀)</text><text x="635" y="480" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#e46a6a" text-anchor="middle">assessed vs EAV / ELV (Directive 2002/44/EC)</text></svg>
+1
.github/images/diagram_human_vibration_es.svg
··· 1 + <svg xmlns="http://www.w3.org/2000/svg" width="900" height="580" viewBox="0 0 900 580"><rect width="900" height="580" fill="#ffffff"/><text x="450.0" y="30" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="26" font-weight="600" fill="#1a1a1a" text-anchor="middle">Cadena de medición de vibración de cuerpo entero (ISO 2631-1 / ISO 8041-1)</text><line x1="40" y1="510.0" x2="350" y2="510.0" stroke="#1a1a1a" stroke-width="2.2" stroke-linecap="round"/><line x1="40" y1="510.0" x2="32" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="64" y1="510.0" x2="56" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="88" y1="510.0" x2="80" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="112" y1="510.0" x2="104" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="136" y1="510.0" x2="128" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="160" y1="510.0" x2="152" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="184" y1="510.0" x2="176" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="208" y1="510.0" x2="200" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="232" y1="510.0" x2="224" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="256" y1="510.0" x2="248" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="280" y1="510.0" x2="272" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="304" y1="510.0" x2="296" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><line x1="328" y1="510.0" x2="320" y2="519.0" stroke="#666666" stroke-width="1.1" stroke-linecap="round"/><rect x="118" y="424" width="132" height="18" rx="4" fill="#f0f2f5" stroke="#1a1a1a" stroke-width="2"/><rect x="118" y="336" width="16" height="90" rx="3" fill="#f0f2f5" stroke="#1a1a1a" stroke-width="2"/><line x1="184" y1="442" x2="184" y2="510.0" stroke="#1a1a1a" stroke-width="2.4" stroke-linecap="round"/><line x1="184" y1="506.0" x2="184.0" y2="461.0" stroke="#d62728" stroke-width="2.4" stroke-linecap="round"/><path d="M 184.0 452.0 L 187.6 461.0 L 180.4 461.0 Z" fill="#d62728" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184" y="498.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="17" fill="#d62728" text-anchor="middle" font-style="italic">entrada de vibración</text><circle cx="178" cy="351.6" r="17.6" fill="#666666" stroke="none" stroke-width="1.5"/><line x1="178" y1="369.2" x2="178" y2="430.8" stroke="#666666" stroke-width="3" stroke-linecap="round"/><line x1="178" y1="430.8" x2="230.8" y2="430.8" stroke="#666666" stroke-width="2.4" stroke-linecap="round"/><line x1="230.8" y1="430.8" x2="230.8" y2="510.0" stroke="#666666" stroke-width="2.4" stroke-linecap="round"/><line x1="178" y1="386.8" x2="216.72" y2="413.2" stroke="#666666" stroke-width="2.4" stroke-linecap="round"/><rect x="167.0" y="412.0" width="18" height="16" rx="2" fill="#d62728" stroke="#1a1a1a" stroke-width="1.5"/><line x1="176.0" y1="412.0" x2="176.0" y2="371.0" stroke="#2ca02c" stroke-width="2.0" stroke-linecap="round"/><path d="M 176.0 362.0 L 179.6 371.0 L 172.4 371.0 Z" fill="#2ca02c" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184.0" y="366.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#2ca02c" text-anchor="start" font-weight="600">z</text><line x1="185.0" y1="420.0" x2="229.0" y2="420.0" stroke="#2ca02c" stroke-width="2.0" stroke-linecap="round"/><path d="M 238.0 420.0 L 229.0 423.6 L 229.0 416.4 Z" fill="#2ca02c" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="242.0" y="425.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#2ca02c" text-anchor="start" font-weight="600">x</text><line x1="169.0" y1="426.0" x2="139.17665747042085" y2="448.5690159683302" stroke="#2ca02c" stroke-width="2.0" stroke-linecap="round"/><path d="M 132.0 454.0 L 137.0 445.7 L 141.3 451.4 Z" fill="#2ca02c" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="124.0" y="464.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#2ca02c" text-anchor="end" font-weight="600">y</text><text x="150" y="544.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#1a1a1a" text-anchor="middle">Interfaz asiento/cuerpo</text><rect x="490.0" y="96.0" width="320.0" height="72.0" rx="12" fill="#f0f2f5" stroke="#1f77b4" stroke-width="2"/><text x="650.0" y="127.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#1a1a1a" text-anchor="middle" font-weight="600">Acelerómetro triaxial</text><text x="650.0" y="152.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#666666" text-anchor="middle">a_x , a_y , a_z (m/s²)</text><rect x="490.0" y="206.0" width="320.0" height="72.0" rx="12" fill="#f0f2f5" stroke="#1f77b4" stroke-width="2"/><text x="650.0" y="237.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#1a1a1a" text-anchor="middle" font-weight="600">Limitación de banda + Wk / Wd</text><text x="650.0" y="262.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#666666" text-anchor="middle">ponderación (ISO 8041-1)</text><rect x="490.0" y="316.0" width="320.0" height="72.0" rx="12" fill="#f0f2f5" stroke="#1f77b4" stroke-width="2"/><text x="650.0" y="347.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#1a1a1a" text-anchor="middle" font-weight="600">a_w eficaz ponderada y VDV</text><text x="650.0" y="372.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#666666" text-anchor="middle">(ISO 2631-1)</text><line x1="650.0" y1="168" x2="650.0" y2="197.0" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 206.0 L 646.4 197.0 L 653.6 197.0 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="278" x2="650.0" y2="307.0" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 316.0 L 646.4 307.0 L 653.6 307.0 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="252" y1="420.0" x2="478.3540341591462" y2="139.00878518174954" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 484.0 132.0 L 481.2 141.3 L 475.6 136.8 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="388" x2="650.0" y2="415.0" stroke="#1a1a1a" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 424.0 L 646.4 415.0 L 653.6 415.0 Z" fill="#1a1a1a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><rect x="400" y="424" width="470" height="78" rx="12" fill="none" stroke="#d62728" stroke-width="2" stroke-dasharray="6,5"/><text x="635" y="452" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="20" fill="#1a1a1a" text-anchor="middle" font-weight="600">a_v = √(Σ k_j² a_wj²) → A(8) = a_v·√(T/T₀)</text><text x="635" y="480" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#d62728" text-anchor="middle">evaluada frente a EAV / ELV (Directiva 2002/44/CE)</text></svg>
+1
.github/images/diagram_human_vibration_es_dark.svg
··· 1 + <svg xmlns="http://www.w3.org/2000/svg" width="900" height="580" viewBox="0 0 900 580"><rect width="900" height="580" fill="#0d1117"/><text x="450.0" y="30" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="26" font-weight="600" fill="#e6e6e6" text-anchor="middle">Cadena de medición de vibración de cuerpo entero (ISO 2631-1 / ISO 8041-1)</text><line x1="40" y1="510.0" x2="350" y2="510.0" stroke="#e6e6e6" stroke-width="2.2" stroke-linecap="round"/><line x1="40" y1="510.0" x2="32" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="64" y1="510.0" x2="56" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="88" y1="510.0" x2="80" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="112" y1="510.0" x2="104" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="136" y1="510.0" x2="128" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="160" y1="510.0" x2="152" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="184" y1="510.0" x2="176" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="208" y1="510.0" x2="200" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="232" y1="510.0" x2="224" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="256" y1="510.0" x2="248" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="280" y1="510.0" x2="272" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="304" y1="510.0" x2="296" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><line x1="328" y1="510.0" x2="320" y2="519.0" stroke="#9a9a9a" stroke-width="1.1" stroke-linecap="round"/><rect x="118" y="424" width="132" height="18" rx="4" fill="#1c2128" stroke="#e6e6e6" stroke-width="2"/><rect x="118" y="336" width="16" height="90" rx="3" fill="#1c2128" stroke="#e6e6e6" stroke-width="2"/><line x1="184" y1="442" x2="184" y2="510.0" stroke="#e6e6e6" stroke-width="2.4" stroke-linecap="round"/><line x1="184" y1="506.0" x2="184.0" y2="461.0" stroke="#e46a6a" stroke-width="2.4" stroke-linecap="round"/><path d="M 184.0 452.0 L 187.6 461.0 L 180.4 461.0 Z" fill="#e46a6a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184" y="498.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="17" fill="#e46a6a" text-anchor="middle" font-style="italic">entrada de vibración</text><circle cx="178" cy="351.6" r="17.6" fill="#9a9a9a" stroke="none" stroke-width="1.5"/><line x1="178" y1="369.2" x2="178" y2="430.8" stroke="#9a9a9a" stroke-width="3" stroke-linecap="round"/><line x1="178" y1="430.8" x2="230.8" y2="430.8" stroke="#9a9a9a" stroke-width="2.4" stroke-linecap="round"/><line x1="230.8" y1="430.8" x2="230.8" y2="510.0" stroke="#9a9a9a" stroke-width="2.4" stroke-linecap="round"/><line x1="178" y1="386.8" x2="216.72" y2="413.2" stroke="#9a9a9a" stroke-width="2.4" stroke-linecap="round"/><rect x="167.0" y="412.0" width="18" height="16" rx="2" fill="#e46a6a" stroke="#e6e6e6" stroke-width="1.5"/><line x1="176.0" y1="412.0" x2="176.0" y2="371.0" stroke="#5abf5a" stroke-width="2.0" stroke-linecap="round"/><path d="M 176.0 362.0 L 179.6 371.0 L 172.4 371.0 Z" fill="#5abf5a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="184.0" y="366.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#5abf5a" text-anchor="start" font-weight="600">z</text><line x1="185.0" y1="420.0" x2="229.0" y2="420.0" stroke="#5abf5a" stroke-width="2.0" stroke-linecap="round"/><path d="M 238.0 420.0 L 229.0 423.6 L 229.0 416.4 Z" fill="#5abf5a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="242.0" y="425.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#5abf5a" text-anchor="start" font-weight="600">x</text><line x1="169.0" y1="426.0" x2="139.17665747042085" y2="448.5690159683302" stroke="#5abf5a" stroke-width="2.0" stroke-linecap="round"/><path d="M 132.0 454.0 L 137.0 445.7 L 141.3 451.4 Z" fill="#5abf5a" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><text x="124.0" y="464.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#5abf5a" text-anchor="end" font-weight="600">y</text><text x="150" y="544.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#e6e6e6" text-anchor="middle">Interfaz asiento/cuerpo</text><rect x="490.0" y="96.0" width="320.0" height="72.0" rx="12" fill="#1c2128" stroke="#4da3d8" stroke-width="2"/><text x="650.0" y="127.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#e6e6e6" text-anchor="middle" font-weight="600">Acelerómetro triaxial</text><text x="650.0" y="152.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#9a9a9a" text-anchor="middle">a_x , a_y , a_z (m/s²)</text><rect x="490.0" y="206.0" width="320.0" height="72.0" rx="12" fill="#1c2128" stroke="#4da3d8" stroke-width="2"/><text x="650.0" y="237.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#e6e6e6" text-anchor="middle" font-weight="600">Limitación de banda + Wk / Wd</text><text x="650.0" y="262.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#9a9a9a" text-anchor="middle">ponderación (ISO 8041-1)</text><rect x="490.0" y="316.0" width="320.0" height="72.0" rx="12" fill="#1c2128" stroke="#4da3d8" stroke-width="2"/><text x="650.0" y="347.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="21" fill="#e6e6e6" text-anchor="middle" font-weight="600">a_w eficaz ponderada y VDV</text><text x="650.0" y="372.0" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#9a9a9a" text-anchor="middle">(ISO 2631-1)</text><line x1="650.0" y1="168" x2="650.0" y2="197.0" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 206.0 L 646.4 197.0 L 653.6 197.0 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="278" x2="650.0" y2="307.0" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 316.0 L 646.4 307.0 L 653.6 307.0 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="252" y1="420.0" x2="478.3540341591462" y2="139.00878518174954" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 484.0 132.0 L 481.2 141.3 L 475.6 136.8 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><line x1="650.0" y1="388" x2="650.0" y2="415.0" stroke="#e6e6e6" stroke-width="2.0" stroke-linecap="round"/><path d="M 650.0 424.0 L 646.4 415.0 L 653.6 415.0 Z" fill="#e6e6e6" stroke="none" stroke-width="1.5" stroke-linejoin="round"/><rect x="400" y="424" width="470" height="78" rx="12" fill="none" stroke="#e46a6a" stroke-width="2" stroke-dasharray="6,5"/><text x="635" y="452" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="20" fill="#e6e6e6" text-anchor="middle" font-weight="600">a_v = √(Σ k_j² a_wj²) → A(8) = a_v·√(T/T₀)</text><text x="635" y="480" font-family="Segoe UI, Helvetica, Arial, sans-serif" font-size="18" fill="#e46a6a" text-anchor="middle">evaluada frente a EAV / ELV (Directiva 2002/44/CE)</text></svg>
.github/images/vibration_weighting.png

This is a binary file and will not be displayed.

.github/images/vibration_weighting_dark.png

This is a binary file and will not be displayed.

.github/images/vibration_weighting_es.png

This is a binary file and will not be displayed.

.github/images/vibration_weighting_es_dark.png

This is a binary file and will not be displayed.

.github/images/weighted_acceleration.png

This is a binary file and will not be displayed.

.github/images/weighted_acceleration_dark.png

This is a binary file and will not be displayed.

.github/images/weighted_acceleration_es.png

This is a binary file and will not be displayed.

.github/images/weighted_acceleration_es_dark.png

This is a binary file and will not be displayed.

+16 -1
docs/CONFORMANCE.md
··· 15 15 16 16 ## Numerical conformance report 17 17 18 - &#9989; **60/60 conformance checks pass** across 12 domains and 37 standards - filters class 1 - weightings within IEC 61672-1 class 1. 18 + &#9989; **67/67 conformance checks pass** across 13 domains and 41 standards - filters class 1 - weightings within IEC 61672-1 class 1. 19 19 20 20 ### Numerical validation - filters &amp; weightings 21 21 ··· 196 196 197 197 </details> 198 198 199 + <details> 200 + <summary>&#9989; <b>Human vibration (ISO 8041 / 2631 / 5349)</b> — 100% (7/7)</summary> 201 + 202 + | Standard | Quantity | Expected (norm) | Computed | &#916; | Status | 203 + |:---|:---|:---|:---|:---|:---:| 204 + | ISO 8041-1:2017 Table B.8 | Wk design-goal factor at 6,31 Hz | 1.054 (+/-0.1%) | 1.0544 | 0 | &#9989; | 205 + | ISO 8041-1:2017 Table B.9 | Wm design-goal factor at 1,585 Hz | 0.9342 (+/-0.1%) | 0.9342 | 0 | &#9989; | 206 + | ISO 8041-1:2017 Table 1 | Wh factor at the 500 rad/s reference | 0.202 (+/-0.15%) | 0.202 | 0 | &#9989; | 207 + | ISO 5349-2:2001 Example E.2.1 | Single-tool daily exposure A(8) | 4.1 m/s^2 (+/-0.05 m/s^2) | 4.14 m/s^2 | 0.037 m/s^2 | &#9989; | 208 + | ISO 5349-2:2001 Example E.3 | Forestry three-task A(8) | 3.6 m/s^2 (+/-0.05 m/s^2) | 3.61 m/s^2 | 0.01 m/s^2 | &#9989; | 209 + | ISO 5349-1:2001 Eq. (C.1) | VWF 10 % lifetime Dy at A(8)=7 | 4 yr (+/-0.1 yr) | 4.04 yr | 0.042 yr | &#9989; | 210 + | Directive 2002/44/EC Art. 3 | HAV/WBV action & limit values | HAV 2.5/5.0, WBV 0.5/1.15 m/s^2 | HAV 2.5/5.0, WBV 0.5/1.15 m/s^2 | 0 | &#9989; | 211 + 212 + </details> 213 +
+1
docs/README.md
··· 17 17 - [Sound Power](sound-power.md) — sound power level by enveloping surface (ISO 3744/3746), reverberation room (ISO 3741), intensity scanning (ISO 9614-2), and the precision grades in an anechoic room (ISO 3745) and by precision intensity scanning (ISO 9614-3) 18 18 - [Acoustic Materials](materials.md) — sound-absorption rating α_w and classes (ISO 11654), airflow resistance static and alternating methods (ISO 9053-1/-2), and impedance-tube measurement of absorption, surface impedance and transmission loss (ISO 10534-1/-2, ASTM E2611) 19 19 - [Surface Scattering, Diffusion and In-situ Absorption](surface-scattering.md) — random-incidence scattering (ISO 17497-1), free-field diffusion coefficient (ISO 17497-2), and in-situ road-surface absorption by the extended-surface subtraction technique (ISO 13472-1) and the spot method (ISO 13472-2) 20 + - [Human Vibration](human-vibration.md) — whole-body and hand-arm frequency weightings (ISO 8041-1), weighted r.m.s. acceleration, running r.m.s./MTVV/VDV and crest factor (ISO 2631-1), vibration in buildings (ISO 2631-2), vibration total value and daily exposure A(8) (ISO 5349-1/-2), and the exposure action/limit values of Directive 2002/44/EC 20 21 - [Calibration and dBFS](calibration.md) — physical SPL and digital analysis 21 22 - [Block Processing](block-processing.md) — stateful real-time workflows 22 23 - [Multichannel](multichannel.md) — vectorized multichannel analysis
+271
docs/human-vibration.md
··· 1 + ← [Documentation index](README.md) 2 + 3 + # Human Vibration — Whole-Body and Hand-Arm Exposure 4 + 5 + Vibration transmitted to a person is evaluated with the same measurement chain 6 + whatever its origin: the acceleration is **frequency-weighted** to reflect how 7 + the body responds at each frequency, reduced to a **weighted r.m.s.** 8 + acceleration (with dose measures for shocks and long records), combined across 9 + axes into a **vibration total value**, and finally normalised to an **8-hour 10 + daily exposure** `A(8)` that is compared against the action and limit values of 11 + the European directive. 12 + 13 + The weightings themselves are defined once, in **ISO 8041-1:2017**, as a cascade 14 + of analog filters; **ISO 2631-1** applies them to whole-body vibration, 15 + **ISO 2631-2** to vibration in buildings, **ISO 2631-4** to rail ride comfort, 16 + and **ISO 5349-1/-2** to hand-transmitted vibration. This page covers the whole 17 + chain. 18 + 19 + <picture><source media="(prefers-color-scheme: dark)" srcset="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/diagram_human_vibration_dark.svg"><img src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/diagram_human_vibration.svg" alt="Whole-body vibration measurement chain: a triaxial accelerometer at the seat/body interface of a seated person measures the x, y and z acceleration; each axis is band-limited and frequency-weighted (Wk vertical, Wd horizontal) per ISO 8041-1, reduced to a weighted r.m.s. a_w and VDV per ISO 2631-1, combined into the vibration total value a_v, normalised to the daily exposure A(8) and assessed against the EAV and ELV of Directive 2002/44/EC" width="94%"></picture> 20 + 21 + ## 1. Frequency weightings (ISO 8041-1) 22 + 23 + Every human-vibration weighting is the product of four analog stages evaluated 24 + at `s = j·2πf` (ISO 8041-1 Formulae (1)–(5)): a second-order Butterworth 25 + **high-pass** and **low-pass** band limiting, an **acceleration–velocity 26 + transition** carrying the overall gain `K`, and an **upward step**: 27 + 28 + $$ 29 + H(s) = H_h(s)\,H_l(s)\,H_t(s)\,H_s(s). 30 + $$ 31 + 32 + A single Table 3 parameter set `(f_1,\,Q_1,\,f_2,\,Q_2,\,f_3,\,f_4,\,Q_4,\, 33 + f_5,\,Q_5,\,f_6,\,Q_6,\,K)` realises all nine weightings — `Wb, Wc, Wd, We, Wf, 34 + Wh, Wj, Wk, Wm` — a corner set to infinity collapsing its stage to unity. The 35 + principal whole-body weighting is `Wk` (vertical, seat surface); `Wd` is the 36 + horizontal weighting, and `Wh` the hand-arm weighting. 37 + 38 + ```python 39 + import phonometry as ph 40 + 41 + # The overall weighting response at any frequencies (ISO 8041-1 Formula (5)). 42 + resp = ph.frequency_weighting("Wk", [1.0, 6.3096, 20.0]) 43 + print(resp.magnitude.round(3)) # [0.482 1.054 0.636] (factors) 44 + print(resp.magnitude_db.round(2)) # [-6.33 0.46 -3.93] (dB) 45 + ``` 46 + 47 + The factors reproduce the ISO 8041-1 Annex B design-goal tables to their four 48 + significant figures: `Wk` plateaus near −6 dB below 2 Hz, peaks at +0.46 dB near 49 + 6.3 Hz and rolls off above. 50 + 51 + <picture><source media="(prefers-color-scheme: dark)" srcset="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/vibration_weighting_dark.png"><img src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/vibration_weighting.png" alt="The whole-body vertical weighting Wk in decibels over 0.4 to 100 Hz: a plateau near -6 dB below 2 Hz, a small +0.5 dB peak near 6 Hz and a roll-off to about -21 dB at 100 Hz" width="88%"></picture> 52 + 53 + <details> 54 + <summary>Show the code for this figure</summary> 55 + 56 + ```python 57 + import numpy as np 58 + import matplotlib.pyplot as plt 59 + import phonometry as ph 60 + 61 + result = ph.frequency_weighting("Wk", np.geomspace(0.4, 100.0, 240)) 62 + 63 + # One line: 64 + result.plot() 65 + plt.show() 66 + 67 + # By hand, from the result's fields, mirroring what WeightingResponse.plot() draws: 68 + fig, ax = plt.subplots() 69 + ax.semilogx(result.frequencies, result.magnitude_db, color="#1f77b4") 70 + ax.set_xlabel("Frequency [Hz]") 71 + ax.set_ylabel("Weighting factor [dB]") 72 + ax.set_title("Whole-body vertical weighting Wk (ISO 8041-1)") 73 + plt.show() 74 + ``` 75 + 76 + </details> 77 + 78 + To weight a time signal, `apply_weighting` applies the exact complex response in 79 + the frequency domain (so magnitude *and* phase match the standard), which the 80 + time-domain dose metrics below then consume. 81 + 82 + ## 2. Weighted acceleration and dose measures (ISO 2631-1) 83 + 84 + The basic evaluation is the **weighted r.m.s. acceleration**. From a 85 + one-third-octave spectrum it is (ISO 2631-1 Eq. (9); the identical construction 86 + gives the hand-arm `a_hw` of ISO 5349-1 Eq. (A.1)): 87 + 88 + $$ 89 + a_w = \sqrt{\sum_i \left(W_i\,a_i\right)^2}, 90 + $$ 91 + 92 + with `W_i` the weighting factor at band centre `i` and `a_i` the measured band 93 + acceleration. 94 + 95 + ```python 96 + import numpy as np 97 + import phonometry as ph 98 + 99 + # A measured vertical seat spectrum (r.m.s. per one-third octave, m/s^2). 100 + freqs = np.array([1.0, 2.0, 4.0, 8.0, 16.0, 31.5, 63.0]) 101 + accel = np.array([0.20, 0.45, 0.42, 0.25, 0.12, 0.05, 0.02]) 102 + 103 + result = ph.weighted_acceleration(accel, freqs, "Wk") 104 + print(round(result.overall, 3)) # 0.555 m/s^2 (a_w) 105 + print(result.weighted.round(3)) # W_i * a_i per band 106 + ``` 107 + 108 + <picture><source media="(prefers-color-scheme: dark)" srcset="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/weighted_acceleration_dark.png"><img src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/weighted_acceleration.png" alt="A measured vehicle-seat acceleration spectrum (grey) and its Wk-weighted contribution (blue) over the one-third octaves from 1 to 80 Hz: the weighting attenuates the low and high bands but leaves the 4 to 8 Hz range nearly unchanged, giving a weighted r.m.s. a_w of about 1.03 m/s^2" width="90%"></picture> 109 + 110 + <details> 111 + <summary>Show the code for this figure</summary> 112 + 113 + ```python 114 + import numpy as np 115 + import matplotlib.pyplot as plt 116 + import phonometry as ph 117 + 118 + freqs = np.array([1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 119 + 12.5, 16.0, 20.0, 25.0, 31.5, 40.0, 63.0, 80.0]) 120 + accel = np.array([0.18, 0.24, 0.33, 0.46, 0.52, 0.55, 0.48, 0.39, 0.31, 0.26, 121 + 0.21, 0.17, 0.13, 0.10, 0.078, 0.060, 0.045, 0.028, 0.020]) 122 + result = ph.weighted_acceleration(accel, freqs, "Wk") 123 + 124 + # One line: 125 + result.plot() 126 + plt.show() 127 + 128 + # By hand, mirroring what WeightedSpectrum.plot() draws: 129 + pos = np.arange(freqs.size) 130 + fig, ax = plt.subplots() 131 + ax.bar(pos - 0.2, result.band_accelerations, 0.4, color="#bbbbbb", 132 + label="Unweighted $a_i$") 133 + ax.bar(pos + 0.2, result.weighted, 0.4, color="#1f77b4", 134 + label="Weighted $W_i a_i$ (Wk)") 135 + ax.set_xticks(pos) 136 + ax.set_xticklabels([f"{f:g}" for f in freqs], rotation=45, ha="right") 137 + ax.set_xlabel("Frequency [Hz]") 138 + ax.set_ylabel(r"r.m.s. acceleration [m/s$^2$]") 139 + ax.set_title(f"Weighted acceleration ($a_w$ = {result.overall:.3f} m/s²)") 140 + ax.legend() 141 + plt.show() 142 + ``` 143 + 144 + </details> 145 + 146 + When the r.m.s. value understates an intermittent or shock-laden exposure, 147 + ISO 2631-1 adds dose measures computed on the weighted time signal: the 148 + **running r.m.s.** and its maximum, the **maximum transient vibration value** 149 + `MTVV` (Eq. (4), a 1 s running r.m.s.); the fourth-power **vibration dose 150 + value** `VDV = (∫ a_w^4\,dt)^{1/4}` (Eq. (5)); the **motion sickness dose value** 151 + `MSDV = (∫ a_w^2\,dt)^{1/2}`; and the **crest factor** (peak / r.m.s.), whose 152 + value above 9 signals that the basic method is inadequate. 153 + 154 + ```python 155 + import numpy as np 156 + import phonometry as ph 157 + 158 + fs = 1000.0 159 + raw = np.random.default_rng(0).standard_normal(int(60 * fs)) # 60 s record 160 + a_w = ph.apply_weighting(raw, fs, "Wk") # weighted signal 161 + 162 + print(round(ph.vibration_dose_value(a_w, fs), 3)) # VDV [m/s^1.75] 163 + print(round(ph.mtvv(a_w, fs), 3)) # MTVV [m/s^2] 164 + print(round(ph.crest_factor(a_w), 2)) # crest factor 165 + ``` 166 + 167 + ## 3. Vibration total value and daily exposure `A(8)` 168 + 169 + Across the three axes the **vibration total value** combines the axis-weighted 170 + r.m.s. accelerations with the posture multiplying factors `k_j` (ISO 2631-1 171 + Eq. (10); for hand-arm, ISO 5349-1 Eq. (1) with every `k = 1`): 172 + 173 + $$ 174 + a_v = \sqrt{\sum_j k_j^2\,a_{wj}^2}. 175 + $$ 176 + 177 + ```python 178 + import phonometry as ph 179 + 180 + # Health, seated: k = 1.4 / 1.4 / 1.0 (ISO 2631-1, 7.2.3). 181 + a_v = ph.vibration_total_value([0.35, 0.28, 0.62], k=[1.4, 1.4, 1.0]) 182 + print(round(a_v, 3)) # 0.882 m/s^2 183 + ``` 184 + 185 + The **daily exposure** normalises the total value to a reference 8-hour day 186 + (`T_0 = 28 800 s`). For a single operation `A(8) = a_v·\sqrt{T/T_0}`; several 187 + operations combine through their partial exposures `A_i(8) = a_{vi}·\sqrt{T_i/T_0}` 188 + as `A(8) = \sqrt{\sum_i A_i(8)^2}` (ISO 5349-1 Eqs. (2)/(3); ISO 5349-2 189 + Eqs. (1)–(3)). 190 + 191 + `daily_vibration_exposure` builds the partial exposures, combines them and 192 + assesses the result against **Directive 2002/44/EC** — hand-arm action value 193 + `A(8) = 2.5` and limit value `5` m/s²; whole-body action `0.5` and limit `1.15` 194 + m/s² (or a VDV of `9.1` / `21` m/s¹·⁷⁵): 195 + 196 + ```python 197 + import phonometry as ph 198 + 199 + # ISO 5349-2 Annex E.3: a forestry worker's three chain-saw tasks. 200 + result = ph.daily_vibration_exposure( 201 + total_values=[4.6, 6.0, 3.6], # a_hv per task, m/s^2 202 + durations_s=[2 * 3600, 1 * 3600, 2 * 3600], # exposure time per task 203 + kind="hav", 204 + labels=["brush-saw", "felling", "stripping"], 205 + ) 206 + print(result.partials.round(2)) # [2.3 2.12 1.8 ] A_i(8) 207 + print(round(result.a8, 2)) # 3.61 m/s^2 208 + print(result.assessment.zone) # 'action' (2.5 <= A(8) < 5.0) 209 + ``` 210 + 211 + <picture><source media="(prefers-color-scheme: dark)" srcset="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/daily_vibration_exposure_dark.png"><img src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/daily_vibration_exposure.png" alt="A bar chart of the three partial hand-arm exposures (about 2.3, 2.1 and 1.8 m/s^2) and the combined A(8) of 3.61 m/s^2, with the Directive 2002/44/EC exposure action value at 2.5 and exposure limit value at 5.0 m/s^2 marked as horizontal lines; the daily exposure sits in the action zone between them" width="82%"></picture> 212 + 213 + <details> 214 + <summary>Show the code for this figure</summary> 215 + 216 + ```python 217 + import numpy as np 218 + import matplotlib.pyplot as plt 219 + import phonometry as ph 220 + 221 + result = ph.daily_vibration_exposure( 222 + [4.6, 6.0, 3.6], [2 * 3600, 1 * 3600, 2 * 3600], kind="hav", 223 + labels=["brush-saw", "felling", "stripping"], 224 + ) 225 + 226 + # One line: 227 + result.plot() 228 + plt.show() 229 + 230 + # By hand, mirroring what DailyVibrationExposure.plot() draws: 231 + labels = [*result.labels, "A(8)"] 232 + values = [*result.partials.tolist(), result.a8] 233 + a = result.assessment 234 + fig, ax = plt.subplots() 235 + ax.bar(range(len(values)), values, 236 + color=["#bbbbbb"] * result.partials.size + ["#1f77b4"]) 237 + ax.axhline(a.action_value, color="#2ca02c", ls="--", label=f"EAV = {a.action_value:g}") 238 + ax.axhline(a.limit_value, color="#d62728", ls="--", label=f"ELV = {a.limit_value:g}") 239 + ax.set_xticks(range(len(values))) 240 + ax.set_xticklabels(labels, rotation=30, ha="right") 241 + ax.set_ylabel(r"Daily exposure A(8) [m/s$^2$]") 242 + ax.legend() 243 + plt.show() 244 + ``` 245 + 246 + </details> 247 + 248 + ## 4. Exposure–response guidance 249 + 250 + For hand-transmitted vibration, ISO 5349-1 Annex C relates the daily exposure to 251 + the group-mean lifetime `D_y` (in years) that produces vibration-white-finger in 252 + 10 % of an exposed group, `D_y = 31.8\,A(8)^{-1.06}` (Eq. (C.1)): 253 + 254 + ```python 255 + import phonometry as ph 256 + 257 + print(round(ph.hav_vwf_lifetime_years(7.0), 1)) # 4.0 years (Table C.1) 258 + ``` 259 + 260 + The standards deliberately define no safe limit — `A(8)` and the directive's 261 + action and limit values are the basis for any exposure criterion. For whole-body 262 + exposure, `energy_equivalent_acceleration` gives the ISO 2631-1 Eq. (B.3) 263 + energy-equivalent magnitude across periods of different magnitude and duration. 264 + 265 + --- 266 + 267 + **Standards.** ISO 8041-1:2017 (weighting definitions and tolerances); 268 + ISO 2631-1:1997 (whole-body evaluation); ISO 2631-2:2003 (buildings, `Wm`); 269 + ISO 2631-4:2001 (rail ride comfort, `Wb`); ISO 5349-1:2001 and ISO 5349-2:2001 270 + (hand-transmitted vibration); Directive 2002/44/EC (exposure action and limit 271 + values).
+70
scripts/conformance_report.py
··· 1267 1267 1268 1268 1269 1269 # =========================================================================== 1270 + # Human vibration (ISO 8041-1 / ISO 2631 / ISO 5349 / Directive 2002/44/EC) 1271 + # =========================================================================== 1272 + _HUMAN_VIB = "Human vibration (ISO 8041 / 2631 / 5349)" 1273 + 1274 + 1275 + def _true_centre(n: int) -> float: 1276 + """True IEC 61260 one-third-octave centre ``10^(n/10)`` Hz.""" 1277 + return float(10.0 ** (n / 10.0)) 1278 + 1279 + 1280 + @register(_HUMAN_VIB, "ISO 8041-1:2017 Table B.8", "Wk design-goal factor at 6,31 Hz") 1281 + def _chk_iso8041_wk_annex_b() -> Outcome: 1282 + factor = float(ph.weighting_factors("Wk", _true_centre(8))[0]) 1283 + return numeric(ref.ISO8041_1_WK_FACTOR_6P31HZ, factor, 1e-3, rel=True, places=4) 1284 + 1285 + 1286 + @register(_HUMAN_VIB, "ISO 8041-1:2017 Table B.9", "Wm design-goal factor at 1,585 Hz") 1287 + def _chk_iso8041_wm_annex_b() -> Outcome: 1288 + factor = float(ph.weighting_factors("Wm", _true_centre(2))[0]) 1289 + return numeric(ref.ISO8041_1_WM_FACTOR_1P585HZ, factor, 1e-3, rel=True, places=4) 1290 + 1291 + 1292 + @register(_HUMAN_VIB, "ISO 8041-1:2017 Table 1", "Wh factor at the 500 rad/s reference") 1293 + def _chk_iso8041_wh_reference() -> Outcome: 1294 + factor = float(ph.weighting_factors("Wh", ref.ISO8041_1_WH_REF_FREQ_HZ)[0]) 1295 + return numeric(ref.ISO8041_1_WH_REF_FACTOR, factor, 1.5e-3, rel=True, places=4) 1296 + 1297 + 1298 + @register(_HUMAN_VIB, "ISO 5349-2:2001 Example E.2.1", "Single-tool daily exposure A(8)") 1299 + def _chk_iso5349_e21() -> Outcome: 1300 + a8 = ph.daily_exposure(7.4, 2.5 * 3600.0) 1301 + return numeric(ref.ISO5349_2_E21_A8, a8, 0.05, unit="m/s^2", places=2) 1302 + 1303 + 1304 + @register(_HUMAN_VIB, "ISO 5349-2:2001 Example E.3", "Forestry three-task A(8)") 1305 + def _chk_iso5349_e3() -> Outcome: 1306 + a8 = ph.hav_daily_exposure( 1307 + [4.6, 6.0, 3.6], [2 * 3600.0, 1 * 3600.0, 2 * 3600.0] 1308 + ) 1309 + return numeric(ref.ISO5349_2_E3_A8, a8, 0.05, unit="m/s^2", places=2) 1310 + 1311 + 1312 + @register(_HUMAN_VIB, "ISO 5349-1:2001 Eq. (C.1)", "VWF 10 % lifetime Dy at A(8)=7") 1313 + def _chk_iso5349_vwf() -> Outcome: 1314 + dy = ph.hav_vwf_lifetime_years(ref.ISO5349_1_VWF_A8) 1315 + return numeric(ref.ISO5349_1_VWF_DY_YEARS, dy, 0.1, unit="yr", places=2) 1316 + 1317 + 1318 + @register(_HUMAN_VIB, "Directive 2002/44/EC Art. 3", "HAV/WBV action & limit values") 1319 + def _chk_directive_2002_44() -> Outcome: 1320 + hav = ph.exposure_assessment(1.0, kind="hav") 1321 + wbv = ph.exposure_assessment(0.1, kind="wbv") 1322 + ok = ( 1323 + hav.action_value == ref.DIRECTIVE_2002_44_HAV_EAV 1324 + and hav.limit_value == ref.DIRECTIVE_2002_44_HAV_ELV 1325 + and wbv.action_value == ref.DIRECTIVE_2002_44_WBV_EAV 1326 + and wbv.limit_value == ref.DIRECTIVE_2002_44_WBV_ELV 1327 + ) 1328 + exp = ( 1329 + f"HAV {ref.DIRECTIVE_2002_44_HAV_EAV}/{ref.DIRECTIVE_2002_44_HAV_ELV}, " 1330 + f"WBV {ref.DIRECTIVE_2002_44_WBV_EAV}/{ref.DIRECTIVE_2002_44_WBV_ELV} m/s^2" 1331 + ) 1332 + got = ( 1333 + f"HAV {hav.action_value}/{hav.limit_value}, " 1334 + f"WBV {wbv.action_value}/{wbv.limit_value} m/s^2" 1335 + ) 1336 + return Outcome(expected=exp, computed=got, delta="0", passed=ok) 1337 + 1338 + 1339 + # =========================================================================== 1270 1340 # Markdown rendering 1271 1341 # =========================================================================== 1272 1342 def _snap(value: float, eps: float = 5e-4) -> float:
+65
scripts/generate_diagrams.py
··· 44 44 _ES: dict[str, str] = { 45 45 "Calibration chain — from calibrator to physical units": 46 46 "Cadena de calibración — del calibrador a unidades físicas", 47 + # Human vibration (ISO 2631-1 / ISO 8041-1 / 2002-44-EC) 48 + "Whole-body vibration measurement chain (ISO 2631-1 / ISO 8041-1)": 49 + "Cadena de medición de vibración de cuerpo entero (ISO 2631-1 / ISO 8041-1)", 50 + "vibration input": "entrada de vibración", 51 + "Seat/body interface": "Interfaz asiento/cuerpo", 52 + "Triaxial accelerometer": "Acelerómetro triaxial", 53 + "Band limiting + Wk / Wd": "Limitación de banda + Wk / Wd", 54 + "weighting (ISO 8041-1)": "ponderación (ISO 8041-1)", 55 + "Weighted r.m.s. a_w & VDV": "a_w eficaz ponderada y VDV", 56 + "(ISO 2631-1)": "(ISO 2631-1)", 57 + "assessed vs EAV / ELV (Directive 2002/44/EC)": 58 + "evaluada frente a EAV / ELV (Directiva 2002/44/CE)", 47 59 "Sound calibrator": "Calibrador acústico", 48 60 "Microphone +": "Micrófono +", 49 61 "preamplifier": "preamplificador", ··· 1891 1903 s.text(450, y, txt, 19 if bold else 18, col, bold=bold) 1892 1904 1893 1905 1906 + def _d_human_vibration(s: SVG, th: Theme) -> None: 1907 + """Whole-body vibration measurement chain (ISO 2631-1 / ISO 8041-1).""" 1908 + gy = 510.0 1909 + # --- Left: a seated person on a vibrating seat, triaxial accelerometer --- 1910 + s.ground(gy, 40, 350) 1911 + # Seat: cushion, backrest and support leg. 1912 + s.rect(118, 424, 132, 18, th.panel, th.fg, rx=4, sw=2) # cushion 1913 + s.rect(118, 336, 16, 90, th.panel, th.fg, rx=3, sw=2) # backrest 1914 + s.line(184, 442, 184, gy, th.fg, 2.4) # pedestal 1915 + # A wavy "vibration" arrow rising into the seat base. 1916 + s.arrow(184, gy - 4, 184, 452, th.secondary, 2.4) 1917 + s.text(184, gy - 12, "vibration input", 17, th.secondary, "middle", italic=True) 1918 + s.person(178, gy, 176, seated=True) 1919 + # Triaxial accelerometer at the seat/body interface with its x, y, z axes. 1920 + ox, oy = 176.0, 420.0 1921 + s.rect(ox - 9, oy - 8, 18, 16, th.secondary, th.fg, rx=2, sw=1.5) 1922 + s.arrow(ox, oy - 8, ox, oy - 58, th.accent, 2.0) # z (vertical) 1923 + s.text(ox + 8, oy - 54, "z", 18, th.accent, "start", bold=True) 1924 + s.arrow(ox + 9, oy, ox + 62, oy, th.accent, 2.0) # x (fore-aft) 1925 + s.text(ox + 66, oy + 5, "x", 18, th.accent, "start", bold=True) 1926 + s.arrow(ox - 7, oy + 6, ox - 44, oy + 34, th.accent, 2.0) # y (lateral) 1927 + s.text(ox - 52, oy + 44, "y", 18, th.accent, "end", bold=True) 1928 + s.text(150, gy + 34, "Seat/body interface", 18, th.fg, "middle") 1929 + 1930 + # --- Right: the vertical signal-processing chain --- 1931 + cx, bw, bh = 650.0, 320.0, 72.0 1932 + x0 = cx - bw / 2 1933 + chain = [ 1934 + (96.0, "Triaxial accelerometer", "a_x , a_y , a_z (m/s²)"), 1935 + (206.0, "Band limiting + Wk / Wd", "weighting (ISO 8041-1)"), 1936 + (316.0, "Weighted r.m.s. a_w & VDV", "(ISO 2631-1)"), 1937 + ] 1938 + for by, l1, l2 in chain: 1939 + s.rect(x0, by, bw, bh, th.panel, th.primary, rx=12, sw=2) 1940 + s.text(cx, by + 31, l1, 21, th.fg, "middle", bold=True) 1941 + s.text(cx, by + 56, l2, 18, th.muted, "middle") 1942 + s.arrow(cx, 168, cx, 206, th.fg, 2.0) 1943 + s.arrow(cx, 278, cx, 316, th.fg, 2.0) 1944 + # Feed the setup into the chain. 1945 + s.arrow(252, oy, x0 - 6, 132, th.fg, 2.0) 1946 + 1947 + # --- Bottom: vector sum, daily exposure and the Directive assessment --- 1948 + s.arrow(cx, 388, cx, 424, th.fg, 2.0) 1949 + s.rect(400, 424, 470, 78, "none", th.secondary, rx=12, sw=2, dash="6,5") 1950 + s.text(635, 452, "a_v = √(Σ k_j² a_wj²) → A(8) = a_v·√(T/T₀)", 1951 + 20, th.fg, "middle", bold=True) 1952 + s.text(635, 480, "assessed vs EAV / ELV (Directive 2002/44/EC)", 1953 + 18, th.secondary, "middle") 1954 + 1955 + 1894 1956 DIAGRAMS = { 1895 1957 "diagram_calibration_setup": (_d1, "Calibration chain — from calibrator to physical units", 560), 1896 1958 "diagram_env_measurement": (_d2, "Environmental noise measurement positions (ISO 1996-2)", 560), ··· 1937 1999 "diagram_intensity_scan": ( 1938 2000 _d_intensity_scan, 1939 2001 "Precision sound intensity scanning (ISO 9614-3)", 600), 2002 + "diagram_human_vibration": ( 2003 + _d_human_vibration, 2004 + "Whole-body vibration measurement chain (ISO 2631-1 / ISO 8041-1)", 580), 1940 2005 } 1941 2006 1942 2007
+140
scripts/generate_graphs.py
··· 122 122 "In-situ road-surface absorption (ISO 13472-1)": 123 123 "Absorción in situ de pavimentos (ISO 13472-1)", 124 124 "Absorption coefficient alpha": "Coeficiente de absorción alpha", 125 + # Human vibration (ISO 8041-1 / ISO 2631 / ISO 5349 / 2002/44/EC) 126 + "Whole-body vertical weighting Wk (ISO 8041-1)": 127 + "Ponderación vertical de cuerpo entero Wk (ISO 8041-1)", 128 + "Weighting factor [dB]": "Factor de ponderación [dB]", 129 + "r.m.s. acceleration [m/s$^2$]": "Aceleración eficaz [m/s$^2$]", 130 + "Unweighted $a_i$": "Sin ponderar $a_i$", 131 + "Weighted $W_i\\,a_i$ (Wk)": "Ponderada $W_i\\,a_i$ (Wk)", 132 + "Daily exposure A(8) [m/s$^2$]": "Exposición diaria A(8) [m/s$^2$]", 133 + "brush-saw": "desbrozadora", 134 + "felling": "tala", 135 + "stripping": "descortezado", 125 136 # Precision sound power (ISO 3745 / ISO 9614-3) 126 137 "Sound power level LW [dB]": "Nivel de potencia sonora LW [dB]", 127 138 "Non-applicable band": "Banda no aplicable", ··· 374 385 r"Potencia sonora de precisión (ISO 3745) LWA = \1 dB(A)"), 375 386 (r"^Precision intensity scanning \(ISO 9614-3\) LWA = (.+) dB\(A\)$", 376 387 r"Barrido de intensidad de precisión (ISO 9614-3) LWA = \1 dB(A)"), 388 + # Human-vibration dynamic titles (numeric a_w / A(8)) 389 + (r"^Weighted seat acceleration \(ISO 2631-1\) (.+)$", 390 + r"Aceleración ponderada del asiento (ISO 2631-1) \1"), 391 + (r"^Hand-arm daily exposure \(ISO 5349 / 2002-44-EC\) (.+)$", 392 + r"Exposición diaria mano-brazo (ISO 5349 / 2002-44-EC) \1"), 377 393 ] 378 394 379 395 ··· 3069 3085 plt.close() 3070 3086 3071 3087 3088 + def generate_vibration_weighting(output_dir: str) -> None: 3089 + """ISO 8041-1: the whole-body vertical weighting Wk over its band.""" 3090 + print("Generating vibration_weighting.png...") 3091 + from phonometry import frequency_weighting 3092 + 3093 + # A user evaluates the principal ISO 2631-1 weighting Wk on a fine 3094 + # frequency grid across the whole-body band (0,4-100 Hz). The result is the 3095 + # ISO 8041-1 cascade H(f): a gentle +0,5 dB peak near 6 Hz, a band-limiting 3096 + # roll-off below 0,4 Hz and above ~16 Hz. 3097 + freqs = np.geomspace(0.4, 100.0, 240) 3098 + result = frequency_weighting("Wk", freqs) 3099 + 3100 + fig, ax = plt.subplots(figsize=(10, 6.3)) 3101 + ax.semilogx(result.frequencies, result.magnitude_db, color=COLOR_PRIMARY, 3102 + linewidth=1.9, zorder=3) 3103 + ax.axhline(0.0, color=COLOR_FG, linewidth=0.8, alpha=0.4, zorder=1) 3104 + ax.set_title("Whole-body vertical weighting Wk (ISO 8041-1)", 3105 + fontweight="bold", pad=12) 3106 + ax.set_xlabel(LABEL_FREQ_HZ) 3107 + ax.set_ylabel("Weighting factor [dB]") 3108 + ax.set_xlim(0.4, 100.0) 3109 + ax.set_ylim(-40.0, 5.0) 3110 + from matplotlib.ticker import NullFormatter 3111 + ax.set_xticks([0.5, 1, 2, 5, 10, 20, 50, 100]) 3112 + # Explicit string labels install a FixedFormatter so the Spanish pass can 3113 + # apply the decimal comma (a log-axis ScalarFormatter would not be caught). 3114 + ax.set_xticklabels(["0.5", "1", "2", "5", "10", "20", "50", "100"]) 3115 + ax.xaxis.set_minor_formatter(NullFormatter()) 3116 + ax.grid(which="major", color=COLOR_GRID, linestyle="-", alpha=0.5) 3117 + ax.set_axisbelow(True) 3118 + plt.tight_layout() 3119 + plt.savefig(themed_path(output_dir, "vibration_weighting.png")) 3120 + plt.close() 3121 + 3122 + 3123 + def generate_weighted_acceleration(output_dir: str) -> None: 3124 + """ISO 2631-1: measured seat spectrum weighted to a_w (Eq. (9)).""" 3125 + print("Generating weighted_acceleration.png...") 3126 + from phonometry import weighted_acceleration 3127 + 3128 + # A measured vertical seat-pan acceleration spectrum (r.m.s. per one-third 3129 + # octave, m/s^2) from a vehicle seat: energy concentrated in the 2-8 Hz 3130 + # whole-body range. Weighting it with Wk gives the health-relevant a_w. 3131 + freqs = np.array([1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 3132 + 10.0, 12.5, 16.0, 20.0, 25.0, 31.5, 40.0, 63.0, 80.0]) 3133 + accel = np.array([0.18, 0.24, 0.33, 0.46, 0.52, 0.55, 0.48, 0.39, 0.31, 3134 + 0.26, 0.21, 0.17, 0.13, 0.10, 0.078, 0.060, 0.045, 3135 + 0.028, 0.020]) 3136 + result = weighted_acceleration(accel, freqs, "Wk") 3137 + 3138 + positions = np.arange(freqs.size, dtype=float) 3139 + width = 0.4 3140 + fig, ax = plt.subplots(figsize=(10.5, 6.3)) 3141 + ax.bar(positions - width / 2, result.band_accelerations, width, 3142 + color=COLOR_GRID, edgecolor=COLOR_FG, linewidth=0.5, 3143 + label="Unweighted $a_i$", zorder=2) 3144 + ax.bar(positions + width / 2, result.weighted, width, color=COLOR_PRIMARY, 3145 + edgecolor=COLOR_FG, linewidth=0.5, label="Weighted $W_i\\,a_i$ (Wk)", 3146 + zorder=3) 3147 + ax.set_xticks(positions) 3148 + ax.set_xticklabels([f"{f:g}" for f in freqs], rotation=45, ha="right") 3149 + ax.set_title( 3150 + f"Weighted seat acceleration (ISO 2631-1) $a_w$ = {result.overall:.3f} " 3151 + "m/s$^2$", fontweight="bold", pad=12) 3152 + ax.set_xlabel(LABEL_FREQ_HZ) 3153 + ax.set_ylabel("r.m.s. acceleration [m/s$^2$]") 3154 + ax.legend(loc="upper right", fontsize=9) 3155 + ax.grid(axis="y", color=COLOR_GRID, linestyle="--", alpha=0.5, zorder=0) 3156 + ax.set_axisbelow(True) 3157 + plt.tight_layout() 3158 + plt.savefig(themed_path(output_dir, "weighted_acceleration.png")) 3159 + plt.close() 3160 + 3161 + 3162 + def generate_daily_vibration_exposure(output_dir: str) -> None: 3163 + """ISO 5349 + Directive 2002/44/EC: A(8) vs the EAV/ELV thresholds.""" 3164 + print("Generating daily_vibration_exposure.png...") 3165 + from phonometry import daily_vibration_exposure 3166 + 3167 + # A forestry worker's day across three chain-saw tasks (the ISO 5349-2 3168 + # Annex E.3 worked example): each task's a_hv and duration give a partial 3169 + # exposure A_i(8); they combine to A(8) = 3,6 m/s^2, assessed against the 3170 + # hand-arm action (2,5) and limit (5,0) values of Directive 2002/44/EC. 3171 + result = daily_vibration_exposure( 3172 + [4.6, 6.0, 3.6], 3173 + [2 * 3600.0, 1 * 3600.0, 2 * 3600.0], 3174 + kind="hav", 3175 + labels=["brush-saw", "felling", "stripping"], 3176 + ) 3177 + 3178 + labels = [*result.labels, "A(8)"] 3179 + values = [*result.partials.tolist(), result.a8] 3180 + positions = np.arange(len(values), dtype=float) 3181 + colors = [COLOR_GRID] * result.partials.size + [COLOR_PRIMARY] 3182 + fig, ax = plt.subplots(figsize=(9.5, 6.3)) 3183 + ax.bar(positions, values, width=0.62, color=colors, edgecolor=COLOR_FG, 3184 + linewidth=0.6, zorder=3) 3185 + eav = result.assessment.action_value 3186 + elv = result.assessment.limit_value 3187 + ax.axhline(eav, color=COLOR_TERTIARY, linestyle="--", linewidth=1.6, 3188 + label=f"EAV = {eav:g} m/s$^2$", zorder=2) 3189 + ax.axhline(elv, color=COLOR_SECONDARY, linestyle="--", linewidth=1.6, 3190 + label=f"ELV = {elv:g} m/s$^2$", zorder=2) 3191 + ax.set_xticks(positions) 3192 + ax.set_xticklabels(labels, rotation=30, ha="right") 3193 + ax.set_ylabel("Daily exposure A(8) [m/s$^2$]") 3194 + ax.set_ylim(0.0, elv * 1.2) 3195 + ax.set_title( 3196 + f"Hand-arm daily exposure (ISO 5349 / 2002-44-EC) A(8) = " 3197 + f"{result.a8:.2f} m/s$^2$", fontweight="bold", pad=12) 3198 + ax.legend(loc="upper left", fontsize=9) 3199 + ax.grid(axis="y", color=COLOR_GRID, linestyle="--", alpha=0.5, zorder=0) 3200 + ax.set_axisbelow(True) 3201 + plt.tight_layout() 3202 + plt.savefig(themed_path(output_dir, "daily_vibration_exposure.png")) 3203 + plt.close() 3204 + 3205 + 3072 3206 def generate_all(img_dir: str) -> None: 3073 3207 """Generate every documentation figure for the currently active theme.""" 3074 3208 generate_filter_type_comparison(img_dir) ··· 3131 3265 generate_insitu_absorption(img_dir) 3132 3266 generate_precision_anechoic_power(img_dir) 3133 3267 generate_intensity_scan_power(img_dir) 3268 + 3269 + # Human vibration (ISO 8041-1, ISO 2631-1/-2/-4, ISO 5349-1/-2, 3270 + # Directive 2002/44/EC): frequency weighting, weighted a_w, daily A(8) 3271 + generate_vibration_weighting(img_dir) 3272 + generate_weighted_acceleration(img_dir) 3273 + generate_daily_vibration_exposure(img_dir) 3134 3274 3135 3275 # Psychoacoustics / open-plan plots (sharpness weighting, spatial decay) 3136 3276 generate_sharpness_weighting(img_dir)
+67
src/phonometry/__init__.py
··· 153 153 spot_microphone_spacing_bounds, 154 154 spot_tube_upper_frequency, 155 155 ) 156 + from .human_vibration import ( 157 + HAV_EAV_A8, 158 + HAV_ELV_A8, 159 + REFERENCE_ACCELERATION, 160 + REFERENCE_DURATION_S, 161 + WBV_EAV_A8, 162 + WBV_EAV_VDV, 163 + WBV_ELV_A8, 164 + WBV_ELV_VDV, 165 + WEIGHTING_NAMES, 166 + DailyVibrationExposure, 167 + ExposureAssessment, 168 + HumanVibrationWarning, 169 + WeightedSpectrum, 170 + WeightingResponse, 171 + apply_weighting, 172 + combine_partial_exposures, 173 + crest_factor, 174 + daily_exposure, 175 + daily_vibration_exposure, 176 + energy_equivalent_acceleration, 177 + exposure_assessment, 178 + frequency_weighting, 179 + hav_daily_exposure, 180 + hav_vwf_lifetime_years, 181 + motion_sickness_dose_value, 182 + mtvv, 183 + partial_exposure, 184 + running_rms, 185 + vibration_dose_value, 186 + vibration_total_value, 187 + weighted_acceleration, 188 + weighting_factors, 189 + ) 156 190 from .room_acoustics import ( 157 191 DecayCurve, 158 192 RoomAcousticsResult, ··· 499 533 "SPOT_FREQUENCY_RANGE", 500 534 "SPOT_NARROW_BAND_RANGE", 501 535 "RoadAbsorptionWarning", 536 + # ISO 8041-1 / ISO 2631 / ISO 5349 / Directive 2002/44/EC human vibration 537 + "frequency_weighting", 538 + "weighting_factors", 539 + "apply_weighting", 540 + "WeightingResponse", 541 + "WEIGHTING_NAMES", 542 + "weighted_acceleration", 543 + "WeightedSpectrum", 544 + "running_rms", 545 + "mtvv", 546 + "vibration_dose_value", 547 + "motion_sickness_dose_value", 548 + "crest_factor", 549 + "vibration_total_value", 550 + "daily_exposure", 551 + "partial_exposure", 552 + "combine_partial_exposures", 553 + "hav_daily_exposure", 554 + "energy_equivalent_acceleration", 555 + "hav_vwf_lifetime_years", 556 + "exposure_assessment", 557 + "ExposureAssessment", 558 + "daily_vibration_exposure", 559 + "DailyVibrationExposure", 560 + "HumanVibrationWarning", 561 + "REFERENCE_ACCELERATION", 562 + "REFERENCE_DURATION_S", 563 + "HAV_EAV_A8", 564 + "HAV_ELV_A8", 565 + "WBV_EAV_A8", 566 + "WBV_ELV_A8", 567 + "WBV_EAV_VDV", 568 + "WBV_ELV_VDV", 502 569 "airborne_insulation", 503 570 "AirborneInsulationResult", 504 571 "impact_insulation",
+117
src/phonometry/_plotting.py
··· 1047 1047 ax.set_title("In-situ road-surface absorption (ISO 13472-1)") 1048 1048 ax.grid(True, axis="y", alpha=0.3) 1049 1049 return ax 1050 + 1051 + 1052 + # --------------------------------------------------------------------------- 1053 + # Human vibration (ISO 8041-1 / ISO 2631 / ISO 5349 / Directive 2002/44/EC) 1054 + # --------------------------------------------------------------------------- 1055 + 1056 + 1057 + def plot_vibration_weighting( 1058 + result: Any, ax: Axes | None = None, **kwargs: Any 1059 + ) -> Axes: 1060 + """Frequency-weighting factor (dB) versus frequency (ISO 8041-1). 1061 + 1062 + :param result: A 1063 + :class:`~phonometry.human_vibration.WeightingResponse` exposing 1064 + ``name``, ``frequencies`` and ``magnitude_db``. 1065 + :param ax: Existing axes, or ``None`` to create a figure. 1066 + :return: The axes. 1067 + """ 1068 + ax = ax if ax is not None else _new_axes() 1069 + freqs = np.asarray(result.frequencies, dtype=np.float64) 1070 + mag_db = np.asarray(result.magnitude_db, dtype=np.float64) 1071 + kwargs.setdefault("color", "#1f77b4") 1072 + ax.semilogx(freqs, mag_db, **kwargs) 1073 + ax.set_xlabel("Frequency [Hz]") 1074 + ax.set_ylabel("Weighting factor [dB]") 1075 + ax.set_title(f"Frequency weighting {result.name} (ISO 8041-1)") 1076 + ax.grid(True, which="both", alpha=0.3) 1077 + return ax 1078 + 1079 + 1080 + def plot_weighted_spectrum( 1081 + result: Any, ax: Axes | None = None, **kwargs: Any 1082 + ) -> Axes: 1083 + """Unweighted vs weighted one-third-octave acceleration spectrum. 1084 + 1085 + Draws the measured band accelerations and, overlaid, the weighted band 1086 + contributions ``W_i*a_i``; the overall ``a_w`` is annotated in the title. 1087 + 1088 + :param result: A 1089 + :class:`~phonometry.human_vibration.WeightedSpectrum` exposing 1090 + ``frequencies``, ``band_accelerations``, ``weighted``, ``overall`` and 1091 + ``weighting_name``. 1092 + :param ax: Existing axes, or ``None`` to create a figure. 1093 + :return: The axes. 1094 + """ 1095 + ax = ax if ax is not None else _new_axes() 1096 + freqs = np.asarray(result.frequencies, dtype=np.float64) 1097 + raw = np.asarray(result.band_accelerations, dtype=np.float64) 1098 + weighted = np.asarray(result.weighted, dtype=np.float64) 1099 + positions = np.arange(freqs.size, dtype=np.float64) 1100 + width = 0.4 1101 + # The weighted bars are the primary artist; forward user kwargs there. 1102 + kwargs.setdefault("color", "#1f77b4") 1103 + ax.bar( 1104 + positions - width / 2, raw, width, color="#bbbbbb", label="Unweighted $a_i$" 1105 + ) 1106 + ax.bar( 1107 + positions + width / 2, 1108 + weighted, 1109 + width, 1110 + label=f"Weighted $W_i a_i$ ({result.weighting_name})", 1111 + **kwargs, 1112 + ) 1113 + ax.set_xticks(positions) 1114 + ax.set_xticklabels([_format_freq(f) for f in freqs], rotation=45, ha="right") 1115 + ax.set_xlabel("Frequency [Hz]") 1116 + ax.set_ylabel(r"r.m.s. acceleration [m/s$^2$]") 1117 + ax.set_title( 1118 + f"Weighted acceleration spectrum ($a_w$ = {float(result.overall):.3f} " 1119 + r"m/s$^2$)" 1120 + ) 1121 + ax.legend(loc="best", fontsize="small") 1122 + ax.grid(True, axis="y", alpha=0.3) 1123 + return ax 1124 + 1125 + 1126 + def plot_daily_exposure( 1127 + result: Any, ax: Axes | None = None, **kwargs: Any 1128 + ) -> Axes: 1129 + """Partial daily exposures against the EAV / ELV (Directive 2002/44/EC). 1130 + 1131 + Draws one bar per operation (its partial exposure ``A_i(8)``), a combined 1132 + ``A(8)`` bar, and the exposure action and limit value as horizontal lines. 1133 + 1134 + :param result: A 1135 + :class:`~phonometry.human_vibration.DailyVibrationExposure` exposing 1136 + ``labels``, ``partials``, ``a8`` and ``assessment``. 1137 + :param ax: Existing axes, or ``None`` to create a figure. 1138 + :return: The axes. 1139 + """ 1140 + ax = ax if ax is not None else _new_axes() 1141 + partials = np.asarray(result.partials, dtype=np.float64) 1142 + labels = [*result.labels, "A(8)"] 1143 + values = [*partials.tolist(), float(result.a8)] 1144 + positions = np.arange(len(values), dtype=np.float64) 1145 + colors = ["#bbbbbb"] * partials.size + ["#1f77b4"] 1146 + kwargs.setdefault("color", colors) 1147 + ax.bar(positions, values, **kwargs) 1148 + ax.set_xticks(positions) 1149 + ax.set_xticklabels(labels, rotation=45, ha="right") 1150 + ax.set_ylabel(r"Vibration exposure A(8) [m/s$^2$]") 1151 + 1152 + assessment = result.assessment 1153 + eav = float(assessment.action_value) 1154 + elv = float(assessment.limit_value) 1155 + ax.axhline(eav, color="#ff7f0e", ls="--", label=f"EAV = {eav:g}") 1156 + ax.axhline(elv, color="#d62728", ls="--", label=f"ELV = {elv:g}") 1157 + top = max(elv, float(np.max(values))) * 1.15 1158 + ax.set_ylim(0.0, top) 1159 + kind = str(assessment.kind).upper() 1160 + ax.set_title( 1161 + f"Daily {kind} exposure (A(8) = {float(result.a8):.2f} " 1162 + rf"m/s$^2$, {assessment.zone})" 1163 + ) 1164 + ax.legend(loc="best", fontsize="small") 1165 + ax.grid(True, axis="y", alpha=0.3) 1166 + return ax
+851
src/phonometry/human_vibration.py
··· 1 + # Copyright (c) 2026. Jose M. Requena-Plens 2 + """ 3 + Human exposure to whole-body and hand-transmitted vibration. 4 + 5 + The measurement chain of the ISO human-vibration family is implemented from the 6 + standards' own analog definitions, clean-room: 7 + 8 + * **ISO 8041-1:2017** - the authoritative *master* definition of every 9 + frequency weighting. A single cascade ``H(s) = Hh(s)*Hl(s)*Ht(s)*Hs(s)`` 10 + (Formula (5)) of second-order band-limiting Butterworth sections 11 + (Formulae (1)/(2)), an acceleration-velocity transition (Formula (3)) and an 12 + upward step (Formula (4)) realises all nine weightings from the one Table 3 13 + parameter set: ``Wb, Wc, Wd, We, Wf, Wh, Wj, Wk, Wm``. The tabulated 14 + design-goal factors of Annex B (Tables B.1-B.9) are reproduced to their 15 + four-significant-figure precision. 16 + 17 + * **ISO 2631-1:1997** - whole-body vibration: the weighted r.m.s. acceleration 18 + ``a_w`` (Eq. (1)/(9)), the vibration total value ``a_v`` with axis 19 + multiplying factors ``k`` (Eq. (10)), the running r.m.s. and maximum 20 + transient vibration value ``MTVV`` (Eqs. (2)-(4)), the vibration dose value 21 + ``VDV`` (Eq. (5)), the crest factor (6.2.1) and the energy-equivalent 22 + magnitude relations of Annex B (Eqs. (B.1)-(B.3)). 23 + 24 + * **ISO 2631-2:2003** - whole-body vibration in buildings: the direction- 25 + independent weighting ``Wm`` (Annex A). 26 + 27 + * **ISO 2631-4:2001** - ride comfort in fixed-guideway (rail) transport: the 28 + vertical weighting ``Wb`` (Annex A) with the ``Wk`` axis multiplying 29 + factors. ``Wb`` is realised here from its exact ISO 8041-1 Table 3 30 + parameters (of which ISO 2631-4 Table A.1 is the two-decimal rounding). 31 + 32 + * **ISO 5349-1:2001 / ISO 5349-2:2001** - hand-transmitted vibration: the 33 + vibration total value ``a_hv`` (5349-1 Eq. (1)), the daily exposure ``A(8)`` 34 + for single and multiple operations (5349-1 Eqs. (2)/(3); 5349-2 Eqs. (1)-(3)) 35 + and the vibration-white-finger dose relation of 5349-1 Annex C (Eq. (C.1)). 36 + 37 + * **Directive 2002/44/EC** - the daily exposure action and limit values 38 + (Article 3) that ISO does not fix: hand-arm ``A(8)`` EAV ``2,5`` / 39 + ELV ``5`` m/s2; whole-body ``A(8)`` EAV ``0,5`` / ELV ``1,15`` m/s2 (or VDV 40 + EAV ``9,1`` / ELV ``21`` m/s^1,75). 41 + 42 + The band (spectrum) method and the exposure arithmetic carry the standards' 43 + worked-example oracles; the time-domain metrics operate on a weighted 44 + acceleration signal, which :func:`apply_weighting` produces from a raw record 45 + by applying the exact analog response of ISO 8041-1 in the frequency domain. 46 + """ 47 + 48 + from __future__ import annotations 49 + 50 + import math 51 + import warnings 52 + from dataclasses import dataclass 53 + from typing import TYPE_CHECKING, Any 54 + 55 + import numpy as np 56 + from numpy.typing import ArrayLike, NDArray 57 + from scipy import signal as sig 58 + 59 + if TYPE_CHECKING: # pragma: no cover - typing only 60 + from matplotlib.axes import Axes 61 + 62 + Real = NDArray[np.float64] 63 + Complex = NDArray[np.complex128] 64 + 65 + __all__ = [ 66 + "HAV_ELV_A8", 67 + "HAV_EAV_A8", 68 + "REFERENCE_ACCELERATION", 69 + "REFERENCE_DURATION_S", 70 + "WBV_EAV_A8", 71 + "WBV_EAV_VDV", 72 + "WBV_ELV_A8", 73 + "WBV_ELV_VDV", 74 + "WEIGHTING_NAMES", 75 + "DailyVibrationExposure", 76 + "ExposureAssessment", 77 + "HumanVibrationWarning", 78 + "WeightedSpectrum", 79 + "WeightingResponse", 80 + "apply_weighting", 81 + "combine_partial_exposures", 82 + "crest_factor", 83 + "daily_exposure", 84 + "daily_vibration_exposure", 85 + "energy_equivalent_acceleration", 86 + "exposure_assessment", 87 + "frequency_weighting", 88 + "hav_daily_exposure", 89 + "hav_vwf_lifetime_years", 90 + "motion_sickness_dose_value", 91 + "mtvv", 92 + "partial_exposure", 93 + "running_rms", 94 + "vibration_dose_value", 95 + "vibration_total_value", 96 + "weighted_acceleration", 97 + "weighting_factors", 98 + ] 99 + 100 + # --------------------------------------------------------------------------- 101 + # Reference constants. 102 + # --------------------------------------------------------------------------- 103 + #: Reference acceleration ``a0 = 10^-6 m/s2`` for vibration levels 104 + #: (ISO 8041-1:2017, 3.1.2.2, after ISO 1683). 105 + REFERENCE_ACCELERATION = 1e-6 106 + 107 + #: Reference duration ``T0 = 8 h = 28 800 s`` of the daily exposure ``A(8)`` 108 + #: (ISO 5349-1:2001, 3.2; ISO 2631-1:1997, B.1). 109 + REFERENCE_DURATION_S = 28800.0 110 + 111 + #: Daily hand-arm exposure action value ``A(8) = 2,5 m/s2`` (Directive 112 + #: 2002/44/EC, Article 3(1)(a)). 113 + HAV_EAV_A8 = 2.5 114 + #: Daily hand-arm exposure limit value ``A(8) = 5 m/s2`` (Directive 115 + #: 2002/44/EC, Article 3(1)(b)). 116 + HAV_ELV_A8 = 5.0 117 + #: Daily whole-body exposure action value ``A(8) = 0,5 m/s2`` (Directive 118 + #: 2002/44/EC, Article 3(2)(a)). 119 + WBV_EAV_A8 = 0.5 120 + #: Daily whole-body exposure limit value ``A(8) = 1,15 m/s2`` (Directive 121 + #: 2002/44/EC, Article 3(2)(b)). 122 + WBV_ELV_A8 = 1.15 123 + #: Whole-body VDV exposure action value ``9,1 m/s^1,75`` (Directive 124 + #: 2002/44/EC, Article 3(2)(a), alternative dose metric). 125 + WBV_EAV_VDV = 9.1 126 + #: Whole-body VDV exposure limit value ``21 m/s^1,75`` (Directive 127 + #: 2002/44/EC, Article 3(2)(b), alternative dose metric). 128 + WBV_ELV_VDV = 21.0 129 + 130 + _Q_BUTTERWORTH = 1.0 / math.sqrt(2.0) 131 + 132 + 133 + class HumanVibrationWarning(UserWarning): 134 + """Advisory for out-of-range human-vibration measurement conditions.""" 135 + 136 + 137 + # --------------------------------------------------------------------------- 138 + # Frequency weightings (ISO 8041-1:2017, Table 3 + Formulae (1)-(5)). 139 + # --------------------------------------------------------------------------- 140 + @dataclass(frozen=True) 141 + class _WParams: 142 + """Master parameters of one ISO 8041-1 Table 3 weighting (Hz / gain).""" 143 + 144 + f1: float # band-limiting high-pass corner 145 + q1: float 146 + f2: float # band-limiting low-pass corner 147 + q2: float 148 + f3: float # a-v transition zero (``inf`` -> absent) 149 + f4: float # a-v transition pole (``inf`` -> absent) 150 + q4: float 151 + f5: float # upward-step zero (``inf`` -> step absent) 152 + q5: float 153 + f6: float # upward-step pole (``inf`` -> step absent) 154 + q6: float 155 + k: float # overall gain K 156 + 157 + 158 + _INF = math.inf 159 + 160 + #: ISO 8041-1:2017, Table 3 - the master parameters of the nine frequency 161 + #: weightings. ``Q1 = Q2 = 1/sqrt(2)`` (band-limiting Butterworth) throughout; 162 + #: ``inf`` corners collapse the corresponding stage to unity (Table 3 NOTEs). 163 + #: ``Wh`` uses the exact ISO 8041-1 values (Table 3, NOTE 2), of which the 164 + #: ISO 5349-1 Table A.1 figures are the five-significant-figure rounding. 165 + _WEIGHTINGS: dict[str, _WParams] = { 166 + # Whole-body, vertical z (ISO 2631-4); K = 1,024. 167 + "Wb": _WParams(0.4, _Q_BUTTERWORTH, 100.0, _Q_BUTTERWORTH, 168 + 16.0, 16.0, 0.55, 2.5, 0.9, 4.0, 0.95, 1.024), 169 + # Whole-body, horizontal x seat-back (ISO 2631-1). 170 + "Wc": _WParams(0.4, _Q_BUTTERWORTH, 100.0, _Q_BUTTERWORTH, 171 + 8.0, 8.0, 0.63, _INF, 1.0, _INF, 1.0, 1.0), 172 + # Whole-body, horizontal x/y (ISO 2631-1). 173 + "Wd": _WParams(0.4, _Q_BUTTERWORTH, 100.0, _Q_BUTTERWORTH, 174 + 2.0, 2.0, 0.63, _INF, 1.0, _INF, 1.0, 1.0), 175 + # Whole-body, rotational (ISO 2631-1). 176 + "We": _WParams(0.4, _Q_BUTTERWORTH, 100.0, _Q_BUTTERWORTH, 177 + 1.0, 1.0, 0.63, _INF, 1.0, _INF, 1.0, 1.0), 178 + # Whole-body, motion sickness, vertical z (ISO 2631-1). 179 + "Wf": _WParams(0.08, _Q_BUTTERWORTH, 0.63, _Q_BUTTERWORTH, 180 + _INF, 0.25, 0.86, 0.0625, 0.80, 0.10, 0.80, 1.0), 181 + # Hand-arm, all directions (ISO 5349-1); exact ISO 8041-1 corners. 182 + "Wh": _WParams(10.0**0.8, _Q_BUTTERWORTH, 10.0**3.1, _Q_BUTTERWORTH, 183 + 100.0 / (2.0 * math.pi), 100.0 / (2.0 * math.pi), 0.64, 184 + _INF, 1.0, _INF, 1.0, 1.0), 185 + # Whole-body, vertical head, recumbent x (ISO 2631-1). 186 + "Wj": _WParams(0.4, _Q_BUTTERWORTH, 100.0, _Q_BUTTERWORTH, 187 + _INF, _INF, 1.0, 3.75, 0.91, 5.32, 0.91, 1.0), 188 + # Whole-body, vertical z, the principal ISO 2631-1 weighting. 189 + "Wk": _WParams(0.4, _Q_BUTTERWORTH, 100.0, _Q_BUTTERWORTH, 190 + 12.5, 12.5, 0.63, 2.37, 0.91, 3.35, 0.91, 1.0), 191 + # Whole-body in buildings, all directions (ISO 2631-2). 192 + "Wm": _WParams(10.0**-0.1, _Q_BUTTERWORTH, 100.0, _Q_BUTTERWORTH, 193 + 1.0 / (0.028 * 2.0 * math.pi), 1.0 / (0.028 * 2.0 * math.pi), 194 + 0.5, _INF, 1.0, _INF, 1.0, 1.0), 195 + } 196 + 197 + #: The nine ISO 8041-1 frequency-weighting names, in Table 3 order. 198 + WEIGHTING_NAMES: tuple[str, ...] = tuple(_WEIGHTINGS) 199 + 200 + 201 + def _params(name: str) -> _WParams: 202 + """Return the Table 3 parameters for ``name`` or raise ``ValueError``.""" 203 + try: 204 + return _WEIGHTINGS[name] 205 + except KeyError: 206 + raise ValueError( 207 + f"Unknown weighting {name!r}; choose from {', '.join(WEIGHTING_NAMES)}." 208 + ) from None 209 + 210 + 211 + def _weighting_response(name: str, freq: Real) -> Complex: 212 + """Complex ``H(j*2*pi*f)`` of weighting ``name`` (ISO 8041-1, Formula (5)). 213 + 214 + The four cascaded stages of Formulae (1)-(4) are evaluated directly. A 215 + corner at infinity collapses its stage to unity, exactly as the Table 3 216 + NOTEs prescribe. ``f <= 0`` returns ``0`` (the high-pass blocks DC). 217 + """ 218 + p = _params(name) 219 + out = np.zeros(freq.shape, dtype=np.complex128) 220 + positive = freq > 0.0 221 + if not np.any(positive): 222 + return out 223 + s = 1j * 2.0 * math.pi * freq[positive].astype(np.float64) 224 + 225 + w1 = 2.0 * math.pi * p.f1 226 + w2 = 2.0 * math.pi * p.f2 227 + # Formula (1) high-pass and Formula (2) low-pass band limiting. 228 + hh = 1.0 / (1.0 + w1 / (p.q1 * s) + (w1 / s) ** 2) 229 + hl = 1.0 / (1.0 + s / (p.q2 * w2) + (s / w2) ** 2) 230 + 231 + # Formula (3) acceleration-velocity transition, gain K. 232 + ones = np.ones_like(s) 233 + num_t = ones if math.isinf(p.f3) else 1.0 + s / (2.0 * math.pi * p.f3) 234 + if math.isinf(p.f4): 235 + den_t = ones 236 + else: 237 + w4 = 2.0 * math.pi * p.f4 238 + den_t = 1.0 + s / (p.q4 * w4) + (s / w4) ** 2 239 + ht = p.k * num_t / den_t 240 + 241 + # Formula (4) upward step. 242 + if math.isinf(p.f5) or math.isinf(p.f6): 243 + hs = ones 244 + else: 245 + w5 = 2.0 * math.pi * p.f5 246 + w6 = 2.0 * math.pi * p.f6 247 + hs = (1.0 + s / (p.q5 * w5) + (s / w5) ** 2) / ( 248 + 1.0 + s / (p.q6 * w6) + (s / w6) ** 2 249 + ) * (w5 / w6) ** 2 250 + 251 + out[positive] = np.asarray(hh * hl * ht * hs, dtype=np.complex128) 252 + return out 253 + 254 + 255 + @dataclass(frozen=True) 256 + class WeightingResponse: 257 + """A frequency-weighting magnitude response (ISO 8041-1, Formula (5)). 258 + 259 + :ivar name: Weighting name (one of :data:`WEIGHTING_NAMES`). 260 + :ivar frequencies: Frequencies at which the response was evaluated, in Hz. 261 + :ivar response: Complex weighting ``H(j*2*pi*f)`` per frequency. 262 + :ivar magnitude: Weighting factor ``|H|`` per frequency. 263 + :ivar magnitude_db: ``20*log10(|H|)`` per frequency, in decibels. 264 + """ 265 + 266 + name: str 267 + frequencies: Real 268 + response: Complex 269 + magnitude: Real 270 + magnitude_db: Real 271 + 272 + def plot(self, ax: Axes | None = None, **kwargs: Any) -> Axes: 273 + """Plot the weighting factor (dB) versus frequency. 274 + 275 + Requires matplotlib (``pip install phonometry[plot]``); returns the 276 + :class:`~matplotlib.axes.Axes` and never calls ``plt.show``. 277 + """ 278 + from ._plotting import plot_vibration_weighting 279 + 280 + return plot_vibration_weighting(self, ax=ax, **kwargs) 281 + 282 + 283 + def frequency_weighting(name: str, frequencies: ArrayLike) -> WeightingResponse: 284 + """Frequency-weighting response ``H(f)`` (ISO 8041-1:2017, Formula (5)). 285 + 286 + Evaluates the overall weighting - band limiting, acceleration-velocity 287 + transition and upward step - of weighting ``name`` at ``frequencies``. 288 + 289 + :param name: Weighting name (one of :data:`WEIGHTING_NAMES`). 290 + :param frequencies: Frequencies at which to evaluate, in hertz (> 0). 291 + :return: A :class:`WeightingResponse` with ``.plot()``. 292 + :raises ValueError: if ``name`` is unknown or ``frequencies`` is empty. 293 + """ 294 + freq = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) 295 + if freq.ndim != 1 or freq.size == 0: 296 + raise ValueError("'frequencies' must be a non-empty 1-D array.") 297 + resp = _weighting_response(name, freq) 298 + mag = np.abs(resp) 299 + with np.errstate(divide="ignore"): 300 + mag_db = 20.0 * np.log10(mag) 301 + return WeightingResponse( 302 + name=name, 303 + frequencies=freq, 304 + response=resp, 305 + magnitude=mag, 306 + magnitude_db=mag_db, 307 + ) 308 + 309 + 310 + def weighting_factors(name: str, frequencies: ArrayLike) -> Real: 311 + """Weighting factors ``|H(f)|`` of weighting ``name`` (ISO 8041-1). 312 + 313 + Convenience wrapper over :func:`frequency_weighting` returning only the 314 + magnitude array (the ``W_i`` of ISO 2631-1 Eq. (9) / ISO 5349-1 Eq. (A.1)). 315 + 316 + :param name: Weighting name (one of :data:`WEIGHTING_NAMES`). 317 + :param frequencies: Band centre frequencies, in hertz. 318 + :return: Weighting factor per frequency. 319 + """ 320 + return frequency_weighting(name, frequencies).magnitude 321 + 322 + 323 + def apply_weighting(signal: ArrayLike, fs: float, name: str) -> Real: 324 + """Apply frequency weighting ``name`` to a time signal (ISO 8041-1). 325 + 326 + The exact analog response :func:`frequency_weighting` is applied in the 327 + frequency domain (real FFT), so the weighted signal reproduces both the 328 + magnitude and phase of the standard's cascade without bilinear warping. 329 + The multiplication is circular, so the record wraps at its ends; apply it 330 + to a record long enough (or pre-tapered) that the boundary transient is 331 + negligible, as for any block frequency-domain filtering. 332 + 333 + :param signal: Unweighted acceleration time history (1-D), in m/s2. 334 + :param fs: Sampling frequency, in hertz (> 0). 335 + :param name: Weighting name (one of :data:`WEIGHTING_NAMES`). 336 + :return: The frequency-weighted acceleration signal, same length as input. 337 + :raises ValueError: if ``signal`` is not 1-D, ``fs`` is not positive, or 338 + ``name`` is unknown. 339 + """ 340 + x = np.asarray(signal, dtype=np.float64) 341 + if x.ndim != 1 or x.size == 0: 342 + raise ValueError("'signal' must be a non-empty 1-D array.") 343 + fs = _positive_fs(fs) 344 + n = x.size 345 + freqs = np.fft.rfftfreq(n, d=1.0 / fs).astype(np.float64) 346 + resp = _weighting_response(name, freqs) 347 + weighted = np.fft.irfft(np.fft.rfft(x) * resp, n=n) 348 + return np.asarray(weighted, dtype=np.float64) 349 + 350 + 351 + # --------------------------------------------------------------------------- 352 + # Band (spectrum) method (ISO 2631-1 Eq. (9); ISO 5349-1 Eq. (A.1)). 353 + # --------------------------------------------------------------------------- 354 + @dataclass(frozen=True) 355 + class WeightedSpectrum: 356 + """A weighted one-third-octave acceleration spectrum and its ``a_w``. 357 + 358 + :ivar frequencies: Band centre frequencies, in hertz. 359 + :ivar band_accelerations: Unweighted r.m.s. acceleration per band, in m/s2. 360 + :ivar weighting_name: Weighting applied (one of :data:`WEIGHTING_NAMES`). 361 + :ivar weighting_factors: Weighting factor ``W_i`` per band. 362 + :ivar weighted: Weighted band contribution ``W_i*a_i``, in m/s2. 363 + :ivar overall: Overall weighted r.m.s. acceleration ``a_w``, in m/s2. 364 + """ 365 + 366 + frequencies: Real 367 + band_accelerations: Real 368 + weighting_name: str 369 + weighting_factors: Real 370 + weighted: Real 371 + overall: float 372 + 373 + def plot(self, ax: Axes | None = None, **kwargs: Any) -> Axes: 374 + """Plot the unweighted and weighted band spectra with ``a_w``. 375 + 376 + Requires matplotlib (``pip install phonometry[plot]``); returns the 377 + :class:`~matplotlib.axes.Axes` and never calls ``plt.show``. 378 + """ 379 + from ._plotting import plot_weighted_spectrum 380 + 381 + return plot_weighted_spectrum(self, ax=ax, **kwargs) 382 + 383 + 384 + def weighted_acceleration( 385 + band_accelerations: ArrayLike, 386 + frequencies: ArrayLike, 387 + weighting: str, 388 + ) -> WeightedSpectrum: 389 + """Weighted r.m.s. acceleration from a band spectrum (ISO 2631-1 Eq. (9)). 390 + 391 + ``a_w = sqrt( sum_i (W_i * a_i)^2 )`` with the per-band weighting factors 392 + ``W_i`` of ISO 8041-1 evaluated at the band centres (ISO 5349-1 Eq. (A.1) 393 + is the identical construction for the hand-arm weighting ``Wh``). 394 + 395 + :param band_accelerations: r.m.s. acceleration ``a_i`` per band, in m/s2. 396 + :param frequencies: Band centre frequencies, in hertz. 397 + :param weighting: Weighting name (one of :data:`WEIGHTING_NAMES`). 398 + :return: A :class:`WeightedSpectrum` with ``.plot()``. 399 + :raises ValueError: if the inputs differ in length or are empty. 400 + """ 401 + accel = np.atleast_1d(np.asarray(band_accelerations, dtype=np.float64)) 402 + freq = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) 403 + if accel.ndim != 1 or accel.size == 0 or accel.shape != freq.shape: 404 + raise ValueError( 405 + "'band_accelerations' and 'frequencies' must be non-empty, 1-D and " 406 + "equal-length." 407 + ) 408 + factors = weighting_factors(weighting, freq) 409 + weighted = factors * accel 410 + overall = float(np.sqrt(np.sum(weighted**2))) 411 + return WeightedSpectrum( 412 + frequencies=freq, 413 + band_accelerations=accel, 414 + weighting_name=weighting, 415 + weighting_factors=factors, 416 + weighted=weighted, 417 + overall=overall, 418 + ) 419 + 420 + 421 + # --------------------------------------------------------------------------- 422 + # Time-domain metrics (ISO 2631-1 clause 6; ISO 8041-1 clause 3.1.2). 423 + # --------------------------------------------------------------------------- 424 + def _weighted_signal(signal: ArrayLike) -> Real: 425 + x = np.asarray(signal, dtype=np.float64) 426 + if x.ndim != 1 or x.size == 0: 427 + raise ValueError("'signal' must be a non-empty 1-D array.") 428 + return x 429 + 430 + 431 + def _positive_fs(fs: float) -> float: 432 + """Return ``fs`` as a positive, finite float or raise ``ValueError``.""" 433 + fs = float(fs) 434 + if not math.isfinite(fs) or fs <= 0.0: 435 + raise ValueError("'fs' must be a positive, finite sampling frequency.") 436 + return fs 437 + 438 + 439 + def running_rms( 440 + signal: ArrayLike, 441 + fs: float, 442 + *, 443 + integration_time: float = 1.0, 444 + method: str = "linear", 445 + ) -> Real: 446 + """Running r.m.s. of a weighted signal (ISO 2631-1 Eqs. (2)/(3)). 447 + 448 + :param signal: Frequency-weighted acceleration signal (1-D), in m/s2. 449 + :param fs: Sampling frequency, in hertz. 450 + :param integration_time: Averaging time ``tau``, in seconds (default 1 s, 451 + the "slow" constant of ISO 2631-1 6.3.1). 452 + :param method: ``"linear"`` (Eq. (2), a sliding rectangular average) or 453 + ``"exponential"`` (Eq. (3), a single-pole average). 454 + :return: The running r.m.s. value ``a_w(t0)`` per sample, in m/s2. 455 + :raises ValueError: for a bad signal, non-positive ``fs``/``tau`` or an 456 + unknown ``method``. 457 + """ 458 + x = _weighted_signal(signal) 459 + fs = _positive_fs(fs) 460 + tau = float(integration_time) 461 + if not math.isfinite(tau) or tau <= 0.0: 462 + raise ValueError("'integration_time' must be positive and finite.") 463 + power = x**2 464 + if method == "linear": 465 + window = max(1, int(round(tau * fs))) 466 + # Causal sliding mean over the trailing ``tau`` seconds, O(n) via a 467 + # prefix sum (early samples average over the samples available, i.e. 468 + # zero-padded at the front, matching the "slow" running r.m.s.). 469 + csum = np.concatenate(([0.0], np.cumsum(power))) 470 + lo = np.maximum(0, np.arange(power.size) + 1 - window) 471 + mean_power = (csum[1:] - csum[lo]) / window 472 + elif method == "exponential": 473 + # Single-pole IIR average: y[i] = (1-alpha) y[i-1] + alpha x[i]. 474 + alpha = 1.0 - math.exp(-1.0 / (tau * fs)) 475 + mean_power = sig.lfilter([alpha], [1.0, -(1.0 - alpha)], power) 476 + else: 477 + raise ValueError("'method' must be 'linear' or 'exponential'.") 478 + return np.sqrt(mean_power) 479 + 480 + 481 + def mtvv(signal: ArrayLike, fs: float, *, integration_time: float = 1.0) -> float: 482 + """Maximum transient vibration value (ISO 2631-1 Eq. (4)). 483 + 484 + ``MTVV = max a_w(t0)``, the peak of the 1 s running r.m.s. value. 485 + 486 + :param signal: Frequency-weighted acceleration signal (1-D), in m/s2. 487 + :param fs: Sampling frequency, in hertz. 488 + :param integration_time: Running-r.m.s. averaging time, in seconds (1 s). 489 + :return: The MTVV, in m/s2. 490 + """ 491 + return float( 492 + np.max(running_rms(signal, fs, integration_time=integration_time)) 493 + ) 494 + 495 + 496 + def vibration_dose_value(signal: ArrayLike, fs: float) -> float: 497 + """Vibration dose value ``VDV`` (ISO 2631-1 Eq. (5)). 498 + 499 + ``VDV = ( integral a_w(t)^4 dt )^(1/4)``, in m/s^1,75; more sensitive to 500 + peaks than the r.m.s. value. 501 + 502 + :param signal: Frequency-weighted acceleration signal (1-D), in m/s2. 503 + :param fs: Sampling frequency, in hertz. 504 + :return: The VDV, in m/s^1,75. 505 + :raises ValueError: for a bad signal or non-positive ``fs``. 506 + """ 507 + x = _weighted_signal(signal) 508 + fs = _positive_fs(fs) 509 + return float(np.sum(x**4 / fs) ** 0.25) 510 + 511 + 512 + def motion_sickness_dose_value(signal: ArrayLike, fs: float) -> float: 513 + """Motion sickness dose value ``MSDV`` (ISO 2631-1 clause 9; 8041-1 3.1.2.5). 514 + 515 + ``MSDV = ( integral a_w(t)^2 dt )^(1/2)``, in m/s^1,5; the ``Wf``-weighted 516 + signal is the intended input. 517 + 518 + :param signal: Frequency-weighted acceleration signal (1-D), in m/s2. 519 + :param fs: Sampling frequency, in hertz. 520 + :return: The MSDV, in m/s^1,5. 521 + :raises ValueError: for a bad signal or non-positive ``fs``. 522 + """ 523 + x = _weighted_signal(signal) 524 + fs = _positive_fs(fs) 525 + return float(np.sum(x**2 / fs) ** 0.5) 526 + 527 + 528 + def crest_factor(signal: ArrayLike) -> float: 529 + """Crest factor of a weighted signal (ISO 2631-1 clause 6.2.1). 530 + 531 + The modulus of the ratio of the peak weighted acceleration to its r.m.s. 532 + value. ISO 2631-1 6.2.2 deems the basic (r.m.s.) method adequate for a 533 + crest factor up to 9. 534 + 535 + Emits a :class:`HumanVibrationWarning` when the crest factor exceeds 9, 536 + the threshold above which ISO 2631-1 6.2.2 deems the basic method 537 + inadequate. 538 + 539 + :param signal: Frequency-weighted acceleration signal (1-D), in m/s2. 540 + :return: The crest factor (dimensionless); ``0`` for an all-zero signal. 541 + """ 542 + x = _weighted_signal(signal) 543 + rms = float(np.sqrt(np.mean(x**2))) 544 + # rms is a root-mean-square, so it is >= 0; <= 0 means an all-zero signal. 545 + if rms <= 0.0: 546 + return 0.0 547 + cf = float(np.max(np.abs(x)) / rms) 548 + if cf > 9.0: 549 + warnings.warn( 550 + f"Crest factor {cf:.1f} exceeds 9; the basic r.m.s. method may be " 551 + "inadequate (ISO 2631-1, 6.2.2) - consider the VDV or running " 552 + "r.m.s. dose measures.", 553 + HumanVibrationWarning, 554 + stacklevel=2, 555 + ) 556 + return cf 557 + 558 + 559 + # --------------------------------------------------------------------------- 560 + # Vector sum and daily exposure (ISO 2631-1 Eq. (10); ISO 5349-1 Eqs. (1)-(3)). 561 + # --------------------------------------------------------------------------- 562 + def vibration_total_value( 563 + components: ArrayLike, 564 + *, 565 + k: ArrayLike | None = None, 566 + ) -> float: 567 + """Vibration total value ``a_v`` / ``a_hv`` (ISO 2631-1 Eq. (10)). 568 + 569 + ``a_v = sqrt( sum_j k_j^2 * a_wj^2 )`` over the (up to three) axis-weighted 570 + r.m.s. accelerations. With ``k = None`` the unweighted vector sum of 571 + ISO 5349-1 Eq. (1) (``a_hv``, all ``k = 1``) is returned. 572 + 573 + :param components: Axis-weighted r.m.s. accelerations ``a_wj``, in m/s2. 574 + :param k: Optional per-axis multiplying factors ``k_j`` (ISO 2631-1 7.2.3); 575 + ``None`` uses unity for every axis. 576 + :return: The vibration total value, in m/s2. 577 + :raises ValueError: if ``components`` is empty or ``k`` differs in length. 578 + """ 579 + comp = np.atleast_1d(np.asarray(components, dtype=np.float64)) 580 + if comp.ndim != 1 or comp.size == 0: 581 + raise ValueError("'components' must be a non-empty 1-D array.") 582 + if k is None: 583 + factors = np.ones_like(comp) 584 + else: 585 + factors = np.atleast_1d(np.asarray(k, dtype=np.float64)) 586 + if factors.shape != comp.shape: 587 + raise ValueError("'k' must have the same length as 'components'.") 588 + return float(np.sqrt(np.sum((factors * comp) ** 2))) 589 + 590 + 591 + def daily_exposure(total_value: float, duration_s: float) -> float: 592 + """Daily exposure ``A(8)`` for one operation (ISO 5349-1 Eq. (2)). 593 + 594 + ``A(8) = a_hv * sqrt(T / T0)`` with ``T0 = 8 h``. The identical form gives 595 + the whole-body ``A(8)`` used by Directive 2002/44/EC. 596 + 597 + :param total_value: Vibration total value ``a_hv`` (or ``a_v``), in m/s2. 598 + :param duration_s: Daily exposure duration ``T``, in seconds. 599 + :return: The daily exposure ``A(8)``, in m/s2. 600 + :raises ValueError: for a negative magnitude or duration. 601 + """ 602 + a = float(total_value) 603 + t = float(duration_s) 604 + if a < 0.0 or t < 0.0: 605 + raise ValueError("'total_value' and 'duration_s' must be non-negative.") 606 + return a * math.sqrt(t / REFERENCE_DURATION_S) 607 + 608 + 609 + #: Alias for :func:`daily_exposure` reading naturally in a hand-arm context. 610 + partial_exposure = daily_exposure 611 + 612 + 613 + def combine_partial_exposures(partials: ArrayLike) -> float: 614 + """Combine partial exposures into ``A(8)`` (ISO 5349-1/-2 Eq. (3)). 615 + 616 + ``A(8) = sqrt( sum_i A_i(8)^2 )``. 617 + 618 + :param partials: Partial exposures ``A_i(8)``, in m/s2. 619 + :return: The combined daily exposure ``A(8)``, in m/s2. 620 + :raises ValueError: if ``partials`` is empty. 621 + """ 622 + parts = np.atleast_1d(np.asarray(partials, dtype=np.float64)) 623 + if parts.ndim != 1 or parts.size == 0: 624 + raise ValueError("'partials' must be a non-empty 1-D array.") 625 + return float(np.sqrt(np.sum(parts**2))) 626 + 627 + 628 + def hav_daily_exposure( 629 + total_values: ArrayLike, 630 + durations_s: ArrayLike, 631 + ) -> float: 632 + """Daily exposure ``A(8)`` for several operations (ISO 5349-1 Eq. (3)). 633 + 634 + ``A(8) = sqrt( (1/T0) * sum_i a_hvi^2 * T_i )``. 635 + 636 + :param total_values: Vibration total value ``a_hvi`` per operation, in m/s2. 637 + :param durations_s: Duration ``T_i`` per operation, in seconds. 638 + :return: The daily exposure ``A(8)``, in m/s2. 639 + :raises ValueError: if the inputs differ in length, are empty or negative. 640 + """ 641 + ahv = np.atleast_1d(np.asarray(total_values, dtype=np.float64)) 642 + t = np.atleast_1d(np.asarray(durations_s, dtype=np.float64)) 643 + if ahv.ndim != 1 or ahv.size == 0 or ahv.shape != t.shape: 644 + raise ValueError( 645 + "'total_values' and 'durations_s' must be non-empty, 1-D and " 646 + "equal-length." 647 + ) 648 + if np.any(ahv < 0.0) or np.any(t < 0.0): 649 + raise ValueError("'total_values' and 'durations_s' must be non-negative.") 650 + return float(np.sqrt(np.sum(ahv**2 * t) / REFERENCE_DURATION_S)) 651 + 652 + 653 + def energy_equivalent_acceleration( 654 + magnitudes: ArrayLike, 655 + durations_s: ArrayLike, 656 + ) -> float: 657 + """Energy-equivalent weighted acceleration (ISO 2631-1 Eq. (B.3)). 658 + 659 + ``a_w,e = sqrt( sum a_wi^2 * T_i / sum T_i )``. 660 + 661 + :param magnitudes: Weighted r.m.s. magnitudes ``a_wi``, in m/s2. 662 + :param durations_s: Duration ``T_i`` per period, in seconds. 663 + :return: The energy-equivalent magnitude ``a_w,e``, in m/s2. 664 + :raises ValueError: if the inputs differ in length, are empty, or the total 665 + duration is zero. 666 + """ 667 + a = np.atleast_1d(np.asarray(magnitudes, dtype=np.float64)) 668 + t = np.atleast_1d(np.asarray(durations_s, dtype=np.float64)) 669 + if a.ndim != 1 or a.size == 0 or a.shape != t.shape: 670 + raise ValueError( 671 + "'magnitudes' and 'durations_s' must be non-empty, 1-D and " 672 + "equal-length." 673 + ) 674 + if np.any(t < 0.0): 675 + raise ValueError("'durations_s' must be non-negative.") 676 + total = float(np.sum(t)) 677 + if total <= 0.0: 678 + raise ValueError("The total duration must be positive.") 679 + return float(np.sqrt(np.sum(a**2 * t) / total)) 680 + 681 + 682 + def hav_vwf_lifetime_years(a8: float) -> float: 683 + """Years to 10 % vibration-white-finger prevalence (ISO 5349-1 Eq. (C.1)). 684 + 685 + ``Dy = 31,8 * A(8)^(-1,06)`` - the group-mean lifetime exposure that 686 + produces finger blanching in 10 % of an exposed group (informative 687 + Annex C). 688 + 689 + :param a8: Daily vibration exposure ``A(8)``, in m/s2 (> 0). 690 + :return: The lifetime exposure duration ``Dy``, in years. 691 + :raises ValueError: if ``a8`` is not positive. 692 + """ 693 + a = float(a8) 694 + if not math.isfinite(a) or a <= 0.0: 695 + raise ValueError("'a8' must be a positive daily exposure.") 696 + return float(31.8 * a**-1.06) 697 + 698 + 699 + # --------------------------------------------------------------------------- 700 + # Exposure assessment against Directive 2002/44/EC action / limit values. 701 + # --------------------------------------------------------------------------- 702 + @dataclass(frozen=True) 703 + class ExposureAssessment: 704 + """A daily exposure assessed against the Directive 2002/44/EC values. 705 + 706 + :ivar value: The assessed daily exposure ``A(8)`` (or VDV), in its unit. 707 + :ivar kind: ``"hav"`` or ``"wbv"``. 708 + :ivar metric: ``"a8"`` (m/s2) or ``"vdv"`` (m/s^1,75). 709 + :ivar action_value: The exposure action value (EAV). 710 + :ivar limit_value: The exposure limit value (ELV). 711 + :ivar exceeds_action: Whether ``value`` reaches or exceeds the EAV. 712 + :ivar exceeds_limit: Whether ``value`` reaches or exceeds the ELV. 713 + :ivar zone: ``"below action"``, ``"action"`` (EAV<=value<ELV) or 714 + ``"limit"`` (value>=ELV). 715 + """ 716 + 717 + value: float 718 + kind: str 719 + metric: str 720 + action_value: float 721 + limit_value: float 722 + exceeds_action: bool 723 + exceeds_limit: bool 724 + zone: str 725 + 726 + 727 + def exposure_assessment( 728 + value: float, 729 + *, 730 + kind: str, 731 + metric: str = "a8", 732 + ) -> ExposureAssessment: 733 + """Assess a daily exposure against Directive 2002/44/EC (Article 3). 734 + 735 + :param value: The daily exposure ``A(8)`` in m/s2 (``metric="a8"``) or the 736 + vibration dose value in m/s^1,75 (``metric="vdv"``, whole-body only). 737 + :param kind: ``"hav"`` (hand-arm) or ``"wbv"`` (whole-body). 738 + :param metric: ``"a8"`` (default) or ``"vdv"``. 739 + :return: An :class:`ExposureAssessment`. 740 + :raises ValueError: for an unknown ``kind``/``metric`` combination or a 741 + negative ``value``. 742 + """ 743 + v = float(value) 744 + if v < 0.0: 745 + raise ValueError("'value' must be non-negative.") 746 + if kind == "hav" and metric == "a8": 747 + eav, elv = HAV_EAV_A8, HAV_ELV_A8 748 + elif kind == "wbv" and metric == "a8": 749 + eav, elv = WBV_EAV_A8, WBV_ELV_A8 750 + elif kind == "wbv" and metric == "vdv": 751 + eav, elv = WBV_EAV_VDV, WBV_ELV_VDV 752 + else: 753 + raise ValueError( 754 + "kind must be 'hav' or 'wbv'; metric 'a8' (both) or 'vdv' (wbv only)." 755 + ) 756 + exceeds_action = v >= eav 757 + exceeds_limit = v >= elv 758 + if exceeds_limit: 759 + zone = "limit" 760 + elif exceeds_action: 761 + zone = "action" 762 + else: 763 + zone = "below action" 764 + return ExposureAssessment( 765 + value=v, 766 + kind=kind, 767 + metric=metric, 768 + action_value=eav, 769 + limit_value=elv, 770 + exceeds_action=exceeds_action, 771 + exceeds_limit=exceeds_limit, 772 + zone=zone, 773 + ) 774 + 775 + 776 + @dataclass(frozen=True) 777 + class DailyVibrationExposure: 778 + """A daily exposure built from several operations, with its assessment. 779 + 780 + :ivar a8: The daily exposure ``A(8)``, in m/s2. 781 + :ivar labels: A label per operation. 782 + :ivar total_values: Vibration total value ``a_hvi`` per operation, in m/s2. 783 + :ivar durations_s: Duration ``T_i`` per operation, in seconds. 784 + :ivar partials: Partial exposure ``A_i(8)`` per operation, in m/s2. 785 + :ivar assessment: The :class:`ExposureAssessment` of ``a8``. 786 + """ 787 + 788 + a8: float 789 + labels: tuple[str, ...] 790 + total_values: Real 791 + durations_s: Real 792 + partials: Real 793 + assessment: ExposureAssessment 794 + 795 + def plot(self, ax: Axes | None = None, **kwargs: Any) -> Axes: 796 + """Plot the partial exposures against the EAV / ELV thresholds. 797 + 798 + Requires matplotlib (``pip install phonometry[plot]``); returns the 799 + :class:`~matplotlib.axes.Axes` and never calls ``plt.show``. 800 + """ 801 + from ._plotting import plot_daily_exposure 802 + 803 + return plot_daily_exposure(self, ax=ax, **kwargs) 804 + 805 + 806 + def daily_vibration_exposure( 807 + total_values: ArrayLike, 808 + durations_s: ArrayLike, 809 + *, 810 + kind: str, 811 + labels: "list[str] | tuple[str, ...] | None" = None, 812 + ) -> DailyVibrationExposure: 813 + """Daily exposure from several operations, assessed (ISO 5349 + Directive). 814 + 815 + Combines the operations via the partial exposures of ISO 5349-1/-2 816 + Eqs. (2)/(3) and assesses the resulting ``A(8)`` against the Directive 817 + 2002/44/EC action and limit values for ``kind``. 818 + 819 + :param total_values: Vibration total value ``a_hvi`` per operation, in m/s2. 820 + :param durations_s: Duration ``T_i`` per operation, in seconds. 821 + :param kind: ``"hav"`` or ``"wbv"`` (selects the EAV/ELV). 822 + :param labels: Optional operation labels; defaults to ``op 1``, ``op 2``, ... 823 + :return: A :class:`DailyVibrationExposure` with ``.plot()``. 824 + :raises ValueError: if the inputs differ in length or ``labels`` mismatches. 825 + """ 826 + ahv = np.atleast_1d(np.asarray(total_values, dtype=np.float64)) 827 + t = np.atleast_1d(np.asarray(durations_s, dtype=np.float64)) 828 + if ahv.ndim != 1 or ahv.size == 0 or ahv.shape != t.shape: 829 + raise ValueError( 830 + "'total_values' and 'durations_s' must be non-empty, 1-D and " 831 + "equal-length." 832 + ) 833 + if labels is None: 834 + labels = tuple(f"op {i + 1}" for i in range(ahv.size)) 835 + else: 836 + labels = tuple(labels) 837 + if len(labels) != ahv.size: 838 + raise ValueError("'labels' must match the number of operations.") 839 + partials = np.array( 840 + [daily_exposure(float(a), float(dt)) for a, dt in zip(ahv, t)], 841 + dtype=np.float64, 842 + ) 843 + a8 = combine_partial_exposures(partials) 844 + return DailyVibrationExposure( 845 + a8=a8, 846 + labels=labels, 847 + total_values=ahv, 848 + durations_s=t, 849 + partials=partials, 850 + assessment=exposure_assessment(a8, kind=kind, metric="a8"), 851 + )
+31 -6
tests/reference_data.py
··· 5 5 suite (``tests/test_*.py``) and the CI conformance report 6 6 (``scripts/conformance_report.py``) import these constants, so the report's 7 7 expected values can never drift from what the tests assert. The PR-B 8 - building-acoustics and PR-E scattering/in-situ/precision-power oracles are the 9 - exception: their test modules re-hardcode the values inline rather than import 10 - them, and dedicated consistency tests 11 - (``test_building_reference_data_matches_published_oracles`` and 12 - ``test_scattering_insitu_precision_reference_data_matches_oracles``) pin this 13 - shared table to those same published results so neither copy can drift. 8 + building-acoustics, PR-E scattering/in-situ/precision-power and PR-F 9 + human-vibration oracles are the exception: their test modules re-hardcode the 10 + values inline rather than import them, and dedicated consistency tests 11 + (``test_building_reference_data_matches_published_oracles``, 12 + ``test_scattering_insitu_precision_reference_data_matches_oracles`` and 13 + ``test_human_vibration_reference_data_matches_oracles``) pin this shared table 14 + to those same published results so neither copy can drift. 14 15 15 16 This module is deliberately dependency-free (stdlib only) so it can be 16 17 imported in the ``pr-comment`` CI job, which installs the runtime ··· 423 424 ISO9614_3_UNIFORM_POWER = 1.0e-4 # radiated power W (W) 424 425 ISO9614_3_UNIFORM_AREAS: tuple[float, ...] = (0.5, 1.0, 0.25, 2.0) 425 426 ISO9614_3_UNIFORM_LW = 80.0 # 10*lg(W/1e-12) (dB) 427 + 428 + # --------------------------------------------------------------------------- 429 + # PR-F human vibration (ISO 8041-1 / ISO 2631 / ISO 5349 / Directive 2002/44/EC). 430 + # The true IEC 61260 one-third-octave centre is 10^(n/10) Hz; the reference 431 + # frequencies of ISO 8041-1 Table 1 are exact (rad/s -> Hz). Design-goal 432 + # factors are from ISO 8041-1:2017 Annex B (Tables B.1-B.9, 4 sig. figs). 433 + # --------------------------------------------------------------------------- 434 + # ISO 8041-1 Annex B design-goal weighting factors at the true band centre. 435 + ISO8041_1_WK_FACTOR_6P31HZ = 1.054 # Table B.8, n = 8 (6,31 Hz) - Wk peak 436 + ISO8041_1_WM_FACTOR_1P585HZ = 0.9342 # Table B.9, n = 2 (1,585 Hz) - Wm 437 + # ISO 8041-1 Table 1 weighting factor at the reference frequency. 438 + ISO8041_1_WH_REF_FREQ_HZ = 500.0 / (2.0 * math.pi) # 79,577 Hz (500 rad/s) 439 + ISO8041_1_WH_REF_FACTOR = 0.2020 # Table 1, Wh @ 500 rad/s 440 + # ISO 5349-2:2001 Annex E worked-example daily exposures A(8), m/s^2. 441 + ISO5349_2_E21_A8 = 4.1 # E.2.1 single tool: 7,4*sqrt(2,5/8) 442 + ISO5349_2_E3_A8 = 3.6 # E.3 forestry three-task combination 443 + # ISO 5349-1:2001 Annex C: Dy = 31,8*A(8)^-1,06; Table C.1 A(8)=7 -> Dy=4 yr. 444 + ISO5349_1_VWF_A8 = 7.0 445 + ISO5349_1_VWF_DY_YEARS = 4.0 446 + # Directive 2002/44/EC Article 3 daily exposure action/limit values. 447 + DIRECTIVE_2002_44_HAV_EAV = 2.5 # A(8) m/s^2, Art. 3(1)(a) 448 + DIRECTIVE_2002_44_HAV_ELV = 5.0 # A(8) m/s^2, Art. 3(1)(b) 449 + DIRECTIVE_2002_44_WBV_EAV = 0.5 # A(8) m/s^2, Art. 3(2)(a) 450 + DIRECTIVE_2002_44_WBV_ELV = 1.15 # A(8) m/s^2, Art. 3(2)(b)
+29
tests/test_conformance_report.py
··· 155 155 ) 156 156 157 157 158 + def test_human_vibration_checks_registered() -> None: 159 + """The PR-F human-vibration checks are wired into their own domain.""" 160 + standards = {c.standard for c in cr.CHECKS} 161 + assert "ISO 8041-1:2017 Table B.8" in standards # Wk design-goal factor 162 + assert "ISO 8041-1:2017 Table 1" in standards # Wh reference-freq factor 163 + assert "ISO 5349-2:2001 Example E.2.1" in standards # A(8) single tool 164 + assert "ISO 5349-2:2001 Example E.3" in standards # A(8) forestry 165 + assert "ISO 5349-1:2001 Eq. (C.1)" in standards # VWF lifetime 166 + assert "Directive 2002/44/EC Art. 3" in standards # EAV/ELV 167 + assert "Human vibration (ISO 8041 / 2631 / 5349)" in cr._domains() 168 + 169 + 170 + def test_human_vibration_reference_data_matches_oracles() -> None: 171 + """Pin the shared PR-F constants to their published values.""" 172 + import reference_data as ref 173 + 174 + assert ref.ISO8041_1_WK_FACTOR_6P31HZ == 1.054 # Table B.8, n=8 175 + assert ref.ISO8041_1_WM_FACTOR_1P585HZ == 0.9342 # Table B.9, n=2 176 + assert ref.ISO8041_1_WH_REF_FACTOR == 0.2020 # Table 1, Wh @ 500 rad/s 177 + assert ref.ISO5349_2_E21_A8 == 4.1 # 7,4*sqrt(2,5/8) 178 + assert ref.ISO5349_2_E3_A8 == 3.6 # forestry three-task 179 + assert ref.ISO5349_1_VWF_A8 == 7.0 180 + assert ref.ISO5349_1_VWF_DY_YEARS == 4.0 # Table C.1 181 + assert ref.DIRECTIVE_2002_44_HAV_EAV == 2.5 182 + assert ref.DIRECTIVE_2002_44_HAV_ELV == 5.0 183 + assert ref.DIRECTIVE_2002_44_WBV_EAV == 0.5 184 + assert ref.DIRECTIVE_2002_44_WBV_ELV == 1.15 185 + 186 + 158 187 @pytest.mark.parametrize("check", cr.CHECKS, ids=lambda c: f"{c.standard} :: {c.quantity}") 159 188 def test_every_check_passes(check: "cr.Check") -> None: 160 189 outcome = check.run()
+458
tests/test_human_vibration.py
··· 1 + # Copyright (c) 2026. Jose M. Requena-Plens 2 + """Tests for :mod:`phonometry.human_vibration`. 3 + 4 + The frequency weightings are validated against the ISO 8041-1:2017 Annex B 5 + design-goal factors (Tables B.1-B.9) and the Table 1 reference-frequency 6 + factors, cross-checked against the ISO 2631-1/-2 and ISO 5349-1 tabulated 7 + weightings. The exposure arithmetic reproduces the ISO 5349-2 Annex E worked 8 + examples, and the assessment follows the Directive 2002/44/EC action/limit 9 + values. 10 + """ 11 + 12 + from __future__ import annotations 13 + 14 + import math 15 + 16 + import numpy as np 17 + import pytest 18 + 19 + from phonometry import human_vibration as hv 20 + 21 + 22 + def _fc(n: int) -> float: 23 + """True IEC 61260 one-third-octave centre ``10^(n/10)`` Hz.""" 24 + return 10.0 ** (n / 10.0) 25 + 26 + 27 + # --------------------------------------------------------------------------- 28 + # Frequency weightings vs ISO 8041-1:2017 Annex B (Tables B.1-B.9). 29 + # --------------------------------------------------------------------------- 30 + # (name, band index n, design-goal factor) sampled across every weighting. 31 + _ANNEX_B = [ 32 + ("Wk", 0, 0.4825), ("Wk", 8, 1.054), ("Wk", 20, 0.08873), ("Wk", -10, 0.03121), 33 + ("Wd", 0, 1.011), ("Wd", 13, 0.1004), ("Wd", 26, 0.0003164), 34 + ("Wc", 1, 1.000), ("Wc", 8, 0.9739), ("Wc", 20, 0.05665), 35 + ("We", 0, 0.8798), ("We", 7, 0.2012), ("We", 20, 0.007071), 36 + ("Wj", 15, 1.000), ("Wj", 0, 0.4844), ("Wj", 20, 0.7075), 37 + ("Wb", 8, 1.054), ("Wb", 0, 0.3853), ("Wb", 20, 0.1154), 38 + ("Wf", -8, 1.004), ("Wf", -10, 0.6951), ("Wf", 0, 0.02352), 39 + ("Wh", 8, 0.7272), ("Wh", 10, 0.9514), ("Wh", 20, 0.1602), ("Wh", 30, 0.01346), 40 + ("Wm", 2, 0.9342), ("Wm", 10, 0.4941), ("Wm", -1, 0.7003), ("Wm", 20, 0.04013), 41 + ] 42 + 43 + 44 + @pytest.mark.parametrize(("name", "n", "expected"), _ANNEX_B) 45 + def test_annex_b_design_goal_factors(name: str, n: int, expected: float) -> None: 46 + """|H| at the true band centre matches the Annex B four-figure factor.""" 47 + got = hv.weighting_factors(name, _fc(n))[0] 48 + # 0,1 % relative + a floor for the deep-attenuation four-figure entries. 49 + assert got == pytest.approx(expected, rel=1e-3, abs=5e-6) 50 + 51 + 52 + # (name, exact reference frequency Hz, Table 1 weighting factor at ref). 53 + _TABLE_1 = [ 54 + ("Wh", 500.0 / (2.0 * math.pi), 0.2020), 55 + ("Wk", 100.0 / (2.0 * math.pi), 0.7718), 56 + ("Wb", 100.0 / (2.0 * math.pi), 0.8126), 57 + ("Wc", 100.0 / (2.0 * math.pi), 0.5145), 58 + ("Wd", 100.0 / (2.0 * math.pi), 0.1261), 59 + ("We", 100.0 / (2.0 * math.pi), 0.06287), 60 + ("Wj", 100.0 / (2.0 * math.pi), 1.019), 61 + ("Wm", 100.0 / (2.0 * math.pi), 0.3362), 62 + ("Wf", 2.5 / (2.0 * math.pi), 0.3888), 63 + ] 64 + 65 + 66 + @pytest.mark.parametrize(("name", "freq", "expected"), _TABLE_1) 67 + def test_table_1_reference_frequency_factors( 68 + name: str, freq: float, expected: float 69 + ) -> None: 70 + """|H| at the Table 1 reference frequency matches the tabulated factor.""" 71 + got = hv.weighting_factors(name, freq)[0] 72 + assert got == pytest.approx(expected, rel=1.5e-3) 73 + 74 + 75 + def test_iso5349_1_wh_third_octave_table_a2() -> None: 76 + """The Wh factors match ISO 5349-1 Table A.2 to its three figures.""" 77 + # (nominal band centre, Whi) from ISO 5349-1:2001 Table A.2. 78 + a2 = { 79 + 6.3: 0.727, 8.0: 0.873, 10.0: 0.951, 12.5: 0.958, 16.0: 0.896, 80 + 20.0: 0.782, 25.0: 0.647, 31.5: 0.519, 40.0: 0.411, 100.0: 0.160, 81 + 250.0: 0.0634, 1000.0: 0.0135, 82 + } 83 + # Evaluate at the *true* centres the table is computed at. 84 + n_of = {6.3: 8, 8.0: 9, 10.0: 10, 12.5: 11, 16.0: 12, 20.0: 13, 25.0: 14, 85 + 31.5: 15, 40.0: 16, 100.0: 20, 250.0: 24, 1000.0: 30} 86 + for nominal, expected in a2.items(): 87 + got = hv.weighting_factors("Wh", _fc(n_of[nominal]))[0] 88 + assert got == pytest.approx(expected, rel=2e-3, abs=5e-4) 89 + 90 + 91 + def test_iso2631_2_wm_table_a1() -> None: 92 + """Wm matches ISO 2631-2 Table A.1 at representative true centres.""" 93 + a1 = {2: 0.934, 3: 0.932, 10: 0.494, 18: 0.0834, -3: 0.368} 94 + for n, expected in a1.items(): 95 + got = hv.weighting_factors("Wm", _fc(n))[0] 96 + assert got == pytest.approx(expected, rel=2e-3, abs=5e-4) 97 + 98 + 99 + def test_weighting_response_fields_and_db() -> None: 100 + resp = hv.frequency_weighting("Wk", [1.0, 6.3096, 100.0]) 101 + assert resp.name == "Wk" 102 + assert resp.frequencies.shape == (3,) 103 + assert np.allclose(resp.magnitude, np.abs(resp.response)) 104 + assert np.allclose(resp.magnitude_db, 20.0 * np.log10(resp.magnitude)) 105 + assert resp.plot # method exists 106 + 107 + 108 + def test_unknown_weighting_raises() -> None: 109 + with pytest.raises(ValueError, match="Unknown weighting"): 110 + hv.frequency_weighting("Wz", [10.0]) 111 + 112 + 113 + def test_frequency_weighting_rejects_empty() -> None: 114 + with pytest.raises(ValueError, match="non-empty"): 115 + hv.frequency_weighting("Wk", []) 116 + 117 + 118 + def test_dc_and_negative_frequencies_are_blocked() -> None: 119 + resp = hv.frequency_weighting("Wk", [0.0, -1.0, 1.0]) 120 + assert resp.magnitude[0] == 0.0 # DC blocked by the high-pass 121 + assert resp.magnitude[1] == 0.0 122 + assert resp.magnitude[2] > 0.0 123 + 124 + 125 + # --------------------------------------------------------------------------- 126 + # apply_weighting (frequency-domain application of the exact response). 127 + # --------------------------------------------------------------------------- 128 + def test_apply_weighting_scales_sine_by_magnitude() -> None: 129 + fs = 2000.0 130 + f0 = 80.0 131 + t = np.arange(int(4 * fs)) / fs 132 + x = np.sqrt(2.0) * np.sin(2.0 * math.pi * f0 * t) # unit-r.m.s. amplitude 133 + y = hv.apply_weighting(x, fs, "Wk") 134 + factor = hv.weighting_factors("Wk", f0)[0] 135 + # Interior r.m.s. (drop edges) equals |H(f0)| * input r.m.s. (~1). 136 + interior = y[int(0.5 * fs):-int(0.5 * fs)] 137 + assert float(np.sqrt(np.mean(interior**2))) == pytest.approx(factor, rel=2e-2) 138 + 139 + 140 + def test_apply_weighting_validates() -> None: 141 + with pytest.raises(ValueError, match="1-D"): 142 + hv.apply_weighting(np.zeros((2, 2)), 1000.0, "Wk") 143 + with pytest.raises(ValueError, match="positive"): 144 + hv.apply_weighting([1.0, 2.0], 0.0, "Wk") 145 + 146 + 147 + # --------------------------------------------------------------------------- 148 + # Band method a_w (ISO 2631-1 Eq. (9) / ISO 5349-1 Eq. (A.1)). 149 + # --------------------------------------------------------------------------- 150 + def test_weighted_acceleration_matches_manual_sum() -> None: 151 + freqs = np.array([8.0, 16.0, 31.5, 63.0]) 152 + accel = np.array([0.5, 0.8, 0.3, 0.1]) 153 + result = hv.weighted_acceleration(accel, freqs, "Wk") 154 + factors = hv.weighting_factors("Wk", freqs) 155 + expected = math.sqrt(np.sum((factors * accel) ** 2)) 156 + assert result.overall == pytest.approx(expected) 157 + assert np.allclose(result.weighted, factors * accel) 158 + assert result.weighting_name == "Wk" 159 + 160 + 161 + def test_weighted_acceleration_single_band_equals_factor_times_level() -> None: 162 + result = hv.weighted_acceleration([1.0], [10.0], "Wh") 163 + assert result.overall == pytest.approx(hv.weighting_factors("Wh", 10.0)[0]) 164 + 165 + 166 + def test_weighted_acceleration_length_mismatch_raises() -> None: 167 + with pytest.raises(ValueError, match="equal-length"): 168 + hv.weighted_acceleration([1.0, 2.0], [8.0], "Wk") 169 + 170 + 171 + # --------------------------------------------------------------------------- 172 + # Time-domain metrics. 173 + # --------------------------------------------------------------------------- 174 + def test_rms_metrics_on_pure_sine() -> None: 175 + fs = 1000.0 176 + t = np.arange(int(10 * fs)) / fs 177 + a = 2.0 178 + x = a * np.sin(2.0 * math.pi * 20.0 * t) # r.m.s. = a/sqrt(2) 179 + rms = a / math.sqrt(2.0) 180 + # Crest factor of a sine is sqrt(2) (discrete peak sampling -> ~0,2 %). 181 + assert hv.crest_factor(x) == pytest.approx(math.sqrt(2.0), rel=5e-3) 182 + # VDV = (integral a_w^4 dt)^(1/4); for a sine of duration T: 183 + # mean of sin^4 = 3/8 -> VDV = a*(3/8*T)^(1/4). 184 + duration = t[-1] + 1.0 / fs 185 + expected_vdv = a * (3.0 / 8.0 * duration) ** 0.25 186 + assert hv.vibration_dose_value(x, fs) == pytest.approx(expected_vdv, rel=1e-2) 187 + # MSDV = (integral a_w^2 dt)^(1/2) = rms * sqrt(T). 188 + assert hv.motion_sickness_dose_value(x, fs) == pytest.approx( 189 + rms * math.sqrt(duration), rel=1e-2 190 + ) 191 + 192 + 193 + def test_running_rms_of_constant_power_signal() -> None: 194 + fs = 500.0 195 + # Steady sine -> running r.m.s. settles to the signal r.m.s. 196 + t = np.arange(int(5 * fs)) / fs 197 + x = math.sqrt(2.0) * np.sin(2.0 * math.pi * 15.0 * t) 198 + for method in ("linear", "exponential"): 199 + r = hv.running_rms(x, fs, integration_time=1.0, method=method) 200 + assert r.shape == x.shape 201 + assert float(r[-1]) == pytest.approx(1.0, rel=5e-2) 202 + assert hv.mtvv(x, fs) == pytest.approx(float(np.max(hv.running_rms(x, fs)))) 203 + 204 + 205 + def test_running_rms_validation() -> None: 206 + with pytest.raises(ValueError, match="method"): 207 + hv.running_rms([1.0, 2.0], 100.0, method="bogus") 208 + with pytest.raises(ValueError, match="integration_time"): 209 + hv.running_rms([1.0, 2.0], 100.0, integration_time=0.0) 210 + 211 + 212 + def test_crest_factor_zero_signal() -> None: 213 + assert hv.crest_factor(np.zeros(10)) == 0.0 214 + 215 + 216 + def test_crest_factor_warns_above_nine() -> None: 217 + # A lone spike gives crest = 10 / (10/sqrt(200)) = sqrt(200) ~ 14.1 > 9. 218 + x = np.zeros(200) 219 + x[0] = 10.0 220 + with pytest.warns(hv.HumanVibrationWarning, match="exceeds 9"): 221 + cf = hv.crest_factor(x) 222 + assert cf > 9.0 223 + 224 + 225 + # --------------------------------------------------------------------------- 226 + # Vector sum and daily exposure (ISO 2631-1 Eq. (10); ISO 5349-1/-2). 227 + # --------------------------------------------------------------------------- 228 + def test_vibration_total_value_unweighted() -> None: 229 + # ISO 5349-1 Eq. (1): a_hv = sqrt(x^2+y^2+z^2). 230 + assert hv.vibration_total_value([3.0, 4.0, 0.0]) == pytest.approx(5.0) 231 + 232 + 233 + def test_vibration_total_value_with_k_factors() -> None: 234 + # ISO 2631-1 7.2.3 health seated: k = 1,4 / 1,4 / 1,0. 235 + got = hv.vibration_total_value([0.5, 0.3, 0.8], k=[1.4, 1.4, 1.0]) 236 + expected = math.sqrt((1.4 * 0.5) ** 2 + (1.4 * 0.3) ** 2 + (1.0 * 0.8) ** 2) 237 + assert got == pytest.approx(expected) 238 + 239 + 240 + def test_iso5349_2_example_e2_1_single_tool() -> None: 241 + """ISO 5349-2 E.2.1: a_hv=7,4 m/s2, T=2,5 h -> A(8)=4,1 m/s2.""" 242 + a8 = hv.daily_exposure(7.4, 2.5 * 3600.0) 243 + assert a8 == pytest.approx(4.1, abs=0.05) 244 + 245 + 246 + def test_iso5349_2_example_e2_4_burst() -> None: 247 + """ISO 5349-2 E.2.4: a_hv=14,6 m/s2, T=4000 s (= 1,1 h) -> A(8)=5,4 m/s2.""" 248 + # T = (1000 nuts/day / 5 nuts measured) * 20 s = 4000 s per Eq. (E.4); the 249 + # standard rounds 4000 s to "1,1 h". A(8) = 14,6*sqrt(4000/28800) = 5,4. 250 + a8 = hv.daily_exposure(14.6, 4000.0) 251 + assert a8 == pytest.approx(5.4, abs=0.05) 252 + 253 + 254 + def test_iso5349_2_example_e3_forestry_multi_tool() -> None: 255 + """ISO 5349-2 E.3: three tasks combine to A(8)=3,6 m/s2.""" 256 + partials = [ 257 + hv.partial_exposure(4.6, 2 * 3600.0), # brush-saw -> 2,3 258 + hv.partial_exposure(6.0, 1 * 3600.0), # felling -> 2,1 259 + hv.partial_exposure(3.6, 2 * 3600.0), # stripping -> 1,8 260 + ] 261 + assert partials[0] == pytest.approx(2.3, abs=0.05) 262 + assert partials[1] == pytest.approx(2.1, abs=0.05) 263 + assert partials[2] == pytest.approx(1.8, abs=0.05) 264 + assert hv.combine_partial_exposures(partials) == pytest.approx(3.6, abs=0.05) 265 + 266 + 267 + def test_iso5349_1_example_multi_operation() -> None: 268 + """ISO 5349-1 5.3 worked example: 1 h/3 h/0,5 h -> A(8)=3,4 m/s2.""" 269 + a8 = hv.hav_daily_exposure( 270 + [2.0, 3.5, 10.0], [3600.0, 3 * 3600.0, 0.5 * 3600.0] 271 + ) 272 + assert a8 == pytest.approx(3.4, abs=0.05) 273 + 274 + 275 + def test_hav_daily_exposure_matches_partial_combination() -> None: 276 + values = [3.0, 5.0] 277 + durations = [2 * 3600.0, 1 * 3600.0] 278 + direct = hv.hav_daily_exposure(values, durations) 279 + partials = [hv.partial_exposure(v, d) for v, d in zip(values, durations)] 280 + assert direct == pytest.approx(hv.combine_partial_exposures(partials)) 281 + 282 + 283 + def test_energy_equivalent_acceleration() -> None: 284 + # ISO 2631-1 Eq. (B.3). 285 + got = hv.energy_equivalent_acceleration([1.0, 2.0], [1.0, 3.0]) 286 + assert got == pytest.approx(math.sqrt((1 + 4 * 3) / 4)) 287 + 288 + 289 + def test_energy_equivalent_rejects_negative_duration() -> None: 290 + with pytest.raises(ValueError, match="non-negative"): 291 + hv.energy_equivalent_acceleration([1.0, 2.0], [1.0, -1.0]) 292 + 293 + 294 + def test_hav_vwf_lifetime_matches_table_c1() -> None: 295 + """ISO 5349-1 Table C.1 / Eq. (C.1): A(8)=7 -> ~4 years; 3,7 -> ~8.""" 296 + assert hv.hav_vwf_lifetime_years(7.0) == pytest.approx(4.0, abs=0.3) 297 + assert hv.hav_vwf_lifetime_years(3.7) == pytest.approx(8.0, abs=0.6) 298 + assert hv.hav_vwf_lifetime_years(14.0) == pytest.approx(2.0, abs=0.2) 299 + 300 + 301 + # --------------------------------------------------------------------------- 302 + # Directive 2002/44/EC assessment. 303 + # --------------------------------------------------------------------------- 304 + def test_exposure_assessment_hav_zones() -> None: 305 + below = hv.exposure_assessment(2.0, kind="hav") 306 + assert below.zone == "below action" and not below.exceeds_action 307 + action = hv.exposure_assessment(3.0, kind="hav") 308 + assert action.zone == "action" and action.exceeds_action 309 + assert not action.exceeds_limit 310 + limit = hv.exposure_assessment(5.5, kind="hav") 311 + assert limit.zone == "limit" and limit.exceeds_limit 312 + assert action.action_value == hv.HAV_EAV_A8 313 + assert action.limit_value == hv.HAV_ELV_A8 314 + 315 + 316 + def test_exposure_assessment_wbv_a8_and_vdv() -> None: 317 + a = hv.exposure_assessment(0.6, kind="wbv") 318 + assert a.action_value == hv.WBV_EAV_A8 and a.limit_value == hv.WBV_ELV_A8 319 + assert a.zone == "action" 320 + v = hv.exposure_assessment(22.0, kind="wbv", metric="vdv") 321 + assert v.action_value == hv.WBV_EAV_VDV and v.limit_value == hv.WBV_ELV_VDV 322 + assert v.zone == "limit" 323 + 324 + 325 + def test_exposure_assessment_invalid() -> None: 326 + with pytest.raises(ValueError, match="kind"): 327 + hv.exposure_assessment(1.0, kind="hav", metric="vdv") 328 + with pytest.raises(ValueError, match="non-negative"): 329 + hv.exposure_assessment(-1.0, kind="hav") 330 + 331 + 332 + def test_daily_vibration_exposure_result() -> None: 333 + result = hv.daily_vibration_exposure( 334 + [4.6, 6.0, 3.6], 335 + [2 * 3600.0, 1 * 3600.0, 2 * 3600.0], 336 + kind="hav", 337 + labels=["brush-saw", "felling", "stripping"], 338 + ) 339 + assert result.a8 == pytest.approx(3.6, abs=0.05) 340 + assert result.labels == ("brush-saw", "felling", "stripping") 341 + assert result.assessment.zone == "action" 342 + assert result.partials.shape == (3,) 343 + assert result.plot # method exists 344 + 345 + 346 + def test_daily_vibration_exposure_default_labels_and_mismatch() -> None: 347 + result = hv.daily_vibration_exposure([2.0], [3600.0], kind="wbv") 348 + assert result.labels == ("op 1",) 349 + with pytest.raises(ValueError, match="labels"): 350 + hv.daily_vibration_exposure([2.0, 3.0], [1.0, 1.0], kind="hav", 351 + labels=["only-one"]) 352 + 353 + 354 + def test_weighting_names_complete() -> None: 355 + assert set(hv.WEIGHTING_NAMES) == { 356 + "Wb", "Wc", "Wd", "We", "Wf", "Wh", "Wj", "Wk", "Wm" 357 + } 358 + 359 + 360 + # --------------------------------------------------------------------------- 361 + # Remaining validation branches and the .plot() renderers. 362 + # --------------------------------------------------------------------------- 363 + def test_all_nonpositive_frequencies_return_zero_response() -> None: 364 + resp = hv.frequency_weighting("Wk", [0.0, -1.0]) 365 + assert np.all(resp.magnitude == 0.0) 366 + 367 + 368 + def test_apply_weighting_rejects_empty_signal() -> None: 369 + with pytest.raises(ValueError, match="non-empty"): 370 + hv.apply_weighting([], 1000.0, "Wk") 371 + 372 + 373 + def test_running_rms_rejects_empty_signal() -> None: 374 + with pytest.raises(ValueError, match="non-empty"): 375 + hv.running_rms([], 1000.0) 376 + 377 + 378 + def test_vibration_total_value_validates() -> None: 379 + with pytest.raises(ValueError, match="non-empty"): 380 + hv.vibration_total_value([]) 381 + with pytest.raises(ValueError, match="same length"): 382 + hv.vibration_total_value([1.0, 2.0], k=[1.4]) 383 + 384 + 385 + def test_daily_exposure_rejects_negative() -> None: 386 + with pytest.raises(ValueError, match="non-negative"): 387 + hv.daily_exposure(-1.0, 3600.0) 388 + 389 + 390 + def test_combine_partial_exposures_rejects_empty() -> None: 391 + with pytest.raises(ValueError, match="non-empty"): 392 + hv.combine_partial_exposures([]) 393 + 394 + 395 + def test_hav_daily_exposure_validates() -> None: 396 + with pytest.raises(ValueError, match="equal-length"): 397 + hv.hav_daily_exposure([2.0, 3.0], [3600.0]) 398 + with pytest.raises(ValueError, match="non-negative"): 399 + hv.hav_daily_exposure([2.0, -3.0], [3600.0, 3600.0]) 400 + 401 + 402 + def test_energy_equivalent_rejects_length_mismatch() -> None: 403 + with pytest.raises(ValueError, match="equal-length"): 404 + hv.energy_equivalent_acceleration([1.0, 2.0], [3600.0]) 405 + 406 + 407 + def test_energy_equivalent_rejects_zero_total_duration() -> None: 408 + with pytest.raises(ValueError, match="total duration"): 409 + hv.energy_equivalent_acceleration([1.0, 2.0], [0.0, 0.0]) 410 + 411 + 412 + def test_hav_vwf_lifetime_rejects_nonpositive() -> None: 413 + with pytest.raises(ValueError, match="positive"): 414 + hv.hav_vwf_lifetime_years(0.0) 415 + 416 + 417 + def test_daily_vibration_exposure_rejects_length_mismatch() -> None: 418 + with pytest.raises(ValueError, match="equal-length"): 419 + hv.daily_vibration_exposure([2.0, 3.0], [3600.0], kind="hav") 420 + 421 + 422 + def test_weighting_response_plot_returns_axes() -> None: 423 + import matplotlib 424 + 425 + matplotlib.use("Agg") 426 + import matplotlib.pyplot as plt 427 + 428 + ax = hv.frequency_weighting("Wk", [1.0, 10.0, 100.0]).plot() 429 + assert isinstance(ax, plt.Axes) 430 + plt.close("all") 431 + 432 + 433 + def test_weighted_spectrum_plot_returns_axes() -> None: 434 + import matplotlib 435 + 436 + matplotlib.use("Agg") 437 + import matplotlib.pyplot as plt 438 + 439 + result = hv.weighted_acceleration( 440 + [0.5, 0.8, 0.3], [16.0, 31.5, 63.0], "Wk" 441 + ) 442 + ax = result.plot() 443 + assert isinstance(ax, plt.Axes) 444 + plt.close("all") 445 + 446 + 447 + def test_daily_vibration_exposure_plot_returns_axes() -> None: 448 + import matplotlib 449 + 450 + matplotlib.use("Agg") 451 + import matplotlib.pyplot as plt 452 + 453 + result = hv.daily_vibration_exposure( 454 + [2.5, 3.0], [3600.0, 1800.0], kind="hav" 455 + ) 456 + ax = result.plot() 457 + assert isinstance(ax, plt.Axes) 458 + plt.close("all")