diff --git a/crates/clawhdf5-format/src/filters.rs b/crates/clawhdf5-format/src/filters.rs
index 2e99d70..0a0225f 100644
--- a/crates/clawhdf5-format/src/filters.rs
+++ b/crates/clawhdf5-format/src/filters.rs
@@ -845,9 +845,33 @@ fn shuffle_decompress(data: &[u8], element_size: usize) -> Result<Vec<u8>, Forma
     let num_elements = data.len() / element_size;
     let mut result = vec![0u8; data.len()];
 
-    for i in 0..num_elements {
-        for j in 0..element_size {
-            result[i * element_size + j] = data[j * num_elements + i];
+    // The shuffled stream is `element_size` byte planes of `num_elements`
+    // bytes each; un-shuffling interleaves them. This is on the read path of
+    // every compressed dataset (shuffle is applied automatically before
+    // compression). The naive `result[i * es + j] = data[j * n + i]` form does
+    // a multiply and two bounds checks per byte and defeats vectorisation;
+    // fixed-width plane arrays sliced to a common length let the compiler
+    // hoist the checks and emit interleaves for the common 4- and 8-byte
+    // element sizes.
+    fn interleave<const W: usize>(data: &[u8], n: usize, out: &mut [u8]) {
+        let planes: [&[u8]; W] = core::array::from_fn(|j| &data[j * n..(j + 1) * n]);
+        for (i, element) in out.chunks_exact_mut(W).enumerate() {
+            for j in 0..W {
+                element[j] = planes[j][i];
+            }
+        }
+    }
+    match element_size {
+        2 => interleave::<2>(data, num_elements, &mut result),
+        4 => interleave::<4>(data, num_elements, &mut result),
+        8 => interleave::<8>(data, num_elements, &mut result),
+        16 => interleave::<16>(data, num_elements, &mut result),
+        _ => {
+            for (i, element) in result.chunks_exact_mut(element_size).enumerate() {
+                for (j, byte) in element.iter_mut().enumerate() {
+                    *byte = data[j * num_elements + i];
+                }
+            }
         }
     }
 
@@ -1848,4 +1872,20 @@ mod tests {
         };
         assert!(decompress_chunk(&data, &pipeline, 16, 1).is_err());
     }
+    #[test]
+    fn unshuffle_inverts_shuffle_for_every_element_size() {
+        for element_size in [1usize, 2, 3, 4, 5, 8, 12, 16, 24] {
+            for elements in [0usize, 1, 2, 7, 64, 1000] {
+                let original: Vec<u8> = (0..element_size * elements)
+                    .map(|i| (i * 31 + 7) as u8)
+                    .collect();
+                let shuffled = shuffle_compress(&original, element_size).unwrap();
+                assert_eq!(
+                    shuffle_decompress(&shuffled, element_size).unwrap(),
+                    original,
+                    "element_size {element_size}, {elements} elements"
+                );
+            }
+        }
+    }
 }