Python中FIR滤波和STFT滤波对比（MNE脑电数据处理）

在脑电数据处理中滤波是很重要的一个步骤，直接影响后面的特征提取等计算流程。在之间写的博客中有过介绍(https://blog.csdn.net/zhoudapeng01/article/details/106124655)，目前在脑电领域应用比较多的滤波方法有FIR,小波，以及STFT（短时傅里叶变换）等。这里主要对比MNE库提供的FIR滤波和STFT方法：

FIR滤波：FIR带通滤波在脑电数据处理中使用的非常多，其本质就是一个带通滤波器，主要用来分离不同频段的脑波数据，用于后续的数据处理工作。其在MNE库中有实现：https://blog.csdn.net/zhoudapeng01/article/details/106124655

STFT：短时傅里变换，是一种时频分析方法，网上有很多相关的介绍，其本质就是将信号按一个时间窗进行频率变换，然后堆叠在一起，可以参考下面的链接。

https://blog.csdn.net/yuelulu0629/article/details/76167229

https://zhuanlan.zhihu.com/p/150709566

博客对应的数据和代码：

https://download.csdn.net/download/zhoudapeng01/12751615

那FIR和STFT这两种方法的滤波效果有什么不同呢？下面分别对比下5种不同频段的滤波效果，代码如下。

# 短时傅里叶变换和FIR滤波效果对比import mne
import matplotlib.pyplot as plt
from scipy import signal, fft
import numpy as np
# 设置MNE库打印Log的级别
mne.set_log_level(False)
# 需要分析的频带及其范围
bandFreqs = [{'name': 'Delta', 'fmin': 1, 'fmax': 3},{'name': 'Theta', 'fmin': 4, 'fmax': 7},{'name': 'Alpha', 'fmin': 8, 'fmax': 13},{'name': 'Beta', 'fmin': 14, 'fmax': 31},{'name': 'Gamma', 'fmin': 31, 'fmax': 40}
]# 定义STFT函数
# epochsData:epochs的数据（mumpy格式）
# sfreq:采样频率
# band:频带类型
def STFT(epochsData, sfreq, band=bandFreqs):# 利用signal包进行STFT变换，f为频率，t为时间，Zxx为变换结果f, t, Zxx = signal.stft(epochsData, fs=sfreq)# 分频带保存STFT后的结果bandResult = []# 单独分析某一个频率范围for iter_freq in band:# 定位有效频率的索引index = np.where((iter_freq['fmin'] < f) & (f < iter_freq['fmax']))# 生成新的参数矩阵，初始化为复数元素为0portion = np.zeros(Zxx.shape, dtype=np.complex_)# 将有效频率赋值给新的参数矩阵portion[:, :, index, :] = Zxx[:, :, index, :]# 进行逆STFT变换，保留目标频率范围的信息_, xrec = signal.istft(portion, fs=sfreq)# 保存滤波后的结果bandResult.append(xrec)return bandResultif __name__ == '__main__':# 加载fif格式的数据epochs = mne.read_epochs(r'F:\BaiduNetdiskDownload\BCICompetition\BCICIV_2a_gdf\Train\Fif\A02T_epo.fif')# 绘图验证结果plt.figure(figsize=(15, 10))# 获取采样频率sfreq = epochs.info['sfreq']# 想要分析的目标频带bandIndex = 4# 想要分析的channelchannelIndex = 0# 想要分析的epochepochIndex = 0# 绘制原始数据plt.plot(epochs.get_data()[epochIndex][channelIndex], label='Raw')# 计算FIR滤波后的数据并绘图（注意这里要使用copy方法，否则会改变原始数据）firFilter = epochs.copy().filter(bandFreqs[bandIndex]['fmin'], bandFreqs[bandIndex]['fmax'])plt.plot(firFilter.get_data()[epochIndex][channelIndex], c=(1, 0, 0), label='FIR_Filter')# 计算STFT滤波后的数据并绘图stft = STFT(epochs.get_data(), sfreq)plt.plot(stft[bandIndex][epochIndex][channelIndex], c=(0, 1, 0), label='STFT_Filter')# 绘制图例和图名plt.legend()plt.title(bandFreqs[bandIndex]['name'])####################################FFT对比两种方法的频谱分布plt.figure(figsize=(15, 10))# 对FIR滤波后的数据进行FFT变换mneFIRFreq = np.abs(fft.fft(firFilter.get_data()[epochIndex][channelIndex]))# 对STFT滤波后的数据进行FFT变换，需要注意STFT变换后数据的点数可能会发生变化，这里截取数据保持一致性pointNum = epochs.get_data()[epochIndex][channelIndex].shape[0]stftFreq = np.abs(fft.fft(stft[bandIndex][epochIndex][channelIndex][:pointNum]))# 想要绘制的点数pointPlot = 500# FIR滤波后x轴对应的频率幅值范围FIR_X = np.linspace(0, sfreq/2, int(mneFIRFreq.shape[0]/2))# STFT滤波后x轴对应的频率幅值范围STFT_X = np.linspace(0, sfreq/2, int(stftFreq.shape[0]/2))# 绘制FIR滤波后的频谱分布plt.plot(FIR_X[:pointPlot], mneFIRFreq[:pointPlot], c=(1, 0, 0), label='FIR_Filter')# 绘制STFT滤波后的频谱分布plt.plot(STFT_X[:pointPlot], stftFreq[:pointPlot], c=(0, 1, 0), label='STFT_FIlter')# 绘制图例和图名plt.legend()plt.title(bandFreqs[bandIndex]['name'])plt.show()

注：STFT滤波后数据长度发生改变，这主要和窗长及计算方式有关，超出原始长度的数据可以不用关心，上面的代码中在进行频谱分析前，也就是计算FFT前对数据的长度进行了处理，这样可以保证分析出频谱数据的一致性。

FIR和STFT滤波在不同频段的效果对比

时域对比：从时域的滤波结果来看，FIR和STFT的趋势基本保持一致，只是其幅值会有些差别。

频域对比：从频谱分析的结果来看，STFT滤波后信号的频带分布范围更加准确，滤波后信号的频谱分布看起来更加符合预期，如下图所示，FIR滤波后的信号频谱分布较广，甚至超出了目标范围。

Delta频段（1Hz-3Hz）对应的结果

Theta频段（4Hz-7Hz）对应的结果

Alpha频段（8Hz-13Hz）对应的结果

Beta频段（14Hz-31Hz）对应的结果

Gamma频段（31Hz-40Hz）对应的结果

博客对应的数据和代码：

https://download.csdn.net/download/zhoudapeng01/12751615