用C# MathNet复现Python spectrogram结果：排查代码偏差原因

我需要用C#的MathNet库复现Python中scipy.signal.spectrogram的计算结果，当前C#结果必须乘以1.5才能和Python结果匹配，希望找出代码错误并修正，避免使用这个魔法数。

Python参考代码 #

import math
from scipy import signal

# Sampling rate [Hz]
fs = 600.0

# FFT window length [sec].
win_sec = 3.0

# Upper frequency [Hz]
fq_lim = 50

# Sampling interval [sec]
dt = 1. / fs

# Data length [points]
n = int(fs * win_sec)

# Power of 2 greater than argument calculated
nfft = 2 ** math.ceil(math.log2(n))

# Frequency signal interval [/sec].
df_f = 1.0 / (nfft * dt)

# Spectrogram calculation Frequency (freq) Time (t) Power spectral density PSD (μV²/Hz) calculation
freq, t, spctgram = signal.spectrogram(sig, fs=fs, nperseg=n, nfft=nfft, noverlap=(win_sec-1)*fs, window='hamming', mode='psd', return_onesided=True)

C#待修正代码 #

public Dictionary<string, float> Analyze(List<float> brainwaveData)
{
    // power of 2 greater than the sampling number
    var fftWindowSize = (int)Math.Pow(2, Math.Ceiling(Math.Log(Const.NumOfSampling, 2)));

    // Fill with zeros until sampling number is next power of 2
    while (brainwaveData.Count < fftWindowSize)
    {
        brainwaveData.Add(0f);
    }

    // Apply Hamming window for FFT
    var coefficientHammingWindow = 0.54f;
    double[] window = Window.HammingPeriodic(fftWindowSize);
    var signalWindow = brainwaveData.Select((d, i) => new Complex32((float)(d * window[i] / coefficientHammingWindow), 0)).ToArray();

    // Execute FFT
    Fourier.Forward(signalWindow, FourierOptions.Matlab);

    // Power spectrum [μV²] calculation Amplitude squared
    var powerSpectrum = signalWindow.Select(d => (d.Real * d.Real) + (d.Imaginary * d.Imaginary)).ToArray();

    // power spectrum normalization [μV²]
    var normalizedPowerSpectrum = powerSpectrum.Select(d => d / fftWindowSize).ToArray();

    // Power spectral density [μV²/Hz] calculation
    var powerSpectrumDensity = normalizedPowerSpectrum.Select(d => d / Const.SamplingFrequencyHz).ToArray();

    // Calculated for each EEG component
    Dictionary<string, float> result = new Dictionary<string, float>();
    foreach (var listBand in listBands)
    {
        result.Add(listBand.BandName, calcFrequencyBandAverage(powerSpectrumDensity, Const.SamplingFrequencyHz, listBand.MinFreqHz, listBand.MaxFreqHz));
    }
    return result;
}

private float calcFrequencyBandAverage(float[] spectrum, int sampleRate, double minHz, double maxHz)
{
    // Find the index for a given frequency band (Hz)
    var minIndex = (int)Math.Ceiling(spectrum.Length * minHz / sampleRate);
    var maxIndex = (int)Math.Floor(spectrum.Length * maxHz / sampleRate);

    // Calculate average
    return spectrum.Skip(minIndex).Take(maxIndex - minIndex).Average();
}

List<BrainwaveFrequencyBands> listBands = new List<BrainwaveFrequencyBands>
{
    new BrainwaveFrequencyBands{ BandName = "Delta", MinFreqHz = 1f, MaxFreqHz = 4f },
    new BrainwaveFrequencyBands{ BandName = "Theta", MinFreqHz = 4f, MaxFreqHz = 8f },
    new BrainwaveFrequencyBands{ BandName = "Alpha1", MinFreqHz = 8f, MaxFreqHz = 10f },
    new BrainwaveFrequencyBands{ BandName = "Alpha2", MinFreqHz = 10f, MaxFreqHz = 12f },
    new BrainwaveFrequencyBands{ BandName = "Beta1", MinFreqHz = 12f, MaxFreqHz = 20f },
    new BrainwaveFrequencyBands{ BandName = "Beta2", MinFreqHz = 20f, MaxFreqHz = 30f },
    new BrainwaveFrequencyBands{ BandName = "Gamma1", MinFreqHz = 30f, MaxFreqHz = 40f },
    new BrainwaveFrequencyBands{ BandName = "Gamma2", MinFreqHz = 40f, MaxFreqHz = 50f },
};

多维度问题分析 #

C#代码核心错误点 #

1. 窗函数错误归一化：MathNet的Window.HammingPeriodic生成的是标准汉明窗（公式：0.54 - 0.46*cos(2πn/N)），代码中额外除以0.54会放大窗的能量，导致后续PSD结果偏高。
2. 数据处理顺序颠倒：Python是先对nperseg长度的分段数据加窗，再补零到nfft；而C#是先补零到fftWindowSize再加窗，改变了加窗的有效数据范围。
3. PSD归一化逻辑缺失：
   • 未考虑窗能量校正因子：scipy会自动用sum(window²)校正窗带来的能量衰减，C#代码完全忽略这一步。
   • 未处理单边谱的能量合并：scipy在return_onesided=True时，会将双边谱的非直流/奈奎斯特频率分量能量加倍，C#代码没有对应操作。
4. 频率索引计算错误：基于全长度谱计算索引，但实际处理的是单边谱，导致频段划分偏移。

Python实现逻辑对照 #

scipy.signal.spectrogram的PSD计算核心公式为：

psd = |FFT(x * window)|² / (fs * sum(window²))

当return_onesided=True时，对实信号的正频率部分（除0和fs/2）乘以2，将双边谱的能量合并到单边，保证能量守恒。

EEG信号处理的一致性要求 #

脑电PSD的物理意义是单位频率的能量（μV²/Hz），窗归一化和能量校正错误会导致结果偏离真实脑电能量分布，影响后续频段分析的可靠性。

修正后的C#代码 #

public Dictionary<string, float> Analyze(List<float> brainwaveData)
{
    int nperseg = Const.NumOfSampling; // 对应Python的n = fs*win_sec
    int fftWindowSize = (int)Math.Pow(2, Math.Ceiling(Math.Log(nperseg, 2)));

    // 1. 确保输入数据长度为nperseg（匹配Python的分段长度）
    if (brainwaveData.Count > nperseg)
    {
        brainwaveData = brainwaveData.Take(nperseg).ToList();
    }
    else while (brainwaveData.Count < nperseg)
    {
        brainwaveData.Add(0f);
    }

    // 2. 应用标准汉明窗（无需额外归一化）
    double[] window = Window.HammingPeriodic(nperseg);
    var signalWindow = brainwaveData.Select((d, i) => new Complex32((float)(d * window[i]), 0)).ToArray();

    // 3. 补零到fftWindowSize长度（对应Python的nfft）
    Array.Resize(ref signalWindow, fftWindowSize);

    // 4. 执行FFT（Matlab选项对齐scipy的FFT输出）
    Fourier.Forward(signalWindow, FourierOptions.Matlab);

    // 5. 计算功率谱（FFT幅值平方）
    var powerSpectrum = signalWindow.Select(d => (d.Real * d.Real) + (d.Imaginary * d.Imaginary)).ToArray();

    // 6. 计算窗的能量和，用于PSD校正
    double windowEnergy = window.Sum(x => x * x);

    // 7. 生成单边PSD，匹配scipy的return_onesided=True逻辑
    var powerSpectrumDensity = new float[fftWindowSize / 2 + 1];
    for (int i = 0; i < powerSpectrumDensity.Length; i++)
    {
        double psdVal = powerSpectrum[i];
        // 单边谱能量合并：非直流/奈奎斯特频率分量乘以2
        if (i != 0 && i != fftWindowSize / 2)
        {
            psdVal *= 2;
        }
        // 应用scipy的PSD公式：校正采样率和窗能量
        psdVal /= (Const.SamplingFrequencyHz * windowEnergy);
        powerSpectrumDensity[i] = (float)psdVal;
    }

    // 8. 计算各频段平均功率
    Dictionary<string, float> result = new Dictionary<string, float>();
    foreach (var band in listBands)
    {
        result.Add(band.BandName, CalculateBandAverage(powerSpectrumDensity, Const.SamplingFrequencyHz, band.MinFreqHz, band.MaxFreqHz));
    }
    return result;
}

private float CalculateBandAverage(float[] spectrum, int sampleRate, double minHz, double maxHz)
{
    // 单边谱的频率分辨率：df = sampleRate / fftWindowSize
    double fftSize = (spectrum.Length - 1) * 2; // 还原完整FFT长度
    double df = (double)sampleRate / fftSize;

    // 计算频段对应的索引
    int minIdx = (int)Math.Ceiling(minHz / df);
    int maxIdx = (int)Math.Floor(maxHz / df);

    // 边界保护，避免索引越界
    minIdx = Math.Max(0, minIdx);
    maxIdx = Math.Min(spectrum.Length - 1, maxIdx);

    if (minIdx >= maxIdx)
    {
        return 0f;
    }

    return spectrum.Skip(minIdx).Take(maxIdx - minIdx).Average();
}

List<BrainwaveFrequencyBands> listBands = new List<BrainwaveFrequencyBands>
{
    new BrainwaveFrequencyBands{ BandName = "Delta", MinFreqHz = 1f, MaxFreqHz = 4f },
    new BrainwaveFrequencyBands{ BandName = "Theta", MinFreqHz = 4f, MaxFreqHz = 8f },
    new BrainwaveFrequencyBands{ BandName = "Alpha1", MinFreqHz = 8f, MaxFreqHz = 10f },
    new BrainwaveFrequencyBands{ BandName = "Alpha2", MinFreqHz = 10f, MaxFreqHz = 12f },
    new BrainwaveFrequencyBands{ BandName = "Beta1", MinFreqHz = 12f, MaxFreqHz = 20f },
    new BrainwaveFrequencyBands{ BandName = "Beta2", MinFreqHz = 20f, MaxFreqHz = 30f },
    new BrainwaveFrequencyBands{ BandName = "Gamma1", MinFreqHz = 30f, MaxFreqHz = 40f },
    new BrainwaveFrequencyBands{ BandName = "Gamma2", MinFreqHz = 40f, MaxFreqHz = 50f },
};

修正说明 #

1. 调整数据流程：先取nperseg长度数据，加窗后再补零，和Python处理顺序一致。
2. 移除错误窗归一化：删除/ coefficientHammingWindow操作，使用标准汉明窗。
3. 添加窗能量校正：用sum(window²)修正窗带来的能量衰减，匹配scipy的归一化逻辑。
4. 单边谱校正：对非直流/奈奎斯特频率分量乘以2，保证能量守恒。
5. 修正频率索引：基于单边谱的频率分辨率计算索引，确保频段划分准确。

内容的提问来源于stack exchange，提问作者Ganessa