`output="complex"` for Welch method should only allow `average=None`

已关闭 适合新手
  • 哪里不对？
#14,495 1 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

维护者通常 1 天内回复

已经有一个关联 PR 被合并了。

  • #14498 来自 @iis-MarkKuang —— 已合并

评估

难度
2/5
预计耗时
1-3 小时
新手友好度
74/100
Issue 类型
缺陷
描述清晰度
描述清楚
活跃度
活跃
技术栈
python
领域
data

调研方向

首先追踪 epochs.compute_psd() 中 Welch 对 output="complex" 和 average 的处理，然后检查 Spectrum.plot() 如何将系数转换为功率。添加检查，在复数输出时拒绝 None 以外的 average 值，并查找覆盖 Welch 输出模式和绘图的现有测试。无效组合会引发明确错误，且 average=None 和功率输出保持预期行为时，即为完成。

由索引模型根据 Issue 内容生成。

描述

BUG
Description of the problem

For Welch method, if average is requested with output="complex", segments are aggregated for complex coefficients, which can mess up the phase information in the resulting coefficients. Subsequently computing power from these coefficients gives very different results to if output="power" was used directly. This also affects the conversion to power from coefficients in Spectrum.plot().

Steps to reproduce
import numpy as np

import mne
from mne.time_frequency.psd import _median_biases

data_path = mne.datasets.sample.data_path()
raw = mne.io.read_raw_fif(
    data_path / "MEG" / "sample" / "sample_audvis_filt-0-40_raw.fif", verbose=False
)
raw.pick("eeg", exclude="bads").crop(tmax=60).load_data(verbose=False)
epochs = mne.make_fixed_length_epochs(raw, duration=4.0, preload=True, verbose=False)
kwargs = dict(method="welch", n_fft=128, verbose=False)  # 4 segments per epoch


def coefs_to_power(spectrum):
    """Convert complex Welch coefficients to one-sided power."""
    coefs = spectrum.get_data()
    freqs, sfreq = spectrum.freqs, spectrum.info["sfreq"]
    # double all bins except DC and Nyquist (scipy.signal.spectrogram does this for
    # mode="psd", but not for mode="complex")
    scale = np.where((freqs == 0) | np.isclose(freqs, sfreq / 2), 1.0, 2.0)
    if coefs.ndim == 4:  # unaggregated: (epoch, channel, freq, segment)
        scale = scale[:, np.newaxis]
    return (coefs * coefs.conj()).real * scale


for average in ("mean", "median"):
    power = epochs.compute_psd(output="power", average=average, **kwargs).get_data()

    # output="complex" with average!=None aggregates the complex coefficients across
    # segments, so power computed afterwards is wrong
    spec = epochs.compute_psd(output="complex", average=average, **kwargs)
    power_agg_first = coefs_to_power(spec)

    # correct: average=None, compute power per segment, then aggregate
    spec = epochs.compute_psd(output="complex", average=None, **kwargs)
    power_seg = coefs_to_power(spec)
    if average == "mean":
        power_agg_last = power_seg.mean(axis=-1)
    else:
        n_segments = power_seg.shape[-1]
        power_agg_last = (
            np.median(power_seg, axis=-1) / _median_biases(n_segments)[n_segments]
        )

    print(f"average={average!r}")
    for name, psd in (
        ("complex -> aggregate -> power", power_agg_first),
        ("complex -> power -> aggregate", power_agg_last),
    ):
        ratio = psd / power
        print(
            f"  {name}: ratio to output='power' median={np.median(ratio):.3f} "
            f"min={ratio.min():.3g} max={ratio.max():.3g} "
            f"allclose={np.allclose(psd, power, rtol=1e-6, atol=0)}"
        )
Link to data

No response

Expected results

average!=None with output="complex" is not allowed, as it gives invalid coefficients.

Actual results

Computing output="complex" with average!=None is allowed, and gives very different results when trying to compute power from these coefficients, compared to the proper approach of unaggregated coefficients → unaggregated power → aggregated power (which is what happens with output="power").

average='mean'
  complex -> aggregate -> power: ratio to output='power' median=0.131 min=2.06e-08 max=0.992 allclose=False
  complex -> power -> aggregate: ratio to output='power' median=1.000 min=1 max=1 allclose=True
average='median'
  complex -> aggregate -> power: ratio to output='power' median=0.326 min=9.42e-08 max=36.2 allclose=False
  complex -> power -> aggregate: ratio to output='power' median=1.000 min=1 max=1 allclose=True
Additional information

AI disclosure: Claude Opus 5.5 initially identified the problem, which I further investigated.

主要语言
Python
星标
3.6k
派生
1.6k
平均合并
1 天 18 小时
30 天内合并 PR
71

环境准备

从这里开始

  1. 先读完整个 Issue，再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库，在一个分支上完成修改。
  4. 提交 Pull Request，并在描述里引用这个 Issue 编号。

mne-tools/mne-python 的其他 Issue

查看 mne-tools/mne-python 的全部 Issue

相似的 Issue

更多 Python Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。