Who are you, Contestant S1?

Before Any Analysis לפני שמתחילים לנתח

DETECTIVE — "Stage 0 only reads the witness. Stage 1 is where you make sure the room itself isn't lying to you before anyone else gets questioned." שלב 0 רק קורא את העד. שלב 1 הוא איפה שמוודאים שהחדר עצמו לא משקר, לפני שמישהו נוסף נחקר.

Stage 0 (see The 9 Stats) only characterizes the raw signal — it changes nothing. Stage 1 is the first stage that actually touches the data: three housekeeping steps, none of which decide what counts as an artifact and none of which touch the raw HDF5 files — each works on a copy. שלב 0 (ראו 9 המדדים) רק מאפיין את האות הגולמי — לא משנה כלום. שלב 1 הוא השלב הראשון שבאמת נוגע בנתונים: שלושה צעדי הכנה, אף אחד מהם לא מחליט מה נחשב ארטיפקט ואף אחד לא נוגע בקובצי ה-HDF5 הגולמיים — כל אחד עובד על עותק.

1 · Normalization — tested, and it didn’t matter

Several standard normalization schemes were compared side by side early on. None of them changed the downstream picture in a meaningful way — so normalization choice is not a sensitive design decision in this pipeline, and the simplest option was kept. כמה שיטות נרמול נבדקו זו מול זו בשלב מוקדם. אף אחת לא שינתה את התמונה בהמשך באופן משמעותי — כך שבחירת הנרמול אינה החלטת עיצוב רגישה כאן, ונשמרה האפשרות הפשוטה ביותר.

2 · Centering — remove the global median

Every channel is shifted so its own global median sits at zero — purely for PSD stability, not a scientific claim about baseline. The median is used (not the mean) because it’s robust to the saturation spikes already known to exist in some channels. כל ערוץ מוזז כך שהחציון הגלובלי שלו יושב באפס — למען יציבות ה-PSD בלבד, לא כטענה מדעית על הבסיס. נבחר חציון (לא ממוצע) כי הוא עמיד מול קפיצות רוויה שכבר ידוע שקיימות בחלק מהערוצים.

3 · The 50/100 Hz problem — Israel’s power grid

Israel’s electrical grid runs at 50 Hz, with its first harmonic at 100 Hz. Both leak into the recordings as sharp spectral lines that have nothing to do with the brain. רשת החשמל בישראל היא 50Hz, וההרמוניה הראשונה יושבת ב-100Hz. שתיהן מדליפות להקלטה כקווים ספקטרליים חדים שאין להם קשר למוח.

The fix: epoch-wise adaptive spectral capping

Instead of a fixed filter, the correction runs per 6-second epoch (\(N=6000\) samples) and only caps a peak that’s actually too loud — it never forces the frequency to zero:

1. FFT the epoch.
2. Around each target \(f_0\in\{50,100\}\,\text{Hz}\), search \(f_0\pm1\,\text{Hz}\) for the real local peak — not just the exact bin.
3. Estimate the local spectral floor from the median magnitude in \(f_0\pm5\,\text{Hz}\), excluding the peak’s own \(\pm1\,\text{Hz}\) zone.
4. If the peak exceeds a factor of the floor (here, \(1.4\times\)), scale down only the bins within \(\pm0.25\,\text{Hz}\) of the peak to that cap.
5. Preserve phase and conjugate symmetry, then inverse-FFT to reconstruct the epoch.

Not every channel needs it

The correction only runs where a locked, manually-reviewed policy says a specific channel of a specific recording actually shows the line-noise problem — e.g. E12_6m channels 1&2, or S9_9m channels 1&2. Most recordings need no correction at all. התיקון מופעל רק היכן שמדיניות נעולה, שנבדקה ידנית, קובעת שערוץ ספציפי בהקלטה ספציפית באמת מציג את בעיית רעש הרשת. רוב ההקלטות לא זקוקות לתיקון בכלל.

Recording length & the exceptions

Each recording runs about 17 minutes at \(f_s=1000\,\text{Hz}\) — nearly all files hold exactly 1,048,575 samples. Three recordings are shorter, for reasons unrelated to the correction above:

Vital Signs: the 9 statistics 9 המדדים

DETECTIVE — "Every channel gets the same nine questions. No vibes. A number, a threshold, a verdict." כל ערוץ מקבל את אותן תשע שאלות. בלי תחושות — מספר, סף, פסק-דין.

S1 — “…” (S1 doesn’t talk. This is the closest thing to a statement it makes — nine numbers, whether it likes it or not.) "..." (S1 לא מדבר. זה הכי קרוב שיש למשפט ממנו — תשעה מספרים, בין אם הוא רוצה או לא.)

For every channel of every recording (\(N\) samples at \(f_s=1000\,\text{Hz}\)) we compute nine metrics. Each targets one failure mode of an LFP signal. Let \(x[n]\) be the signal and \(\Delta x[n]=x[n]-x[n-1]\) its sample-to-sample difference.

1 · RMS — overall energy

The signal’s effective amplitude. Abnormally high RMS flags a noisy or artifact-dominated channel. אמפליטודה אפקטיבית — ערך גבוה מסמן ערוץ רועש.

2 · JumpRate — motion jerks

A robust noise scale is estimated from the differences using the Median Absolute Deviation:

Percentage of samples that jump more than 8 robust-sigmas — the fingerprint of movement. The \(1.4826\) makes MAD a consistent estimator of \(\sigma\) for Gaussian noise. אחוז הדגימות שקופצות מעל 8 סיגמות רובוסטיות — טביעת האצבע של תנועה.

3 · Kurtosis — spikiness

The 4th standardized moment. A Gaussian sits at \(3\); higher values mean heavy tails — rare, sharp spikes riding on the signal. מומנט רביעי — ערך גבוה = פסגות חדות וחריגות.

4 · MaxAbsJump — the worst single jerk

Magnitude of the single largest step — the size of the worst movement event.

5 · Drift (EdgeMedianDiff) — baseline wander

Slow baseline shift across the recording. Medians (not means) make it robust to spikes. סחף בסיס איטי לאורך ההקלטה — חציונים עמידים לחריגים.

6 · Delta-Ratio — low-frequency dominance

Fraction of total power sitting in the slow delta band. Motion and drift dump energy low; a high ratio is a technical red flag (not a biological claim at this stage). שבר האנרגיה בפס האיטי — תנועה/סחף שופכים אנרגיה נמוך.

7–8 · Saturation counts — hitting the wall

How many samples are stuck at the exact minimum / maximum — the signature of ADC clipping at the \(\pm427\,\mu V\) rail. (We subtract 1 because every signal has one natural min and max.) כמה דגימות תקועות בדיוק על המינימום/מקסימום — חתימת רוויית ADC.

9 · Robust range — dynamic span

The 1st-to-99th percentile span — the amplitude range while ignoring the most extreme 1% on each side (outlier-proof).

Each channel gets labelled Failure / High / Moderate / Clean by running its nine numbers against a locked set of thresholds — no single metric decides alone.

DSP & Wavelets — the Analysis Toolkit כלי הניתוח

LUMINA — "Time tells you when. Frequency tells you what. Put them together and the signal has nowhere left to hide." הזמן מספר מתי. התדר מספר מה. יחד — לאות אין איפה להסתתר.

Foundations: sampling & the spectrum

Sampling at \(f_s=1000\,\text{Hz}\) means the highest representable frequency is the Nyquist limit:

We read the signal’s energy across biological bands via its power spectral density \(P(f)\):

Time–frequency: the Wavelet scalogram

A raw trace shows voltage over time; an FFT shows frequency but loses when. Motion artifacts are not always sharp spikes — they can be slow baseline shifts, decays, recoveries. To catch them we need both axes at once. The Continuous Wavelet Transform slides a scaled/shifted wavelet \(\psi\) across the signal:

where \(a\) is scale (∝ 1/frequency) and \(b\) is time. The result — a scalogram — is an energy map: time on X, frequency on Y, colour = energy. This is exactly what Stage 2 actually runs on: a Morlet scalogram across 2–250Hz, Z-scored per frequency row. A brain rhythm lives in one narrow frequency lane; a real artifact lights up many adjacent lanes at once.

TQWT — building the alibi

The Tunable-Q Wavelet Transform (Selesnick) sorts signal by resonance, not just frequency: sustained oscillations (brain rhythms) are high-Q, abrupt transients are low-Q. This project uses that selectivity for construction, not classification — per channel, per 6-second window, it rebuilds the signal from only the sub-bands that fall in one range:

X_ref never decides what's an artifact — Stage 2's rule-based detector already did that; X_ref is just the "clean" stand-in, ready to go there.

Longitudinal fingerprint — is 6-month S1 still S1?

A comparison between 5xFAD and WT only means something if the same mouse, recorded again at a later age, is still comparable to its earlier self — not a different recording setup wearing the same mouse ID. Each mouse gets a spectral fingerprint (its own PSD shape) computed independently at 4, 6, and 9 months, and the fingerprints are correlated pairwise across ages.

Choosing who’s allowed to testify — pre-registered QC

29 of 32 recordings make it into the analysis. The 3 that don’t (see Prep for the exact list) are excluded because they failed a quality threshold derived from Stage 0’s statistics — decided before any group comparison was run, not after looking at results. That ordering is the entire point: a pass/fail rule fixed in advance and documented is reproducible quality control. The same exclusion made after seeing which recordings are inconvenient for the hypothesis is no longer quality control — it’s selecting the answer.

The Correction Pipeline צינור התיקון

DETECTIVE — "Four stages. Not a pile-up — a chain of custody. Detect, build an alibi, substitute it in, then check the receipts." ארבע שלבים. לא ערימה — שרשרת. לזהות, לבנות אליבי, להציב אותו, ואז לבדוק את הקבלות.

Stage 2 — the wavelet detector

A continuous wavelet transform (Morlet, 2–250Hz) turns each 60-second block into a scalogram, Z-scored per frequency row against its own median/MAD. Two severity tiers, plus an amplitude check merged into the milder one:

Whatever RED already claims is carved out of BLUE first — the two never overlap in the final mask. An amplitude jump only counts toward BLUE if at least 2 of the 4 channels register it in the same window — one twitchy wire alone doesn’t qualify. כל מה ש-RED תופס נגזר מ-BLUE מראש — השניים לעולם לא חופפים במסכה הסופית. קפיצת עוצמה נספרת רק אם לפחות 2 מתוך 4 הערוצים מזהים אותה באותו חלון.

Stage 3 — building the alibi

A fully independent, per-channel reference signal — no cross-channel averaging at all. Every 6-second window, TQWT rebuilds the signal from only the sub-bands whose center frequency lands in one range:

In practice the toolbox’s actual sub-band centers resolve to 4.105–95.627Hz, not a clean 4–105. This reference doesn’t decide what’s an artifact — Stage 2 already did that. X_ref is the "clean" stand-in, ready to be substituted wherever a method needs it.

Stage 3.5 — the second alibi

The pure TQWT reference has a gap: it starts at ~4Hz, so it says nothing about the slow delta band. Stage 3.5 fills that gap with a Butterworth bandpass and adds it in:

The result covers 0.5–95.627Hz — the full biologically relevant range, not just 4Hz and up. From here on, every correction stage runs twice: once fed the pure TQWT reference, once fed this wider Unified one — two alibis, compared side by side. מכאן ואילך כל שלב תיקון רץ פעמיים: פעם עם האזין הטהור של TQWT, פעם עם ה-Unified הרחב יותר — שני אליבים, בהשוואה זה לצד זה.

Dead end — the suspect that walked

Before settling on CAR/Kalman/Khorasani, the most promising lead was SSA (Singular Spectrum Analysis): decompose the pre-artifact context, forecast forward into the damaged window, decompose the post-artifact context, forecast backward, then blend the two gradually across it — nicknamed Kalman Blend even though it isn’t a Kalman filter at all. The flaw: a forecast is a plausible statistical continuation, not a recovered measurement. LFP is non-stationary, so “the pre-artifact dynamics keep going through the gap” is a fragile assumption — and inside saturation windows the original information doesn’t exist to forecast toward in the first place. Forecasts were sometimes noisier than the raw signal, especially in short windows, and had to fall back to source.

One question, three answers

Stages 4, 5 and 6 all answer the same question with a different amount of freedom. Model each channel as x(t) = s(t) + w(t)·n(t) — measured = brain signal + weight × shared noise. The whole design space is what you let \(w\) be:

Stage 4 — Static Weighted CAR

The classic idea behind CAR (Common Average Reference): assume each channel’s genuine noise is independent, so subtracting the average across enough channels cancels the noise while leaving the shared signal alone:

With N=16 independent channels (the textbook case) that suppresses buried noise by \(10\log_{10}(1/N)\approx-12\)dB. This project has only 3 independent channels (CH1, CH2, CH4) — the same trick buys just -4.8dB, a real ~7dB gap. Flat averaging isn’t enough, so every channel gets its own weight instead of a blind \(1/N\) — computed once per recording, then applied to every sample:

Weight per channel: \(\displaystyle w_i = \frac{\text{rms}(X_{\text{pre}}-X_{\text{ref}})_i}{\sum_i \text{rms}(X_{\text{pre}}-X_{\text{ref}})_i}\) — one number per channel, for the whole recording. משקל לכל ערוץ לפי rms ההפרש בינו לבין הייחוס — מספר אחד לכל ערוץ, לכל ההקלטה.

Adaptive Kalman קלמן אדפטיבי

LUMINA — "One weight for a whole recording is a blunt instrument. Let it breathe, sample by sample." משקל אחד לכל הקלטה הוא כלי גס. תני לו לנשום, דגימה-דגימה.

Stage 4’s CAR weight is one static number per channel. Stage 5 lets that same weight drift over time — a real Kalman filter, not a metaphor. Model per channel \(i\): the residual \(d_i[k]=X_{\text{pre},i}[k]-X_{\text{ref},i}[k]\) is explained by a time-varying weight on the common mode \(c[k]=\text{mean}_j(X_{\text{pre}}[k,:])\):

Forward pass, then an RTS backward smoother refines the whole weight trajectory using future samples too — not just causal, the full recording is known in advance:

Correction: \(X_{\text{corrected},i}[k] = X_{\text{pre},i}[k] - w_i^{\text{RTS}}[k]\cdot c[k]\). התיקון: מחסרים מכל דגימה את המשקל המוחלק (RTS) כפול המצב המשותף.

Khorasani’s Filter המסנן של קורשאני

DETECTIVE — "We ruled this suspect out on paper. Turns out paper doesn’t run experiments." פסלנו את החשוד הזה על הנייר. מסתבר שהנייר לא מריץ ניסויים.

Khorasani (2019) estimates a dynamic AR(5) weight for the common mode instead of one static number — five taps of memory, not one:

where \(h(t)=[n(t),\,n(t{-}1),\,\dots,\,n(t{-}4)]\) is the common mode’s own recent history, \(\alpha=0.99\), and the Kalman update is standard (predict, gain, correct). Run in 3 input modes — fed raw X_pre, the TQWT reference, or the Unified one:

The Verdict: Which Method Wins פסק הדין — איזו שיטה מנצחת

S1 — “...” (13 method variants, 1,508 channel-recordings, one honest table.)

This is a methods verdict, not the biological one. It answers which correction technique to trust — not whether 5xFAD differs from WT in theta/gamma. That question is still open.

Genotype check (descriptive only, small cohort, no causal claim): 5xFAD and WT show near-identical correction behavior (+1.14dB vs +1.25dB) — the pipeline treats sick and healthy brains the same way. Whether their underlying theta/gamma actually differ is a separate analysis, not yet run.

Epilogue — the cost of finding out

Nine self-contained stages, each reading the previous one’s output from disk and writing its own — no shared memory between stages, so any one of them can be re-run and re-verified in isolation. 29 recordings clear the entire chain in about 25 minutes: 52.6 seconds per recording, all on a MacBook M3 with 8GB — no cluster, no GPU, no cloud.

What got built here isn’t an answer — it’s a laboratory: a documented, re-runnable, arguable chain. The question that started the investigation — does an Alzheimer’s-model brain really look different from a healthy one, once motion is properly removed — is still open. What exists now is the instrument that can finally ask it honestly. מה שנבנה כאן הוא לא תשובה — זו מעבדה: שרשרת מתועדת, ניתנת להרצה חוזרת ולערעור. השאלה שבגללה התחילה החקירה עדיין פתוחה. מה שיש עכשיו הוא המכשיר שסוף-סוף יכול לשאול אותה ביושר.