原文作者：
我们将了解如何使用Welch 和多锥谱估计方法计算特定频率范围内信号的平均功率。本教程主要面向具有一些 EEG 信号处理基础知识的神经科学家/睡眠研究人员。
前言
分析 EEG 数据最广泛使用的方法之一是将信号分解为功能上不同的频带，例如delta（0.5-4 Hz）、 theta（4-8 Hz）、alpha（8-12 Hz）、beta （12-30 Hz）和gamma（30-100 Hz）。
这意味着将 EEG 信号分解为频率分量，这通常是通过傅里叶变换实现的。计算傅里叶变换的几乎总是使用的算法（可以说是最重要的信号处理算法）是快速傅里叶变换 (FFT)，它为每个频率点返回一个复数，然后可以从中轻松提取该特定频率下信号的幅度和相位。在频谱分析中，通常取 FFT 的幅度平方以获得功率谱密度（或功率谱或周期图）的估计值，以（微）伏2表示对于 EEG 数据来说，每赫兹。
虽然可以从功率谱密度进行无数种分析，但我将重点介绍一种非常简单的分析：平均频带功率，即计算一个数字，总结给定频带对信号总功率的贡献。这在机器学习方法中可能特别有用，因为您经常需要从数据中提取一些关键特征（并得到一个可以总结数据特定方面的数字）。
平均频带功率也是睡眠研究的一个非常相关的指标，因为它可以区分不同的睡眠阶段。例如，深度睡眠其特点是慢波占主导地位，频率范围在 0.5 到 4 Hz 之间（即 delta 波段），反映了同步的大脑活动。相反，清醒状态的特点是 delta 活动很少，高频活动更多。因此，如果您要计算深度睡眠和清醒状态的 delta 波段功率，前者会非常高，而后者会非常低。请自行查看以下内容：
睡眠和清醒状态下的典型大脑 (Fz 和 Cz)、眼部 (EOG) 和肌肉 (EMG) 活动。每个框代表同一个人的 10 秒数据。请注意深度睡眠 (N3) 主要由高振幅和低频振荡主导，而清醒状态主要由低振幅和高频振荡主导。
数据加载
为了本教程的目的，请查看以下来自一位年轻人的真实慢波睡眠的 30 秒摘录。采样频率为 100 Hz，通道为 F3。
是时候打开你最喜欢的 Python 编辑器了！如果你是 Python 新手，我强烈建议你使用Jupyter Lab。加载数据相当容易：
import numpy as np
data = np.loadtxt('data.txt')我们来看一下数据：
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.2)
# Define sampling frequency and time vector
sf = 100.
time = np.arange(data.size) / sf
# Plot the signal
fig, ax = plt.subplots(1, 1, figsize=(12, 4))
plt.plot(time, data, lw=1.5, color='k')
plt.xlabel('Time (seconds)')
plt.ylabel('Voltage')
plt.xlim([time.min(), time.max()])
plt.title('N3 sleep EEG data (F3)')
sns.despine()
计算功率谱密度
为了计算 delta 频带中的平均频带功率，我们首先需要计算功率谱密度的估计值。最广泛使用的方法是韦尔奇周期图，它包括对信号小窗口的连续傅里叶变换进行平均，有重叠或无重叠。
Welch 的方法提高了经典周期图的准确性。原因很简单：EEG 数据总是随时间变化的，这意味着如果你查看 30 秒的 EEG 数据，信号看起来非常（非常）不可能像纯正弦波的完美总和。相反，EEG 的频谱内容会随时间而变化，并不断受到头皮下神经元活动的改变。问题是，要返回真实的频谱估计值，经典周期图要求信号的频谱内容在考虑的时间段内保持平稳（即不随时间变化）。由于情况从来都不是这样的，因此周期图通常存在偏差，并且包含太多方差（请参阅本教程的结尾）。通过对在窗口的短段上获得的周期图进行平均，Welch 方法可以大幅减少这种方差。然而，这是以较低的频率分辨率为代价的。事实上，频率分辨率定义为：
Fs是信号的采样频率，N样本总数和，t信号的持续时间（以秒为单位）。换句话说，如果我们要使用整个数据长度（30 秒），我们最终的频率分辨率将是 1/30=0.033 Hz，即每赫兹 30 个频率点。通过使用 4 秒滑动窗口，我们将此频率分辨率降低到每赫兹 4 个频率点，即每步代表 0.25 Hz。
- 从上面的公式可以得出另一个结论：唯一能提高频率分辨率的就是时间。采样频率的变化不会提高频率分辨率，只会提高频率覆盖范围。点击此处了解更多信息。
from scipy import signal
# Define window length (4 seconds)
win = 4 * sf
freqs, psd = signal.welch(data, sf, nperseg=win)
# Plot the power spectrum
sns.set(font_scale=1.2, style='white')
plt.figure(figsize=(8, 4))
plt.plot(freqs, psd, color='k', lw=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (V^2 / Hz)')
plt.ylim([0, psd.max() * 1.1])
plt.title("Welch's periodogram")
plt.xlim([0, freqs.max()])
sns.despine()
该freqs向量包含 x 轴（频率区间），而该psd向量包含 y 轴（功率谱密度）。处理 EEG 数据时，功率谱密度的单位通常是微伏平方每赫兹（μv / Hz）
- 请注意，x 轴的最大值始终是采样频率的一半，这正好是奈奎斯特频率。这就是频率覆盖概念发挥作用的地方：如果我们的信号以 200 Hz 而不是 100 Hz 采样，则 x 轴上的最大值将是 200 / 2 = 100 Hz，而不是 50 Hz。换句话说，增加采样频率会导致更大的频率范围。
定义δ（delta）频带
现在，在计算平均delta带宽功率之前，我们需要找到与delta频率范围相交的频率箱。
# Define delta lower and upper limits
low, high = 0.5, 4
# Find intersecting values in frequency vector
idx_delta = np.logical_and(freqs >= low, freqs <= high)
# Plot the power spectral density and fill the delta area
plt.figure(figsize=(7, 4))
plt.plot(freqs, psd, lw=2, color='k')
plt.fill_between(freqs, psd, where=idx_delta, color='skyblue')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (uV^2 / Hz)')
plt.xlim([0, 10])
plt.ylim([0, psd.max() * 1.1])
plt.title("Welch's periodogram")
sns.despine()
平均频带功率
绝对 delta 功率等于上一个图的蓝色区域。由于没有闭式公式来积分这个区域，我们需要对其进行近似。这通常使用复合辛普森规则来实现。其背后的想法其实非常简单：我们将这个区域分解为几个抛物线，然后对这些抛物线的面积求和。请注意，这也可以使用梯形（梯形规则）或矩形（如 Matlab 的bandpower函数）来完成，但抛物线通常会给出更好的估计。
from scipy.integrate import simps
# Frequency resolution
freq_res = freqs[1] - freqs[0] # = 1 / 4 = 0.25
# Compute the absolute power by approximating the area under the curve
delta_power = simps(psd[idx_delta], dx=freq_res)
print('Absolute delta power: %.3f uV^2' % delta_power)Absolute delta power: 321.064 uV^2
相对功率
在实践中，人们可能不想报告绝对频带功率，而是希望将频带中的功率表示为信号总功率的百分比。这称为相对频带功率。它可以很容易地从上式计算出来：
total_power = simps(psd, dx=freq_res)
delta_rel_power = delta_power / total_power
print('Relative delta power: %.3f' % delta_rel_power)Relative delta power: 0.787
换句话说，78.7% 是delta频带功率中。
概括
下面的函数是上述函数的概括，可用于轻松获取特定频带中的平均绝对功率或相对功率。它与 Matlab bandpower函数非常相似，不同之处在于它使用 Welch 周期图而不是经典周期图，并且它使用抛物线而不是矩形来近似面积。
- 在专用于睡眠分析的 Python 包YASA中实现了该函数的更高级、更灵活的版本。除其他功能外， yasa.bandpower()函数可以直接处理多个通道和多个光谱带。
def bandpower(data, sf, band, window_sec=None, relative=False):
"""Compute the average power of the signal x in a specific frequency band.
Parameters
----------
data : 1d-array
Input signal in the time-domain.
sf : float
Sampling frequency of the data.
band : list
Lower and upper frequencies of the band of interest.
window_sec : float
Length of each window in seconds.
If None, window_sec = (1 / min(band)) * 2
relative : boolean
If True, return the relative power (= divided by the total power of the signal).
If False (default), return the absolute power.
Return
------
bp : float
Absolute or relative band power.
"""
from scipy.signal import welch
from scipy.integrate import simps
band = np.asarray(band)
low, high = band
# Define window length
if window_sec is not None:
nperseg = window_sec * sf
else:
nperseg = (2 / low) * sf
# Compute the modified periodogram (Welch)
freqs, psd = welch(data, sf, nperseg=nperseg)
# Frequency resolution
freq_res = freqs[1] - freqs[0]
# Find closest indices of band in frequency vector
idx_band = np.logical_and(freqs >= low, freqs <= high)
# Integral approximation of the spectrum using Simpson's rule.
bp = simps(psd[idx_band], dx=freq_res)
if relative:
bp /= simps(psd, dx=freq_res)
return bp两个频带之间的比率
报告两个频带之间的比率也很常见。例如，delta/beta 比率是慢波睡眠质量的著名指标。在计算两个频带之间的比率时，重要的是控制两个频带的周期图窗口长度相同。事实上，如果您对两个频带使用不同的窗口长度，这将导致两个不同的周期图，因此比率将毫无意义。
# Define the duration of the window to be 4 seconds
win_sec = 4
# Delta/beta ratio based on the absolute power
db = bandpower(data, sf, [0.5, 4], win_sec) / bandpower(data, sf, [12, 30], win_sec)
# Delta/beta ratio based on the relative power
db_rel = bandpower(data, sf, [0.5, 4], win_sec, True) / bandpower(data, sf, [12, 30], win_sec, True)
print('Delta/beta ratio (absolute): %.3f' % db)
print('Delta/beta ratio (relative): %.3f' % db_rel)Delta/beta ratio (absolute): 42.214
Delta/beta ratio (relative): 42.214
使用多锥度方法
Multitaper 是一种频谱分析方法，由 David J. Thompson 于 1982 年首次开发，旨在克服传统频谱估计技术的一些局限性。它结合了两种方法的优点：高频分辨率和低方差，提供了比传统和 Welch 周期图更稳健的频谱估计。
要了解它的工作原理，我强烈建议您阅读Prerau 等人 (2017) 的这篇论文，该论文很好地解释了什么是 Multitaper，以及睡眠研究和睡眠医学如何从中受益。对我有用的另一个来源是关于多锥光谱分析的 Matlab 文档。
简而言之，多锥化方法首先使用一组最佳带通滤波器（称为斯莱皮序列 (DPSS)）过滤原始信号。此过滤是通过将原始信号与斯莱皮序列进行卷积来完成的。其次，为每个新的过滤（或“锥形”）数据计算经典周期图，然后通过对所有得到的周期图取平均值来获得最终频谱。多锥化方法的真正优势在于斯莱皮序列与所有其他序列正交，因此锥形信号提供了底层频谱的统计独立估计。换句话说，每个经过过滤的信号都会突出显示信号频谱内容的一个特定方面。
使用 Multitaper 的平均带宽功率
多锥谱估计方法在MNE-Python 包中实现。在下面的例子中，我调整了我们之前创建的 bandpower 函数以添加多锥方法。
def bandpower(data, sf, band, method='welch', window_sec=None, relative=False):
"""Compute the average power of the signal x in a specific frequency band.
Requires MNE-Python >= 0.14.
Parameters
----------
data : 1d-array
Input signal in the time-domain.
sf : float
Sampling frequency of the data.
band : list
Lower and upper frequencies of the band of interest.
method : string
Periodogram method: 'welch' or 'multitaper'
window_sec : float
Length of each window in seconds. Useful only if method == 'welch'.
If None, window_sec = (1 / min(band)) * 2.
relative : boolean
If True, return the relative power (= divided by the total power of the signal).
If False (default), return the absolute power.
Return
------
bp : float
Absolute or relative band power.
"""
from scipy.signal import welch
from scipy.integrate import simps
from mne.time_frequency import psd_array_multitaper
band = np.asarray(band)
low, high = band
# Compute the modified periodogram (Welch)
if method == 'welch':
if window_sec is not None:
nperseg = window_sec * sf
else:
nperseg = (2 / low) * sf
freqs, psd = welch(data, sf, nperseg=nperseg)
elif method == 'multitaper':
psd, freqs = psd_array_multitaper(data, sf, adaptive=True,
normalization='full', verbose=0)
# Frequency resolution
freq_res = freqs[1] - freqs[0]
# Find index of band in frequency vector
idx_band = np.logical_and(freqs >= low, freqs <= high)
# Integral approximation of the spectrum using parabola (Simpson's rule)
bp = simps(psd[idx_band], dx=freq_res)
if relative:
bp /= simps(psd, dx=freq_res)
return bp让我们用下面的代码尝试一下我们的新函数。与 Welch 方法相比，Multitaper 方法的一个优点是我们不需要指定窗口持续时间，因为它基本上会计算整个信号的周期图。由于我们使用信号的整个长度，因此多锥估计的频率分辨率将为 1 / 30 = 0.033 Hz。
# Multitaper delta power
bp = bandpower(data, sf, [0.5, 4], 'multitaper')
bp_rel = bandpower(data, sf, [0.5, 4], 'multitaper', relative=True)
print('Absolute delta power: %.3f' % bp)
print('Relative delta power: %.3f' % bp_rel)
# Delta-beta ratio
# One advantage of the multitaper is that we don't need to define a window length.
db = bandpower(data, sf, [0.5, 4], 'multitaper') / bandpower(data, sf, [12, 30], 'multitaper')
# Ratio based on the relative power
db_rel = bandpower(data, sf, [0.5, 4], 'multitaper', relative=True) / \
bandpower(data, sf, [12, 30], 'multitaper', relative=True)
print('Delta/beta ratio (absolute): %.3f' % db)
print('Delta/beta ratio (relative): %.3f' % db_rel)
绝对 delta 功率: 311.559
相对 delta 功率: 0.790
Delta/beta 比率 (绝对) : 41.225
Delta/beta 比率 (相对) : 41.225
结果与使用 Welch 方法获得的结果非常接近。只要您的数据不太嘈杂，这应该是正确的。但是，如果您处理的是嘈杂的数据，多锥体方法将始终提供比 Welch 方法更稳健的频谱估计。
只是为了好玩，让我们比较一下使用经典周期图、Welch 周期图和多锥体方法获得的功率谱密度估计：
def plot_spectrum_methods(data, sf, window_sec, band=None, dB=False):
"""Plot the periodogram, Welch's and multitaper PSD.
Requires MNE-Python >= 0.14.
Parameters
----------
data : 1d-array
Input signal in the time-domain.
sf : float
Sampling frequency of the data.
band : list
Lower and upper frequencies of the band of interest.
window_sec : float
Length of each window in seconds for Welch's PSD
dB : boolean
If True, convert the power to dB.
"""
from mne.time_frequency import psd_array_multitaper
from scipy.signal import welch, periodogram
sns.set(style="white", font_scale=1.2)
# Compute the PSD
freqs, psd = periodogram(data, sf)
freqs_welch, psd_welch = welch(data, sf, nperseg=window_sec*sf)
psd_mt, freqs_mt = psd_array_multitaper(data, sf, adaptive=True,
normalization='full', verbose=0)
sharey = False
# Optional: convert power to decibels (dB = 10 * log10(power))
if dB:
psd = 10 * np.log10(psd)
psd_welch = 10 * np.log10(psd_welch)
psd_mt = 10 * np.log10(psd_mt)
sharey = True
# Start plot
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=sharey)
# Stem
sc = 'slategrey'
ax1.stem(freqs, psd, linefmt=sc, basefmt=" ", markerfmt=" ")
ax2.stem(freqs_welch, psd_welch, linefmt=sc, basefmt=" ", markerfmt=" ")
ax3.stem(freqs_mt, psd_mt, linefmt=sc, basefmt=" ", markerfmt=" ")
# Line
lc, lw = 'k', 2
ax1.plot(freqs, psd, lw=lw, color=lc)
ax2.plot(freqs_welch, psd_welch, lw=lw, color=lc)
ax3.plot(freqs_mt, psd_mt, lw=lw, color=lc)
# Labels and axes
ax1.set_xlabel('Frequency (Hz)')
if not dB:
ax1.set_ylabel('Power spectral density (V^2/Hz)')
else:
ax1.set_ylabel('Decibels (dB / Hz)')
ax1.set_title('Periodogram')
ax2.set_title('Welch')
ax3.set_title('Multitaper')
if band is not None:
ax1.set_xlim(band)
ax1.set_ylim(ymin=0)
ax2.set_ylim(ymin=0)
ax3.set_ylim(ymin=0)
sns.despine()
# Example: plot the 0.5 - 2 Hz band
plot_spectrum_methods(data, sf, 4, [0.5, 2], dB=True) 
- 经典周期图具有良好的频率分辨率（一个频率箱 = 0.033 Hz），但方差太大。
- Welch 周期图的方差较低，但频率分辨率较低（一个频率箱 = 0.25 Hz）。
- 多锥周期图兼具前两种方法的优点：高频率分辨率和低方差。
%timeit bandpower(data, sf, [0.5, 4], method="welch", relative=True) %timeit bandpower(data, sf, [0.5, 4], method="multitaper", relative=True)
每循环 348 µs ± 18.3 µs（7 次运行的平均值 ± 标准差，每次 1000 次循环）
每循环 108 ms ± 2.38 ms（7 次运行的平均值 ± 标准差，每次 10 次循环）
多锥体法比其他方法更耗费计算资源。在上面的例子中，多锥体法比 Welch 法慢 300 倍。使用 100 Hz 采样的 30 秒数据，执行时间差异约为 ~100ms。现在，假设您有几个小时的数据和几个以 1000 Hz 采样的通道；在我们的 30 秒数据上几乎察觉不到的东西可能会变成几个小时的差异（例如，Welch 为 1 分钟，多锥体法为 ~5 小时）。在决定使用哪种方法时，您肯定需要考虑这一点。
其次，尽管多锥体方法总是能提供更稳健的光谱估计，但我认为选择使用哪种技术取决于手头的数据。例如，如果你有从年轻健康个体获得的干净数据，那么多锥体光谱估计与韦尔奇估计相差不大的可能性很大。此外，韦尔奇方法可能是迄今为止使用最广泛的光谱估计技术，并且非常直观易懂。相比之下，多锥体是一种相对较新的方法，在概念上更难掌握。
301
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
