在腦電數據處理中濾波是很重要的一個步驟，直接影響后面的特征提取等計算流程。在之間寫的博客中有過介紹(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

?

總結

以上是生活随笔為你收集整理的Python中FIR滤波和STFT滤波对比（MNE脑电数据处理）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。