# EEG ERSP and ITC measures

**Created:** 11/2/2026 4:23:30  
**Updated:** 11/2/2026 5:39:13  
**Exported:** 11/2/2026 17:20:14  
**Link:** [https://claude.ai/chat/75a74557-2e74-4d25-8583-7a2b763b50fc](https://claude.ai/chat/75a74557-2e74-4d25-8583-7a2b763b50fc)  

## User:

> 11/2/2026 4:23:30

Hi claude, in the eeglab conversation i asked you to remove ersps from recordings to look at them without events. Is this something reasonable or complex/nonsensical to do?

## Assistant:

> 11/2/2026 4:23:43

> Deciding to answer the ERSP question directly.
>
> - Deciding to answer the ERSP question directly.
> - **Done**

It's reasonable, with one catch about what "remove" can mean.

**What an ERSP is:** the *average* change in spectral power around an event, relative to a baseline. It doesn't exist in the raw recording as a separate signal. It's a statistic computed over many events.

**Two workable versions of the idea:**

1. **Just look at the continuous spectrogram.** Compute time-frequency power on the whole recording, ignoring events. You see everything: the ongoing rhythms plus the event responses mixed in. This is simple and always valid.

2. **Subtract the event template, then look at the residual.** Compute the continuous spectrogram in dB. At each event onset, subtract the ERSP (the average event-locked power change). What's left is the ongoing activity plus each trial's deviation from the average. This is a real technique. It's the power-domain cousin of subtracting the average ERP from each trial to isolate "induced" activity.

**Where it gets shaky:**

- Power doesn't add linearly the way voltage does, so the subtraction is only approximate. Doing it in log power (dB) is the usual compromise.
- If events are closer together than the ERSP window, the templates overlap and double-subtract. You'd need a regression approach (fit all overlapping templates at once) instead of plain subtraction.
- The residual isn't "event-free brain." It still holds trial-to-trial variability in the response. That can be the interesting part.

So it's not nonsensical. Version 1 is trivial; version 2 is a legitimate, slightly fiddly analysis. Want runnable code for version 2 on one of the datasets you were exploring?

## User:

> 11/2/2026 4:25:04

Ok so the spectral power residual can be precisely calculated via log subtraction (division?) — but going backward to the time domain would require an advancement in ersps?

## Assistant:

> 11/2/2026 4:25:13

> Weighing whether phase loss allows exact time-domain inversion.
>
> - Weighing whether phase loss allows exact time-domain inversion.
> - **Done**

Yes on the first part: subtracting in dB *is* division in linear power. An ERSP is already a ratio to baseline, so dividing it out is exact arithmetic. The only approximation is the modeling assumption that the event effect is multiplicative.

On going back to the time domain: no advancement needed. You just have to keep the phase instead of throwing it away.

- **Power alone can't be inverted.** It discards phase, so infinitely many signals share the same spectrogram.
- **The complex transform can be inverted.** Keep the STFT coefficients (magnitude and phase). Scale each coefficient's magnitude by `1 / sqrt(ERSP gain)` at the matching time and frequency, leave the phase alone, then run the inverse STFT.

This is the same trick audio denoising uses (spectral subtraction and Wiener-style gain masks). It's standard, not new.

Two caveats:

- **After editing, the coefficient grid may not match any real signal exactly.** The inverse STFT returns the closest real signal in the least-squares sense. In practice that works fine when the gain changes smoothly.
- **Use STFT, not EEGLAB's default wavelets.** `newtimef` wavelets are redundant and don't come with a clean inverse. An STFT with a proper overlap window (e.g. Hann at 50%) inverts cleanly with `scipy.signal.istft`.

So the pipeline is: STFT, divide out the event gain around each event, ISTFT. The result is a time-domain recording with the average spectral response removed.

## User:

> 11/2/2026 4:28:04

How do you correctly take the mean of a signal with phase information?

## Assistant:

> 11/2/2026 4:28:17

> Working out the circular mean of phase angles.
>
> - Working out the circular mean of phase angles.
> - **Done**

It depends on what you want the mean to keep. There are three standard ways, and each answers a different question:

1. **Mean of the complex coefficients** (coherent average). Average the complex values directly. Components with consistent phase across trials survive; random-phase components cancel toward zero. This is the phase-locked (evoked) part, the same thing you'd get from the spectrum of the ERP.

2. **Mean of the power** (`|z|²`, phase discarded). Nothing cancels. This is total power, and its baseline-normalized form is what the ERSP shows.
   - Total minus the coherent part's power gives **induced power**: activity that is event-related but not phase-locked.

3. **Mean of unit phasors** (`z / |z|`). Each trial gets an equal vote regardless of amplitude.
   - The magnitude of the result is **inter-trial coherence** (ITC): 0 means random phase, 1 means perfectly locked.
   - Its angle is the **circular mean phase**.

The common mistake is averaging phase *angles* arithmetically. That breaks at the wraparound: the mean of 359° and 1° should be 0°, not 180°. Always average on the unit circle instead.

```python
import numpy as np

# tf_coefficients: complex STFT values, shape (trials, freqs, times)
rng = np.random.default_rng(0)
trials, freqs, times = 50, 4, 10
locked_phase = 0.5                                   # radians, shared by every trial
random_phase = rng.uniform(-np.pi, np.pi, (trials, freqs, times))
tf_coefficients = (np.exp(1j * locked_phase)         # phase-locked component
                   + np.exp(1j * random_phase))      # random-phase component, same amplitude

# 1. Coherent average: random phases cancel, the locked part survives
evoked_power = np.abs(tf_coefficients.mean(axis=0)) ** 2

# 2. Total power: average |z|^2, so nothing cancels
total_power = (np.abs(tf_coefficients) ** 2).mean(axis=0)
induced_power = total_power - evoked_power           # event-related but not phase-locked

# 3. Unit phasors: amplitude ignored, phase consistency only
unit_phasors = tf_coefficients / np.abs(tf_coefficients)
mean_phasor = unit_phasors.mean(axis=0)
inter_trial_coherence = np.abs(mean_phasor)          # 0 = random phase, 1 = locked
circular_mean_phase = np.angle(mean_phasor)          # correct across the ±π wraparound

print("evoked  ~1:", evoked_power.mean().round(2))
print("total   ~2:", total_power.mean().round(2))
print("induced ~1:", induced_power.mean().round(2))
print("ITC        :", inter_trial_coherence.mean().round(2))
print("phase ~0.5 :", circular_mean_phase.mean().round(2))
```

For the residual pipeline we discussed, the choice matters:
- Dividing out the **power mean (2)** removes the average power change but leaves each trial's phase.
- Subtracting the **complex mean (1)** removes only the phase-locked part, which is the time-domain ERP subtraction.

Doing both, in that order, gives you something close to "ongoing activity."

## User:

> 11/2/2026 4:30:16

It sounds like the only time we have information that is novel after the erp is when we discard phase? Or do you propose unit phasors have it?

## Assistant:

> 11/2/2026 4:30:26

> Comparing induced power and inter-trial coherence against the ERP.
>
> - Comparing induced power and inter-trial coherence against the ERP.
> - **Done**

Unit phasors do carry novel information, so it's not only when you discard phase.

The ERP is the complex mean: amplitude and phase consistency multiplied together. A small ERP could mean either of two things:
- weak activity that is perfectly phase-locked, or
- strong activity with sloppy phase.

ITC separates those two. It measures consistency alone, so it tells you something the ERP can't.

The general pattern: the ERP keeps only the **first moment** of each trial's complex value (its plain average). Anything about how trials differ from each other is lost. New information comes from other statistics of the same data:

| Statistic | What it adds beyond the ERP |
|---|---|
| Mean power `\|z\|²` | Energy that cancels in the average (induced activity) |
| Mean unit phasor (ITC) | Phase consistency, separated from amplitude |
| Full phase distribution | Shape, e.g. two preferred phases, which a mean averages into mush |
| Phase *between* channels or frequencies | Connectivity (phase-locking value, coherence) and cross-frequency coupling |

That last row is the key correction to "discard phase." What gets averaged away is phase **relative to the event**. Phase **relative to another channel** survives fine, and it's a large part of what's left to study after the ERP.

## User:

> 11/2/2026 4:32:31

Help me understand: itc has information lost in erp, and relates to phase difference

## Assistant:

> 11/2/2026 4:32:42

> Showing that ITC and PLV share one computation.
>
> - Showing that ITC and PLV share one computation.
> - **Done**

Every phase you measure is a difference against some reference. With no reference specified, it's the event clock: "where was this oscillation in its cycle at time t after the stimulus?"

**ITC is consistency of that phase difference across trials.** Picture each trial as an arrow:

- **Arrow angle** = the oscillation's phase relative to event onset
- **ERP** = average the arrows as they are, so long arrows count more
- **ITC** = shrink every arrow to length 1 first, then average; the result's length is how tightly the angles cluster

**What the ERP loses:** for amplitude and phase that vary independently across trials, roughly

`|ERP| ≈ mean amplitude × ITC`

The ERP only reports the product. Two very different brains give the same ERP:

```python
import numpy as np
rng = np.random.default_rng(0)
trials = 10000

def summarize(name, amplitude, phase):
    coefficients = amplitude * np.exp(1j * phase)
    erp_magnitude = np.abs(coefficients.mean())
    itc = np.abs(np.exp(1j * phase).mean())          # unit arrows, amplitude ignored
    print(f"{name:28s} |ERP|={erp_magnitude:.2f}  ITC={itc:.2f}  amplitude={amplitude}")

# Weak but perfectly locked: every arrow points the same way
summarize("weak, locked", 1.0, np.zeros(trials))

# Strong but sloppy: phase spread with a von Mises distribution
# (a circular normal; kappa sets concentration, chosen so ITC is ~0.5)
summarize("strong, sloppy", 2.0, rng.vonmises(0.0, 1.16, trials))
```

Both print |ERP| ≈ 1.0, but ITC is 1.0 versus about 0.5. The ERP can't tell "small and reliable" from "big and jittery." ITC can.

**Link to channel-to-channel phase difference:** swap the reference. Instead of each trial's phase relative to the event, use its phase relative to another channel in the same trial. The identical formula is then the **phase-locking value** (PLV): it asks whether two regions hold a consistent phase relationship, whether or not either is locked to the stimulus.

So it's one computation, "how consistent is this phase difference?", with the event or another channel as the reference.

## User:

> 11/2/2026 4:34:28

Can you make me a small in-browser comparison of the information from the different measures, including ersp?

## Assistant:

> 11/2/2026 4:38:55

> - Choosing between an artifact and inline output.
> - Loaded a skill
> - Styling the plot panels and sizing the canvases.
> - Created a file
> - Ran a command
> - Used artifacts
> - **Done**

It's a bench of 60 simulated trials at 10 Hz, with an evoked burst near 0.15 s and an induced burst near 0.5 s. You get six panels: ERP, ERSP, ITC, the evoked/induced power split, live phasor dials, and a readout at the cursor. Each panel carries its equation as code.

Two presets are worth trying first:
- **Weak + locked** and **Strong + sloppy** give nearly the same ERP. The ERSP and ITC tell them apart.
- **Induced only** shows a flat ERP and near-zero ITC, but a clear ERSP blob.

## User:

> 11/2/2026 5:02:15

What is W?

## Assistant:

> 11/2/2026 5:02:24

> Explaining what the complex wavelet coefficient means.
>
> - Explaining what the complex wavelet coefficient means.
> - **Done**

W is the complex wavelet coefficient: one complex number per trial, per frequency, per time point.

You get it by sliding a short "probe" wave along the signal. The probe is a 10 Hz sine inside a Gaussian window (a Morlet wavelet). At each time point you multiply the probe against the signal and sum, which asks "how much 10 Hz is here, right now, and where in its cycle is it?"

- `abs(W)`: the local amplitude of that frequency
- `angle(W)`: the local phase, relative to the event clock
- `abs(W)**2`: the local power

So `W` is the arrow drawn for each trial in the left phasor dial. Every measure on the page is a different way of averaging those arrows across trials.

```python
import numpy as np

def wavelet_coefficient(signal, center_sample, hz, sampling_rate, cycles=4):
    # Gaussian width in seconds; more cycles = sharper frequency, blurrier time
    sigma_seconds = cycles / (2 * np.pi * hz)
    half_length = int(3 * sigma_seconds * sampling_rate)
    lag_seconds = np.arange(-half_length, half_length + 1) / sampling_rate

    window = np.exp(-0.5 * (lag_seconds / sigma_seconds) ** 2)
    probe = window * np.exp(-2j * np.pi * hz * lag_seconds)   # conjugated so angle(W) is the signal's phase

    segment = signal[center_sample - half_length : center_sample + half_length + 1]
    return 2 * np.sum(segment * probe) / np.sum(window)       # scaled so a cosine of amplitude A gives abs(W) == A

# Check: a 10 Hz cosine with amplitude 1.5 and phase 0.7 rad
sampling_rate = 250
times = np.arange(0, 2, 1 / sampling_rate)
signal = 1.5 * np.cos(2 * np.pi * 10 * times + 0.7)
center_sample = 250                                           # t = 1.0 s, so expected phase is 2π·10·1.0 + 0.7 ≡ 0.7
coefficient = wavelet_coefficient(signal, center_sample, 10, sampling_rate)
print(abs(coefficient), np.angle(coefficient))                # ≈ 1.5, ≈ 0.7
```

It's the same object as an STFT coefficient, just with a Gaussian window whose width scales with frequency.

## User:

> 11/2/2026 5:09:18

Add the spectral analysis of the erp to the output to compare information content.

## Assistant:

> 11/2/2026 5:10:29

> Checking the lazy getter's initialization order.
>
> - Adding a wavelet heatmap for the ERP.
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Checking the lazy getter's initialization order.
> - Ran a command
> - **Done**

The difference is float rounding only, so the identity holds empirically.

> - Used artifacts
> - **Done**

Added a "Spectrum of the ERP" heatmap beside the ERSP, on the same baseline and dB scale. The ERP's spectrum also appears as a dotted line on the power-split plot and as a readout row.

The result: the ERP spectrum is exactly the evoked power. The wavelet is linear, so averaging and then transforming gives the same thing as transforming and then averaging. I checked this numerically, and the largest difference was 2e-7, which is just float rounding.

So the ERP's spectrum adds nothing new. Everything the ERSP shows beyond it is induced power. **Induced only** makes that starkest: the ERP map goes blank while the ERSP lights up.

I made this a general `ErpSpectrum` class and gave `drawHeatmap` a `source` argument, so it draws any measure without a special case.

## User:

> 11/2/2026 5:12:35

Ok so the itc and erp are both evoked — their spectra look similar but are plotted differently, and you say that itc has information not derivable from erp?

## Assistant:

> 11/2/2026 5:12:50

> Checking whether one measure can be derived from another.
>
> - Checking whether one measure can be derived from another.
> - **Done**

Both measure phase-locking, but they weight trials differently, so neither can be derived from the other alone.

- **ERP spectrum** `abs(mean(W))**2` gives big trials more say.
- **ITC** `abs(mean(W/abs(W)))` gives every trial one equal vote.

They look alike in the bench because the simulation varies amplitude and phase independently. In that case `|ERP| ≈ mean(|W|) · ITC`, so the ERP **plus the ERSP** recovers ITC approximately, but the ERP alone never can. **Weak + locked** versus **Strong + sloppy** shows this: matching ERP spectra, different ITC.

The approximation breaks when amplitude and phase are correlated. For example, half the trials are large and locked, and the other half are small and random:

```python
import numpy as np
rng = np.random.default_rng(0)
half = 5000

big_locked = 2.0 * np.exp(1j * np.zeros(half))                        # strong trials, same phase
small_random = 0.2 * np.exp(1j * rng.uniform(-np.pi, np.pi, half))    # weak trials, random phase
coefficients = np.concatenate([big_locked, small_random])

erp_magnitude = abs(coefficients.mean())
mean_amplitude = abs(coefficients).mean()
itc = abs((coefficients / abs(coefficients)).mean())

print("ERP / mean amplitude:", round(erp_magnitude / mean_amplitude, 2))   # ~0.91, looks highly locked
print("ITC                 :", round(itc, 2))                               # ~0.50, only half the trials lock
```

The amplitude-weighted view says "almost perfectly locked," because the big trials dominate it. ITC says "half the trials." Only ITC sees that per-trial structure.

The plots also differ in scale. The ERP spectrum is in dB relative to baseline power; ITC is a linear 0–1 consistency with a noise floor near 0.11 at 60 trials.

Want me to add this "big-locked, small-random" case as a preset, plus a panel showing `abs(mean(W)) / mean(abs(W))` beside ITC so the divergence is visible?

## User:

> 11/2/2026 5:14:29

Yeah. Im confused on how its helpful to normalize a raw measure when they’re clearly summed from many mixed signals

## Assistant:

> 11/2/2026 5:15:50

> Weighing when per-trial normalization can actually help.
>
> - Weighing when per-trial normalization can actually help.
> - **Done**

Your skepticism is well-founded. Then I'll add the preset and panel.

**What normalizing can't do:** it doesn't unmix anything. W's phase is the phase of the whole mixture at that frequency, so it follows whichever source is biggest in that trial. A perfectly locked but weak source buried under strong random alpha gives low ITC. ITC of a mixture isn't the ITC of any one source.

**Why people normalize anyway:**
- **Robustness.** One huge artifact trial can't dominate the result.
- **A known null.** Under random phase, ITC's distribution depends only on trial count, not on amplitude or gain. That gives an easy significance test (the Rayleigh test).
- **Comparability.** It doesn't change with electrode impedance, reference, or subject gain.

**The real point of your objection:** for "locked signal plus independent noise," ITC is just a monotonic function of signal-to-noise ratio. That makes it redundant with the ERP plus the ERSP. ITC only adds information when amplitude and phase are correlated across trials. That's also why EEGLAB runs ICA first: unmix into components, then compute ITC per component, where the normalization means something.

> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Edited a file
> - Ran a command
> - Used artifacts
> - **Done**

I added a full-width "ITC vs amplitude-weighted locking" panel and a **Big locked, small random** preset. The preset is driven by two new sliders: locked trial fraction and unlocked trial size. I measured both locking values at the evoked peak for three presets:

| Preset | Weighted locking | ITC |
|---|---|---|
| Both | 0.90 | 0.89 |
| Strong + sloppy | 0.57 | 0.58 |
| Big locked, small random | 0.86 | **0.50** |

The two measures agree whenever amplitude and phase are independent. That's your point: in that case ITC is redundant with the ERP plus the ERSP. They split only when the big trials are the locked ones. Then ITC reports that half the trials lock, which nothing amplitude-weighted can see.

## User:

> 11/2/2026 5:34:41

What popularized ITC?

I’m thinking of alternatives to subtracting an evoked signal in time space, and wondering if it would work to construct a periodic waveform for each frequency bin and signal (or across signals with phase change but that seems harder), giving a waveform or phase for each signal associated with the ersp.  Is the line of thinking a wild goose chase?  If not, is there an existing similar approach?

## Assistant:

> 11/2/2026 5:35:08

> Weighing whether the literature search is needed.
>
> - Weighing whether the literature search is needed.
> - **Done**

**Where ITC came from:** Tallon-Baudry and colleagues (1996, *J Neurosci*) used the same measure under the name "phase-locking factor" to separate phase-locked from non-phase-locked 40 Hz gamma. Makeig's group named it ITC and made it mainstream. Two things did that: EEGLAB (Delorme & Makeig 2004), which plots ITC beside every ERSP by default, and Makeig et al. 2002 in *Science*, which argued that some ERP features come from phase resetting of ongoing rhythms. That claim set off a long debate in which ITC was the central evidence.

**Your idea isn't a wild goose chase**, but the first half already exists, and the second half has one trap.

1. **"A waveform per bin, per trial" is W itself.** `abs(W)` and `angle(W)` over time is exactly a per-trial, per-band oscillation; `real(W)` is the band-passed trace. Because the transform is linear, subtracting the evoked part in this domain is clean:

   ```python
   import numpy as np
   rng = np.random.default_rng(0)
   # tf_coefficients: complex W, shape (trials, times) for one frequency
   tf_coefficients = (np.exp(1j * 0.5) + rng.normal(size=(60, 100)) + 1j * rng.normal(size=(60, 100)))

   evoked = tf_coefficients.mean(axis=0)                # == W(erp), by linearity
   induced_coefficients = tf_coefficients - evoked      # complex subtraction, per trial, phase intact

   induced_power = (abs(induced_coefficients) ** 2).mean(axis=0)
   total_minus_evoked = (abs(tf_coefficients) ** 2).mean(axis=0) - abs(evoked) ** 2
   print(np.allclose(induced_power, total_minus_evoked * 60 / 60))   # True: same split, but you keep per-trial waveforms
   ```

   So you get a phase-bearing induced waveform for each trial, not just a power number.

2. **The trap: fitting each trial's own evoked waveform.** If every trial gets a freely fitted amplitude and phase per bin and you subtract it, you remove *everything* at that bin, induced activity included. The fit can't tell evoked from induced. To separate them you need a constraint: a shared template shape, with only latency, phase, or gain free per trial.

**Existing methods that do the constrained version:**
- **Woody filter** (1967): aligns each trial's latency to a template, then averages.
- **RIDE** (residue iteration decomposition): separates components with different per-trial latencies, then subtracts each trial's own jittered evoked estimate.
- **Phase-sorted ERP images** (EEGLAB `erpimage`): sort trials by phase at a bin to see what phase alignment does.
- **Across channels:** complex PCA, or generalized eigendecomposition (Mike X Cohen's GED tutorials). These find a spatial filter that maximizes phase-locked versus total power, which is roughly your "across signals with phase change" idea.

Your line of thinking lands on a real open problem: removing evoked activity that jitters from trial to trial without eating the induced activity. RIDE is the closest existing answer.

## User:

> 11/2/2026 5:38:46

Sorry, i meant induced. To subtract the ersp. All these other measures and approaches sound important though!

## Assistant:

> 11/2/2026 5:39:13

> Comparing the gain approach with the subtraction model.
>
> - Comparing the gain approach with the subtraction model.
> - **Done**

That fits even better: it's the gain-mask idea from earlier. The phase has to come from each trial, because induced activity has no shared phase to subtract. The ERSP supplies only the amplitude to remove.

So your scheme is: for each trial and bin, keep the trial's own phase, remove the ERSP's share of amplitude, and invert. That's sound. There are two ways to model what the event did to power:

- **Multiplicative (divide):** `W_clean = W / sqrt(gain)`. This assumes the event scales an existing rhythm, as in alpha desynchronization (ERD).
- **Additive (spectral subtraction):** `abs(W_clean)**2 = max(abs(W)**2 - excess_power, 0)`, keeping the phase. This assumes the event added a new burst on top of what was there.

```python
import numpy as np
from scipy.signal import stft, istft

class InducedRemover:
    """Removes the average event-related power change (the ERSP) from each trial, keeping each trial's phase."""

    def __init__(self, sampling_rate, window_samples=64):
        self.sampling_rate = sampling_rate
        self.window_samples = window_samples          # Hann window, 50% overlap: inverts cleanly

    def coefficients(self, trials):
        _, _, tf = stft(trials, self.sampling_rate, nperseg=self.window_samples)
        return tf                                     # shape (trials, freqs, times)

    def remove(self, trials, baseline_columns, model="multiplicative"):
        tf = self.coefficients(trials)
        power = abs(tf) ** 2
        mean_power = power.mean(axis=0)               # ERSP numerator: mean(abs(W)**2)
        baseline = mean_power[:, baseline_columns].mean(axis=1, keepdims=True)

        if model == "multiplicative":
            gain = mean_power / baseline              # ERSP as a ratio (10*log10 of this is the dB map)
            cleaned = tf / np.sqrt(gain)
        else:
            excess = np.maximum(mean_power - baseline, 0)
            cleaned_magnitude = np.sqrt(np.maximum(power - excess, 0))
            cleaned = cleaned_magnitude * np.exp(1j * np.angle(tf))   # trial's own phase kept

        _, restored = istft(cleaned, self.sampling_rate, nperseg=self.window_samples)
        return restored[..., :trials.shape[-1]]

class Demo:
    def run(self):
        rng = np.random.default_rng(0)
        rate, trials, samples = 250, 60, 500
        time = np.arange(samples) / rate
        envelope = np.exp(-0.5 * ((time - 1.2) / 0.15) ** 2)           # induced burst at 1.2 s
        phases = rng.uniform(0, 2 * np.pi, (trials, 1))                 # random phase: no template could subtract it
        data = 2 * envelope * np.cos(2 * np.pi * 10 * time + phases) + rng.normal(size=(trials, samples))

        remover = InducedRemover(rate)
        baseline_columns = slice(1, 6)                                  # STFT columns well before the burst
        for model in ["multiplicative", "additive"]:
            cleaned = remover.remove(data, baseline_columns, model)
            burst = slice(275, 325)
            print(model, "burst power before/after:", (data[:, burst] ** 2).mean().round(2), (cleaned[:, burst] ** 2).mean().round(2))

Demo().run()
```

**Is there an existing approach?** In audio, yes: spectral subtraction (Boll, 1979) and Wiener gain masks are exactly this. In EEG, the usual way to remove an induced rhythm in the time domain is spatial: ICA or GED finds the rhythm's component, and you project it out. That's the "across signals" version you suspected was harder. It is harder, but it doesn't depend on the multiplicative-versus-additive assumption.

The same rule from before applies: the ERSP works as the shared template, and only phase is free per trial. If you let each trial's amplitude be fitted freely too, everything at that bin gets removed.

Powered by Claude Exporter (https://www.ai-chat-exporter.net)