88and a clinician / certified Holter software. Use this only to get a sense of the
99order of magnitude.
1010
11+ Note: the ectopic percentage is derived from R-R timing with a self-scaling
12+ detection threshold, so it is invariant to the ADC gain -- microvolt calibration
13+ does not change these numbers.
14+
1115Pipeline:
12- 1. Reconstruct one channel (CH3 -- sharpest R waves, fewest artifacts) in
13- chunks with a continuous cumulative sum.
16+ 1. Reconstruct one channel in chunks with a continuous cumulative sum.
1417 2. Band-pass 5-25 Hz (QRS enhancement; removes drift/respiration/noise).
1518 3. R-peak detection: adaptive threshold + 250 ms refractory + T-wave rejection.
16- 4. Per-beat classification against the local R-R:
17- - premature (ectopic): RR < PREMATURE * local-median RR, and the next
18- RR > COMPENSATORY * local-median (compensatory pause);
19- - pause: RR > PAUSE_S seconds.
19+ 4. Per-beat classification against the local R-R (premature + compensatory pause).
2020 5. Percentages + per-hour breakdown.
2121"""
22- import numpy as np
22+ import argparse
2323import datetime
24+ import numpy as np
2425from scipy .signal import butter , filtfilt , find_peaks
2526
26- BLK , HDR , NCH = 1024 , 13 , 3
27- SPB = (BLK - HDR ) // NCH # 337
28- FS = 300
29- CH = 2 # CH3 (index 2) -- cleanest for detection
30- START = datetime .datetime (2026 , 7 , 1 , 10 , 12 , 17 )
27+ BLOCK_SIZE = 1024
3128
3229# rhythm-classification thresholds
3330PREMATURE = 0.80 # RR shorter than 80% of local median = premature
3734TWAVE_RATIO = 0.55 # a peak below 55% of the previous one within the window = likely T
3835
3936
40- def bandpass (x , lo = 5 , hi = 25 ):
41- b , a = butter (2 , [lo / (FS / 2 ), hi / (FS / 2 )], "band" )
37+ def bandpass (x , fs , lo = 5 , hi = 25 ):
38+ b , a = butter (2 , [lo / (fs / 2 ), hi / (fs / 2 )], "band" )
4239 return filtfilt (b , a , x )
4340
4441
45- def detect_all (path = "ECG.DAT" ):
46- mm = np .fromfile (path , dtype = np .uint8 ).reshape (- 1 , BLK )
47- nblk = int (np .where (mm .any (axis = 1 ))[0 ].max ()) + 1
48- payload = mm [:nblk , HDR :HDR + SPB * NCH ].reshape (- 1 , NCH )[:, CH ].astype (np .int8 ).astype (np .int32 )
42+ def detect_all (path , fs , channels , header_len , detect_channel ):
43+ mm = np .fromfile (path , dtype = np .uint8 ).reshape (- 1 , BLOCK_SIZE )
44+ nz = np .where (mm .any (axis = 1 ))[0 ]
45+ nblk = int (nz .max ()) + 1 if nz .size else 0
46+ spb = (BLOCK_SIZE - header_len ) // channels
47+ payload = (mm [:nblk , header_len :header_len + spb * channels ]
48+ .reshape (- 1 , channels )[:, detect_channel ].astype (np .int8 ).astype (np .int32 ))
4949 total = payload .size
5050
51- CHUNK , OVER = 300_000 , 3_000 # ~1000 s chunks, overlap to keep edge beats
51+ CHUNK , OVER = 300_000 , 3_000
5252 acc , pos = 0 , 0
5353 peaks , amps = [], []
5454 while pos < total :
5555 a , b = max (0 , pos - OVER ), min (total , pos + CHUNK )
5656 base = acc - int (payload [a :pos ].sum ()) if pos > a else acc
5757 sig = np .cumsum (payload [a :b ]) + base
58- f = bandpass (sig .astype (float ))
58+ f = bandpass (sig .astype (float ), fs )
5959 thr = max (4.0 * np .median (np .abs (f )) / 0.6745 , 1.0 )
60- pk , props = find_peaks (f , height = thr , distance = int (0.25 * FS ))
60+ pk , props = find_peaks (f , height = thr , distance = int (0.25 * fs ))
6161 gpk = pk + a
6262 keep = gpk >= pos
6363 peaks .extend (gpk [keep ].tolist ())
@@ -68,16 +68,15 @@ def detect_all(path="ECG.DAT"):
6868 order = np .argsort (peaks )
6969 peaks , amps = peaks [order ], amps [order ]
7070
71- # T-wave rejection: a peak soon after a much larger one
72- keep = np .ones (peaks .size , bool )
71+ keep = np .ones (peaks .size , bool ) # T-wave rejection
7372 for i in range (1 , peaks .size ):
74- if (peaks [i ] - peaks [i - 1 ]) / FS < TWAVE_WIN and amps [i ] < TWAVE_RATIO * amps [i - 1 ]:
73+ if (peaks [i ] - peaks [i - 1 ]) / fs < TWAVE_WIN and amps [i ] < TWAVE_RATIO * amps [i - 1 ]:
7574 keep [i ] = False
7675 return peaks [keep ], amps [keep ], total
7776
7877
79- def classify (peaks ):
80- rr = np .diff (peaks ) / FS
78+ def classify (peaks , fs ):
79+ rr = np .diff (peaks ) / fs
8180 n = peaks .size
8281 med = np .empty (rr .size )
8382 for i in range (rr .size ):
@@ -92,10 +91,9 @@ def classify(peaks):
9291 return rr , premature , pause
9392
9493
95- def signal_quality_mask (peaks , rr ):
96- """Flag ~10 s windows as unreliable (noise): extreme R-R variability, an
97- implausible sustained rate, or too many physiologically-fast intervals."""
98- WIN = 10 * FS
94+ def signal_quality_mask (peaks , rr , fs ):
95+ """Flag ~10 s windows as unreliable (noise)."""
96+ WIN = 10 * fs
9997 good = np .ones (rr .size , bool )
10098 for e in np .arange (0 , peaks [- 1 ] + WIN , WIN ):
10199 m = (peaks [:- 1 ] >= e ) & (peaks [:- 1 ] < e + WIN )
@@ -109,41 +107,71 @@ def signal_quality_mask(peaks, rr):
109107 return good
110108
111109
112- def main ():
113- print ("== Ectopic-burden estimate (rhythm-based) ==" )
114- peaks , amps , nsamp = detect_all ()
115- dur_h = nsamp / FS / 3600
116- rr , premature , pause = classify (peaks )
117- good = signal_quality_mask (peaks , rr )
110+ def _resolve_start (args ):
111+ if args .start :
112+ return datetime .datetime .fromisoformat (args .start )
113+ try : # read from FC1_PRG.DAT if present
114+ import meditech_decode
115+ s = meditech_decode .parse_prg (args .prg ).get ("recording_start" )
116+ if s :
117+ return s
118+ except Exception :
119+ pass
120+ return datetime .datetime (1970 , 1 , 1 )
121+
122+
123+ def _parse_args (argv = None ):
124+ p = argparse .ArgumentParser (description = "Rough rhythm-based ectopic-burden estimate "
125+ "for a CardioMera Holter recording (NOT a diagnosis)." )
126+ p .add_argument ("ecg" , nargs = "?" , default = "ECG.DAT" )
127+ p .add_argument ("--prg" , default = "FC1_PRG.DAT" , help = "config file for the start timestamp" )
128+ p .add_argument ("--start" , default = None , help = "override start time (ISO 8601)" )
129+ p .add_argument ("--channels" , type = int , default = 3 )
130+ p .add_argument ("--rate" , type = int , default = 300 )
131+ p .add_argument ("--header-len" , type = int , default = 13 )
132+ p .add_argument ("--detect-channel" , type = int , default = 2 , help = "channel index for R detection" )
133+ p .add_argument ("--events-out" , default = "arrhythmia_events.npz" )
134+ return p .parse_args (argv )
135+
136+
137+ def main (argv = None ):
138+ a = _parse_args (argv )
139+ start = _resolve_start (a )
140+ print ("== Ectopic-burden estimate (rhythm-based; NOT a diagnosis) ==" )
141+ peaks , amps , nsamp = detect_all (a .ecg , a .rate , a .channels , a .header_len , a .detect_channel )
142+ dur_h = nsamp / a .rate / 3600
143+ rr , premature , pause = classify (peaks , a .rate )
144+ good = signal_quality_mask (peaks , rr , a .rate )
118145 total = peaks .size
119146 hr = 60 / rr
120-
121147 npre_all = int (premature .sum ())
122148 npre_clean = int ((premature [:- 1 ] & good ).sum ())
123149 nbeat_clean = int (good .sum ())
124150
125151 print (f"Length: { dur_h :.2f} h, total beats: { total :,} " )
126- print (f"Median HR: { np .median (hr ):.0f} bpm (p2 { np .percentile (hr ,2 ):.0f} , p98 { np .percentile (hr ,98 ):.0f} )" )
127- print (f"Premature beats (all): { npre_all :,} = { 100 * npre_all / total :.2f} % of beats" )
128- print (f"Premature beats (quality-gated): { npre_clean :,} = { 100 * npre_clean / max (nbeat_clean ,1 ):.2f} %" )
152+ print (f"Median HR: { np .median (hr ):.0f} bpm (p2 { np .percentile (hr ,2 ):.0f} , "
153+ f"p98 { np .percentile (hr ,98 ):.0f} )" )
154+ print (f"Premature beats (all): { npre_all :,} = { 100 * npre_all / max (total ,1 ):.2f} % of beats" )
155+ print (f"Premature beats (quality-gated): { npre_clean :,} = "
156+ f"{ 100 * npre_clean / max (nbeat_clean ,1 ):.2f} %" )
129157 print (f"Excluded as noise/unreliable: { 100 * (1 - good .mean ()):.1f} % of intervals" )
130158 print (f"Pauses > { PAUSE_S :.0f} s: { int (pause .sum ())} " )
131159
132- hours = (peaks / FS / 3600 ).astype (int )
160+ hours = (peaks / a . rate / 3600 ).astype (int )
133161 print ("\n Per-hour burden (recording hour : beats : ectopics : %):" )
134162 for hh in range (int (dur_h ) + 1 ):
135163 m = hours == hh
136164 tb = int (m .sum ())
137165 if tb == 0 :
138166 continue
139167 eb = int (premature [m ].sum ())
140- ts = (START + datetime .timedelta (hours = hh )).strftime ("%m-%d %H:%M" )
168+ ts = (start + datetime .timedelta (hours = hh )).strftime ("%m-%d %H:%M" )
141169 bar = "#" * min (40 , int (60 * eb / max (tb , 1 )))
142170 print (f" h{ hh :2d} [{ ts } ]: { tb :5d} : { eb :4d} : { 100 * eb / tb :5.2f} % { bar } " )
143171
144- np .savez ("arrhythmia_events.npz" , peaks = peaks , premature = premature ,
145- pause = pause , good = np .append (good , False ), fs = FS )
146- print ("\n Events saved: arrhythmia_events.npz " )
172+ np .savez (a . events_out , peaks = peaks , premature = premature , pause = pause ,
173+ good = np .append (good , False ), fs = a . rate )
174+ print (f "\n Events saved: { a . events_out } " )
147175
148176
149177if __name__ == "__main__" :
0 commit comments