Review index · Download raw source
reference/cloud/timing_metrics.py
1"""Time-domain postprocessing for CloudScat impulse responses (NumPy only).23Units and conventions4---------------------5``dt`` is in seconds; ``h`` is a causal sampled response density, so its mass6is ``h * dt``. Array element n is assigned delay n*dt. Supply *excess arrival7delay*, after subtracting the direct source-to-observer travel time. If a8histogram uses bin centres, account for that origin before calling this module.9No missing/truncated part of h is reconstructed. Timing quantiles are10conditional on the response supplied; integrated response weight is retained.1112The source is exp(-t/3 ms)*(1-exp(-t/tau_r))*H(t), normalized to unit peak.13tau_r is calibrated to the requested rising-edge 10--90% time, relative to14the source's peak. This is deliberately distinct from onset-to-peak time.15The optional electronic response is an ideal first-order low-pass, not a16complete antialias filter. Its discrete kernel is exact for a linearly17interpolated input, preserving grid times without a half-bin leading bias.18"""1920from functools import lru_cache21import json22import math23import numpy as np242526SOURCE_DECAY_S = 3.0e-327SOURCE_CALIBRATION_DT_S = 0.2e-6282930def _validate_dt(dt):31 dt = float(dt)32 if not np.isfinite(dt) or dt <= 0:33 raise ValueError("dt must be finite and positive, in seconds")34 return dt353637def _response(h, dt):38 dt = _validate_dt(dt)39 h = np.asarray(h, dtype=float)40 if h.ndim != 1 or h.size == 0 or not np.all(np.isfinite(h)):41 raise ValueError("h must be a nonempty, finite one-dimensional array")42 scale = float(np.max(np.abs(h)))43 if np.any(h < -max(scale * 1e-12, 1e-300)):44 raise ValueError("h must be nonnegative")45 return np.maximum(h, 0.0), dt464748def _crossing(t, y, level, stop_index):49 """Earliest rising crossing before the selected peak, interpolated."""50 indexes = np.flatnonzero(y[:stop_index + 1] >= level)51 if indexes.size == 0:52 return None53 j = int(indexes[0])54 if j == 0:55 # Source and electronic outputs start at zero; this branch also56 # handles an impulse/step waveform whose first point is nonzero.57 return float(t[0])58 y0, y1 = float(y[j - 1]), float(y[j])59 if y1 == y0:60 return float(t[j])61 return float(t[j - 1] + (level - y0) / (y1 - y0) * (t[j] - t[j - 1]))626364def waveform_metrics(t_s, y, onset_s=0.0):65 """Return rising-edge metrics; onset is supplied, never noise-thresholded."""66 t_s = np.asarray(t_s, dtype=float)67 y = np.asarray(y, dtype=float)68 if t_s.ndim != 1 or t_s.size < 2 or y.shape != t_s.shape:69 raise ValueError("t_s and y must be matching one-dimensional arrays")70 if not np.all(np.isfinite(y)) or not np.all(np.isfinite(t_s)):71 raise ValueError("waveform values must be finite")72 if not np.all(np.diff(t_s) > 0):73 raise ValueError("waveform times must be strictly increasing")74 peak_index = int(np.argmax(y))75 peak = float(y[peak_index])76 if peak <= 0:77 return {"rise_10_90_us": None, "onset_to_peak_us": None,78 "t10_us": None, "t90_us": None, "peak_time_us": None,79 "peak_amplitude": 0.0, "peak_at_window_end": False}80 t10 = _crossing(t_s, y, 0.1 * peak, peak_index)81 t90 = _crossing(t_s, y, 0.9 * peak, peak_index)82 return {83 "rise_10_90_us": None if t10 is None or t90 is None else (t90 - t10) * 1e6,84 "onset_to_peak_us": float((t_s[peak_index] - onset_s) * 1e6),85 "t10_us": None if t10 is None else (t10 - onset_s) * 1e6,86 "t90_us": None if t90 is None else (t90 - onset_s) * 1e6,87 "peak_time_us": float((t_s[peak_index] - onset_s) * 1e6),88 "peak_amplitude": peak,89 "peak_at_window_end": bool(peak_index == y.size - 1),90 }919293def _source_peak(tau_r_s, decay_s):94 t_peak = tau_r_s * math.log1p(decay_s / tau_r_s)95 peak = math.exp(-t_peak / decay_s) * (-math.expm1(-t_peak / tau_r_s))96 return t_peak, peak979899def _source_rise_on_fine_grid(tau_r_s, decay_s):100 t_peak, peak = _source_peak(tau_r_s, decay_s)101 t = np.arange(0.0, t_peak, SOURCE_CALIBRATION_DT_S)102 t = np.append(t, t_peak)103 y = np.exp(-t / decay_s) * (-np.expm1(-t / tau_r_s)) / peak104 t10 = _crossing(t, y, 0.1, y.size - 1)105 t90 = _crossing(t, y, 0.9, y.size - 1)106 return t90 - t10107108109@lru_cache(maxsize=128)110def source_parameters(intrinsic_rise_us, decay_s=SOURCE_DECAY_S):111 """Calibrate tau_r by bisection with 0.2 us crossing interpolation."""112 target = float(intrinsic_rise_us) * 1e-6113 decay_s = float(decay_s)114 if not np.isfinite(target) or target <= 0 or decay_s <= 0:115 raise ValueError("source rise and decay time must be positive")116 lo = max(target * 1e-4, 1e-12)117 hi = max(target, 1e-9)118 for _ in range(30):119 if _source_rise_on_fine_grid(hi, decay_s) >= target:120 break121 hi *= 2.0122 else:123 raise ValueError("requested rise is too long for the selected decay time")124 for _ in range(64):125 mid = (lo + hi) / 2.0126 if _source_rise_on_fine_grid(mid, decay_s) < target:127 lo = mid128 else:129 hi = mid130 tau_r = (lo + hi) / 2.0131 t_peak, _ = _source_peak(tau_r, decay_s)132 return {"intrinsic_rise_us": float(intrinsic_rise_us),133 "tau_r_s": tau_r, "tau_r_us": tau_r * 1e6,134 "decay_s": decay_s, "decay_ms": decay_s * 1e3,135 "analytic_onset_to_peak_us": t_peak * 1e6,136 "calibration_dt_us": SOURCE_CALIBRATION_DT_S * 1e6}137138139def pulse_waveform(dt, intrinsic_rise_us, decay_s=SOURCE_DECAY_S,140 tail_e_folds=12.0):141 """Return (time_seconds, unit-peak source); tail is finite and documented."""142 dt = _validate_dt(dt)143 p = source_parameters(float(intrinsic_rise_us), float(decay_s))144 t = np.arange(int(math.ceil(tail_e_folds * decay_s / dt)) + 1) * dt145 _, peak = _source_peak(p["tau_r_s"], decay_s)146 y = np.exp(-t / decay_s) * (-np.expm1(-t / p["tau_r_s"])) / peak147 return t, y148149150def _fft_convolve(x, z):151 size = len(x) + len(z) - 1152 nfft = 1 << (size - 1).bit_length()153 y = np.fft.irfft(np.fft.rfft(x, nfft) * np.fft.rfft(z, nfft), nfft)[:size]154 return np.maximum(y, 0.0)155156157def one_pole_kernel(dt, fc_hz, tail_tolerance=1e-12):158 """Exact first-order kernel for piecewise-linear input sample values.159160 y[n] = a*y[n-1] + b0*x[n] + b1*x[n-1], zero prehistory.161 This is not a claim that a first-order filter provides adequate antialiasing.162 """163 dt = _validate_dt(dt)164 fc_hz = float(fc_hz)165 if not np.isfinite(fc_hz) or fc_hz <= 0:166 raise ValueError("fc_hz must be positive, or None to omit electronics")167 u = 2.0 * math.pi * fc_hz * dt168 a = math.exp(-u)169 q = -math.expm1(-u) / u170 b0, b1 = 1.0 - q, q - a171 n = max(2, int(math.ceil(-math.log(tail_tolerance) / u)) + 2)172 kernel = np.empty(n)173 kernel[0] = b0174 kernel[1:] = (a * b0 + b1) * np.exp(-u * np.arange(n - 1))175 return kernel176177178def impulse_metrics(h, dt):179 """Timing of the supplied weighted impulse, with total mass unnormalized."""180 h, dt = _response(h, dt)181 mass = h * dt182 total = float(np.sum(mass))183 out = {"integrated_response_weight": total,184 "response_window_end_us": (h.size - 1) * dt * 1e6,185 "quantiles_conditional_on_supplied_response": True,186 "delay_grid_origin_us": 0.0,187 "delay_grid_convention": "mass at delay n*dt"}188 if total == 0:189 out.update({"status": "no_detected_weight", "q10_us": None,190 "q50_us": None, "q90_us": None,191 "spread_q90_q10_us": None,192 "fraction_after_100us": None, "fraction_after_200us": None,193 "fraction_after_3ms": None})194 return out195 t = np.arange(h.size) * dt196 positive = mass > 0197 # Invert the weighted empirical CDF without interpolating across gaps:198 # a direct-light atom plus a separate delayed atom must retain that gap.199 cdf = np.cumsum(mass)[positive] / total200 times = t[positive]201 indexes = np.searchsorted(cdf, [0.1, 0.5, 0.9], side="left")202 quantiles = times[np.minimum(indexes, times.size - 1)] * 1e6203 out.update({"status": "ok", "q10_us": float(quantiles[0]),204 "q50_us": float(quantiles[1]), "q90_us": float(quantiles[2]),205 "spread_q90_q10_us": float(quantiles[2] - quantiles[0])})206 for name, limit in [("100us", 100e-6), ("200us", 200e-6), ("3ms", 3e-3)]:207 out["fraction_after_" + name] = float(np.sum(mass[t > limit]) / total)208 return out209210211def propagate_waveform(h, dt, intrinsic_rise_us, fc_hz=None, return_metadata=False):212 """Return (time_seconds, output); output retains response weight.213214 With return_metadata=True, return (time_seconds, output, metadata).215 """216 h, dt = _response(h, dt)217 _, source = pulse_waveform(dt, intrinsic_rise_us)218 y = _fft_convolve(source, h * dt)219 if fc_hz is not None:220 y = _fft_convolve(y, one_pole_kernel(dt, fc_hz))221 t = np.arange(y.size) * dt222 metadata = {"time_step_us": dt * 1e6,223 "intrinsic_rise_us": float(intrinsic_rise_us),224 "source": dict(source_parameters(float(intrinsic_rise_us))),225 "analog_bandwidth_hz": None if fc_hz is None else float(fc_hz),226 "source_tail_e_folds": 12.0,227 "electronic_kernel": None if fc_hz is None else "exact piecewise-linear input",228 "output_window_end_ms": float(t[-1] * 1e3)}229 return (t, y, metadata) if return_metadata else (t, y)230231232def metrics(h, dt, intrinsic_rise_us, fc_hz=None):233 """JSON-serializable result; use propagate_waveform for waveform arrays."""234 h, dt = _response(h, dt)235 im = impulse_metrics(h, dt)236 if im["integrated_response_weight"] == 0:237 return {"status": "no_detected_weight", "intrinsic_rise_us": float(intrinsic_rise_us),238 "analog_bandwidth_hz": None if fc_hz is None else float(fc_hz),239 "rise_10_90_us": None, "onset_to_peak_us": None, "impulse": im}240 t, y, result = propagate_waveform(h, dt, intrinsic_rise_us, fc_hz, return_metadata=True)241 result.update(waveform_metrics(t, y))242 result.update({"status": "ok", "impulse": im})243 return result244245246def _summary(values):247 a = np.asarray([v for v in values if v is not None], dtype=float)248 if a.size == 0:249 return {k: None for k in ["min", "p05", "median", "p95", "max"]}250 v = np.quantile(a, [0.0, 0.05, 0.5, 0.95, 1.0])251 return dict(zip(["min", "p05", "median", "p95", "max"], map(float, v)))252253254def sampling_phase_metrics(y, dt, fs_hz, phase_count=128, fc_hz=None):255 """Noiseless point-sampling across clock phases, using sampled peak.256257 Each sampled waveform's own peak defines its 10/90% levels. Negative-time258 samples are explicitly zero. Crossing times use linear interpolation.259 These spreads do not include shot noise or detector/ADC noise, and the260 optional fc check is only a necessary Nyquist check, not filter validation.261 """262 dt = _validate_dt(dt)263 y = np.asarray(y, dtype=float)264 if y.ndim != 1 or y.size < 2 or not np.all(np.isfinite(y)):265 raise ValueError("y must be a finite one-dimensional waveform")266 fs_hz = float(fs_hz)267 if not np.isfinite(fs_hz) or fs_hz <= 0 or phase_count < 2:268 raise ValueError("fs_hz must be positive and phase_count at least two")269 period = 1.0 / fs_hz270 t = np.arange(y.size) * dt271 rises, peaks = [], []272 for phase in np.arange(int(phase_count)) * period / phase_count:273 n = int(math.floor((t[-1] - phase) / period))274 sampled_t = phase + np.arange(-1, n + 1) * period275 sampled_y = np.interp(sampled_t, t, y, left=0.0, right=0.0)276 m = waveform_metrics(sampled_t, sampled_y)277 rises.append(m["rise_10_90_us"])278 peaks.append(m["onset_to_peak_us"])279 continuous = waveform_metrics(t, y)280 flags = []281 if fs_hz <= 1.0e4:282 resolution = "invalid_resolution_for_30us_rise"283 flags.append("100us or longer sample intervals cannot resolve a 30us rise")284 elif fs_hz < 1.0e5:285 resolution = "marginal_resolution_for_30us_rise"286 flags.append("fewer than three sample intervals across a 30us rise")287 else:288 resolution = "at_least_three_intervals_per_30us_rise"289 if period < dt:290 flags.append("sampling faster than computational grid: interpolation adds no information")291 if fc_hz is not None and float(fc_hz) >= 0.5 * fs_hz:292 flags.append("analog cutoff reaches/exceeds Nyquist; sampling result is not antialiased")293 return {"sampling_rate_hz": fs_hz, "sample_interval_us": period * 1e6,294 "phase_count": int(phase_count), "resolution_status": resolution,295 "flags": flags, "noise_included": False,296 "continuous_rise_10_90_us": continuous["rise_10_90_us"],297 "sampled_rise_10_90_us": _summary(rises),298 "sampled_onset_to_peak_us": _summary(peaks),299 "peak_reference": "each sampled waveform's own peak"}300301302def selfcheck():303 """Compare a delta cloud response with independent continuous expressions."""304 dt = 2e-6305 delta = np.array([1.0 / dt])306 # Independent continuous convolution of two source exponentials with RC.307 expected_rises = {(30, 10000): 48.03517444, (30, 50000): 31.31325808,308 (80, 10000): 90.28856458, (80, 50000): 80.40709689,309 (200, 10000): 204.42401224, (200, 50000): 200.03885856}310 rows = []311 for (rise, fc), expected in expected_rises.items():312 m = metrics(delta, dt, rise, fc)313 error = abs(m["rise_10_90_us"] - expected)314 if error > 0.6:315 raise AssertionError("RC/source analytical check failed: " + str((rise, fc, error)))316 rows.append({"intrinsic_rise_us": rise, "fc_hz": fc,317 "computed_rise_us": m["rise_10_90_us"],318 "continuous_reference_rise_us": expected,319 "absolute_error_us": error,320 "computed_onset_to_peak_us": m["onset_to_peak_us"]})321 im = impulse_metrics(delta, dt)322 assert im["integrated_response_weight"] == 1.0 and im["q90_us"] == 0.0323 t, y = propagate_waveform(delta, dt, 30, 10000)324 sampled = sampling_phase_metrics(y, dt, 10000, fc_hz=10000)325 assert sampled["resolution_status"] == "invalid_resolution_for_30us_rise"326 assert metrics(np.zeros(10), dt, 30)["status"] == "no_detected_weight"327 return {"status": "passed", "dt_us": dt * 1e6, "analytic_checks": rows,328 "example_sampling_10kSps": sampled}329330331if __name__ == "__main__":332 print(json.dumps(selfcheck(), indent=2))333