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

程序员文章站 2022-09-21 08:39:22

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

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

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

博客对应的数据和代码：

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']))
        # 生成新的参数矩阵，初始化为复数元素为0
        portion = np.zeros(Zxx.shape, dtype=np.complex_)
        # 将有效频率赋值给新的参数矩阵
        portion[:, :, index, :] = Zxx[:, :, index, :]
        # 进行逆STFT变换，保留目标频率范围的信息
        _, xrec = signal.istft(portion, fs=sfreq)
        # 保存滤波后的结果
        bandResult.append(xrec)
    return bandResult



if __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
    # 想要分析的channel
    channelIndex = 0
    # 想要分析的epoch
    epochIndex = 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前对数据的长度进行了处理，这样可以保证分析出频谱数据的一致性。

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