1"""Download single-channel ephys traces from the benchcompress datasets and
2cache them as .npy, together with the bandpass+requantized "filtered" variant
3that benchcompress benchmarks.
5Usage: python fetch.py [outdir]
6"""
7import os
8import sys
10import numpy as np
12CACHE = sys.argv[1] if len(sys.argv) > 1 else os.path.join(
13 os.path.dirname(os.path.abspath(__file__)), "cache")
15NUM_SAMPLES = 500_000
16RATE = 30000.0
18# (name, dandi asset url, dataset path in the nwb, channel)
19SOURCES = [
20 ("ecephys-000876-ch45",
21 "https://api.dandiarchive.org/api/assets/7e1de06d-d478-40e2-9b64-9dd04eafaa4c/download/",
22 "/acquisition/ElectricalSeriesAP/data", 45),
23 ("ecephys-000409-ch101",
24 "https://api.dandiarchive.org/api/assets/c04f6b30-82bf-40e1-9210-34f0bcd8be24/download/",
25 "/acquisition/ElectricalSeriesAp/data", 101),
26 ("ecephys-001290-ch0",
27 "https://api.dandiarchive.org/api/assets/78c99d23-da88-4ecd-9086-c488a126eac5/download/",
28 "/acquisition/ElectricalSeriesAPImec/data", 0),
29]
32def bandpass(x, lowcut, highcut, rate):
33 from scipy.signal import butter, lfilter
34 nyq = 0.5 * rate
35 b, a = butter(5, [lowcut / nyq, highcut / nyq], btype="band")
36 return lfilter(b, a, x)
39def highpass(x, lowcut, rate):
40 from scipy.signal import butter, lfilter
41 nyq = 0.5 * rate
42 b, a = butter(5, lowcut / nyq, btype="high")
43 return lfilter(b, a, x)
46def noise_level(x, rate):
47 xf = highpass(x, 300.0, rate)
48 return float(np.median(np.abs(xf - np.median(xf))) / 0.6745)
51def filtered_variant(x, rate=RATE, v=0.25, lowcut=300.0, highcut=6000.0):
52 """benchcompress's `-filtered` transform: bandpass, normalize by the MAD
53 noise level, requantize at step v (so the noise std is ~1/v = 4 steps)."""
54 xf = bandpass(x - np.median(x), lowcut, highcut, rate)
55 nl = noise_level(xf, rate)
56 return np.round(xf / nl / v).astype(np.int16)
59def main():
60 os.makedirs(CACHE, exist_ok=True)
61 for name, url, path, ch in SOURCES:
62 raw_path = os.path.join(CACHE, f"{name}.raw.npy")
63 filt_path = os.path.join(CACHE, f"{name}.filtered.npy")
64 if os.path.exists(raw_path) and os.path.exists(filt_path):
65 print(f"{name}: cached")
66 continue
67 print(f"{name}: downloading {NUM_SAMPLES} samples ...", flush=True)
68 import lindi
69 h5f = lindi.LindiH5pyFile.from_hdf5_file(url)
70 ds = h5f[path]
71 raw = np.asarray(ds[:NUM_SAMPLES, ch : ch + 1]).flatten().astype(np.int16)
72 np.save(raw_path, raw)
73 np.save(filt_path, filtered_variant(raw.astype(np.float64)))
74 print(f"{name}: raw std={raw.std():.1f} "
75 f"filtered std={np.load(filt_path).std():.2f}", flush=True)
78if __name__ == "__main__":
79 main()