diff --git a/src/bids/modeling/transformations/compute.py b/src/bids/modeling/transformations/compute.py
index 8e5068b7..06cf01c2 100644
--- a/src/bids/modeling/transformations/compute.py
+++ b/src/bids/modeling/transformations/compute.py
@@ -89,13 +89,29 @@ def _transform(self, var, model='spm', derivative=False, dispersion=False,
             res=min_interval)
         # Bound the effective sampling rate between min_freq and max_freq
         effective_sr = max(min_freq, min(safety / required_resolution, max_freq))
+        # convert fir_delays from scans (BOLD sampling rate) to resample_frames sampling rate
+        if fir_delays is not None:
+            tr = var.run_info[0].tr
+            dt = resample_frames[1] - resample_frames[0]
+            fir_delays_resample_adjusted = np.round(np.array(fir_delays) * tr / dt).astype(int)
+            adjusted_to_original_mapping = dict(zip(fir_delays_resample_adjusted, fir_delays))
+        else:
+            fir_delays_resample_adjusted = None
+            adjusted_to_original_mapping = {}
+
         convolved = hrf.compute_regressor(
-            vals, model, resample_frames, fir_delays=fir_delays, min_onset=0,
+            vals, model, resample_frames, fir_delays=fir_delays_resample_adjusted, min_onset=0,
             oversampling=int(np.ceil(effective_sr / sampling_rate))
-            )
-
+            )        
         results = []
         arr, names = convolved
+
+        # If FIR model, rename columns to original delays
+        if model == 'fir' and adjusted_to_original_mapping:
+            names = [
+                f"{name.rsplit('_', 1)[0]}_{adjusted_to_original_mapping[int(name.rsplit('_', 1)[1])]}"
+                for name in names
+            ]
         for conv, name in zip(np.split(arr, arr.shape[1], axis=1), names):
             new_name = '_'.join([var.name, name.split('_')[-1]]) if '_' in name else var.name
             results.append(