diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index 94f2fb59c99ed2bd76e1860b1c6fa78bdf1f4d19..9b657eb3898c9399feb47768e181bc305bfc14d8 100644
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -57,9 +57,11 @@ class Instrument:
 
     def __init__(self,
                  # Sampling parameters
-                 size=2592000, dt=1/4, t0='orbits', telemetry_dt=86400,
+                 size=2592000, dt=1/4, t0='orbits',
                  # Physics simulation sampling and filtering
                  physics_upsampling=4, aafilter=('kaiser', 240, 1.1, 2.9),
+                 # Telemetry sampling
+                 telemetry_downsampling=86400 * 4, telemetry_t0='orbits',
                  # Inter-spacecraft propagation
                  orbits='static', neglect_tps=False,
                  gws=None, interpolation=('lagrange', 31),
@@ -98,6 +100,8 @@ class Instrument:
                 or None for no filter; to design a filter from a Kaiser window, use a tuple
                 ('kaiser', attenuation [dB], f1 [Hz], f2 [Hz]) with f1 < f2 the frequencies defining
                 the transition band
+            telemetry_downsampling: ratio of sampling frequencies for measurements vs. telemetry events
+            telemetry_t0: time of initial telemetry event [s], or 'orbits' to match that of the orbits
             orbits: path to orbit file, dictionary of constant PPRs for static arms, 'static'
                 for a set of static PPRs corresponding to a fit of Keplerian orbits around t = 0,
                 or dictionary of PPR time series
@@ -157,26 +161,39 @@ class Instrument:
         self.simulated = False
 
         # Measurement sampling
-        self.size = int(size)
-        self.dt = float(dt)
         if t0 == 'orbits':
             if isinstance(orbits, str) and orbits != 'static':
                 logger.debug("Reading initial time from orbit file '%s'", orbits)
                 with h5py.File(orbits, 'r') as orbitf:
-                    self.t0 = float(orbitf.attrs['tau0'])
+                    attr = 't0' if neglect_tps else 'tau0'
+                    self.t0 = float(orbitf.attrs[attr])
             else:
                 self.t0 = 0.0
         else:
             self.t0 = float(t0)
+        self.size = int(size)
+        self.dt = float(dt)
         self.fs = 1 / self.dt
         self.duration = self.dt * self.size
         logger.info("Computing measurement time vector (size=%s, dt=%s)", self.size, self.dt)
         self.t = self.t0 + numpy.arange(self.size, dtype=numpy.float64) * self.dt
 
-        self.telemetry_dt = telemetry_dt
-        self.telemetry_fs = 1 / telemetry_dt
-        self.telemetry_size = int(numpy.ceil(size / telemetry_dt * self.dt))
-        self.telemetry_t = self.t0 + numpy.arange(self.telemetry_size, dtype=numpy.float64) * self.telemetry_dt
+        # Sampling of telemetry events given in TCB
+        if telemetry_t0 == 'orbits':
+            if isinstance(orbits, str) and orbits != 'static':
+                logger.debug("Reading telemetry initial time from orbit file '%s'", orbits)
+                with h5py.File(orbits, 'r') as orbitf:
+                    self.telemetry_t0 = float(orbitf.attrs['t0'])
+            else:
+                self.telemetry_t0 = 0.0
+        else:
+            self.telemetry_t0 = float(telemetry_t0)
+        self.telemetry_downsampling = int(telemetry_downsampling)
+        self.telemetry_dt = self.dt * self.telemetry_downsampling
+        self.telemetry_fs = 1 / self.telemetry_dt
+        self.telemetry_size = int(numpy.ceil(size / self.telemetry_dt * self.dt))
+        self.telemetry_t = \
+            self.telemetry_t0 + numpy.arange(self.telemetry_size, dtype=numpy.float64) * self.telemetry_dt
 
         # Physics sampling
         self.physics_upsampling = int(physics_upsampling)
@@ -1516,6 +1533,7 @@ class Instrument:
 
         self.write_metadata(hdf5)
         hdf5['physics_t'] = self.physics_t
+        hdf5['telemetry_t'] = self.telemetry_t
 
         if keep_all: