diff --git a/pysp2/util/__init__.py b/pysp2/util/__init__.py
index 68bba37..0a224ee 100644
--- a/pysp2/util/__init__.py
+++ b/pysp2/util/__init__.py
@@ -18,3 +18,4 @@
 from .particle_properties import calc_diams_masses, process_psds
 from .deadtime import deadtime
 from .leo_fit import beam_shape,leo_fit
+from .normalized_derivative_method import central_difference
diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py
new file mode 100644
index 0000000..b89ae66
--- /dev/null
+++ b/pysp2/util/normalized_derivative_method.py
@@ -0,0 +1,66 @@
+import numpy as np
+from scipy.optimize import curve_fit
+import xarray as xr
+
+def central_difference(S, num_records=None):
+
+    """ 
+    Compute normalized derivative S'(t) / S(t) using the 
+    central difference scheme (Moteki & Kondo, Eq. A.2). 
+    Interior points:
+    - Fourth-order central difference 
+    Edge cases: 
+    - First two points: fourth-order forward difference  
+    - Last two points: fourth-order backward difference 
+
+    Parameters 
+    ---------- 
+    S: xarray Dataset 
+    The scattering signal dataset. 
+    num_records: int or None 
+    Only process first num_records datapoints. 
+    Set to None to process all records. 
+
+    Returns 
+    ------- 
+    dSdt : xarray Dataset 
+    Fourth-order numerical derivative S'(t). 
+    """
+
+    if num_records is None:
+        num_records = S.sizes['event_index']
+
+    dt = 200e-9
+    channels = ['Data_ch0', 'Data_ch4']
+    dSdt = {}
+
+    for ch in channels:
+        y = S[ch].isel(event_index=slice(0, num_records)).values
+        d = np.full_like(y, np.nan, dtype=np.float64)
+
+        # Interior points (vectorized)
+        d[:, 2:-2] = (
+            -y[:, 4:] + 8*y[:, 3:-1] - 8*y[:, 1:-3] + y[:, 0:-4]
+        ) / (12 * dt)
+
+        # Forward edges
+        d[:, 0] = (-25*y[:, 0] + 48*y[:, 1] - 36*y[:, 2] + 16*y[:, 3] - 3*y[:, 4]) / (12*dt)
+        d[:, 1] = (-25*y[:, 1] + 48*y[:, 2] - 36*y[:, 3] + 16*y[:, 4] - 3*y[:, 5]) / (12*dt)
+
+        # Backward edges
+        n = y.shape[1]
+        d[:, n-2] = (25*y[:, n-2] - 48*y[:, n-3] + 36*y[:, n-4] - 16*y[:, n-5] + 3*y[:, n-6]) / (12*dt)
+        d[:, n-1] = (25*y[:, n-1] - 48*y[:, n-2] + 36*y[:, n-3] - 16*y[:, n-4] + 3*y[:, n-5]) / (12*dt)
+
+        dSdt[ch] = xr.DataArray(
+            d,
+            dims=('event_index', 'time'),
+            coords={
+                'event_index': S['event_index'][:num_records],
+                'time': S['time']
+            },
+            name=f'd{ch}_dt'
+        )
+
+    return xr.Dataset(dSdt)
+
diff --git a/tests/test_ndm.py b/tests/test_ndm.py
new file mode 100644
index 0000000..8867890
--- /dev/null
+++ b/tests/test_ndm.py
@@ -0,0 +1,16 @@
+import pysp2
+import numpy as np
+
+def test_central_difference():
+    my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B)
+    my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI)
+    my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False)
+    dSdt = pysp2.util.central_difference(my_binary)
+    
+    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=0).item(), 
+                                   8.3333333333e6, decimal=2)
+    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=99).item(), 
+                                   7.166666666e7, decimal=2)
+    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=19).item(), 
+                                   1.5e7, decimal=2)
+    assert np.isfinite(dSdt).all()
\ No newline at end of file