diff --git a/src/mrinufft/io/nsp.py b/src/mrinufft/io/nsp.py
index e1a196c7..ef8518d1 100644
--- a/src/mrinufft/io/nsp.py
+++ b/src/mrinufft/io/nsp.py
@@ -404,6 +404,7 @@ def read_trajectory(
 def read_arbgrad_rawdat(
     filename: str,
     removeOS: bool = False,
+    doAverage: bool = True,
     squeeze: bool = True,
     slice_num: int | None = None,
     contrast_num: int | None = None,
@@ -417,6 +418,8 @@ def read_arbgrad_rawdat(
         The path to the Siemens MRI file.
     removeOS : bool, optional
         Whether to remove the oversampling, by default False.
+    doAverage : bool, optional
+        Whether to average the data acquired along NAve dimension, by default True.
     squeeze : bool, optional
         Whether to squeeze the dimensions of the data, by default True.
     slice_num : int, optional
@@ -447,6 +450,7 @@ def read_arbgrad_rawdat(
     data, hdr, twixObj = read_siemens_rawdat(
         filename=filename,
         removeOS=removeOS,
+        doAverage=doAverage,
         squeeze=squeeze,
         slice_num=slice_num,
         contrast_num=contrast_num,
diff --git a/src/mrinufft/io/siemens.py b/src/mrinufft/io/siemens.py
index 325fcbf0..dd8c96d7 100644
--- a/src/mrinufft/io/siemens.py
+++ b/src/mrinufft/io/siemens.py
@@ -8,6 +8,7 @@
 def read_siemens_rawdat(
     filename: str,
     removeOS: bool = False,
+    doAverage: bool = True,
     squeeze: bool = True,
     return_twix: bool = True,
     slice_num: int | None = None,
@@ -21,6 +22,8 @@ def read_siemens_rawdat(
         The path to the Siemens MRI file.
     removeOS : bool, optional
         Whether to remove the oversampling, by default False.
+    doAverage : bool, option
+        Whether to average the data acquired along NAve dimension.
     squeeze : bool, optional
         Whether to squeeze the dimensions of the data, by default True.
     data_type : str, optional
@@ -61,12 +64,14 @@ def read_siemens_rawdat(
     if isinstance(twixObj, list):
         twixObj = twixObj[-1]
     twixObj.image.flagRemoveOS = removeOS
+    twixObj.image.flagDoAverage = doAverage
     hdr = {
         "n_coils": int(twixObj.image.NCha),
-        "n_shots": int(twixObj.image.NLin),
+        "n_shots": int(twixObj.image.NLin) * int(twixObj.image.NPar),
         "n_contrasts": int(twixObj.image.NSet),
         "n_adc_samples": int(twixObj.image.NCol),
         "n_slices": int(twixObj.image.NSli),
+        "n_average": int(twixObj.image.NAve),
     }
     if slice_num is not None and hdr["n_slices"] < slice_num:
         raise ValueError("The slice number is out of bounds.")
@@ -92,6 +97,7 @@ def read_siemens_rawdat(
         hdr["n_shots"] * hdr["n_adc_samples"],
         hdr["n_slices"] if slice_num is None else 1,
         hdr["n_contrasts"] if contrast_num is None else 1,
+        hdr["n_average"] if hdr["n_average"] > 1 and not doAverage else 1,
     )
     if return_twix:
         return data, hdr, twixObj