小波變換有信號顯微鏡之稱，在EEG分析中也有廣泛的應用，印象中小波算法是來源于地球物理解釋的。之前有介紹過小波的一些資料和實現：https://blog.csdn.net/zhoudapeng01/article/details/107025901可以參考下，這里主要分析小波和FIR濾波效果的對比。

博客對應的代碼和數據

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

# 短時傅里葉變換和FIR濾波效果對比import mne import matplotlib.pyplot as plt from scipy import signal, fft import numpy as np import pywt # 設置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} ]def __CalcWP(data, sfreq, wavelet, maxlevel, band):# 如果maxlevel太小部分波段分析不到wp = pywt.WaveletPacket(data=data, wavelet=wavelet, mode='symmetric', maxlevel=maxlevel)# 頻譜由低到高的對應關系，這里需要注意小波變換的頻帶排列默認并不是順序排列，所以這里需要使用’freq‘排序。freqTree = [node.path for node in wp.get_level(maxlevel, 'freq')]# 計算maxlevel最小頻段的帶寬，采樣頻率的一半freqBand = (sfreq/2) / (2 ** maxlevel)bandResult = []#######################根據實際情況計算頻譜對應關系，這里要注意系數的順序for iter_freq in band:# 構造空的小波包new_wp = pywt.WaveletPacket(data=None, wavelet=wavelet, mode='symmetric', maxlevel=maxlevel)for i in range(len(freqTree)):# 第i個頻段的最小頻率bandMin = i * freqBand# 第i個頻段的最大頻率bandMax = bandMin + freqBand# 判斷第i個頻段是否在要分析的范圍內if (iter_freq['fmin'] <= bandMin and iter_freq['fmax'] >= bandMax):# 給新構造的小波包參數賦值# print('freq',bandMin, bandMax,'fmin',iter_freq['fmin'],'fmax',iter_freq['fmax'])new_wp[freqTree[i]] = wp[freqTree[i]].data# 計算對應頻率的數據bandResult.append(new_wp.reconstruct(update=True))return bandResult########################################小波包變換-重構造分析不同頻段的特征(注意maxlevel，如果太小可能會導致部分頻段分析不到)######################### # 定義WP函數 # epochsData:epochs的數據（mumpy格式） # sfreq:采樣頻率 # wavelet:小波類型 # maxlevel:小波層數 # band:頻帶類型def WP(epochsData, sfreq, wavelet='db4', maxlevel=8, band=bandFreqs):# 輸出的維度順序為 頻率->epoch->channel->timeDataresult = []for epochData in epochsData:channel = []for channelData in epochData:# print('channel:')channel.append(__CalcWP(channelData, sfreq, wavelet=wavelet, maxlevel=maxlevel, band=band))result.append(channel)return np.array(result).transpose((2, 0, 1, 3))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 = 0# 想要分析的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')# 計算小波包濾波后的數據并繪圖wpFilter = WP(epochs.get_data(), sfreq)plt.plot(wpFilter[bandIndex][epochIndex][channelIndex], c=(0, 1, 0), label='WP_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]))# 對小波包濾波后的數據進行FFT變換，需要注意小波包變換后數據的點數可能會發生變化，這里截取數據保持一致性pointNum = epochs.get_data()[epochIndex][channelIndex].shape[0]wpFreq = np.abs(fft.fft(wpFilter[bandIndex][epochIndex][channelIndex][:pointNum]))# 想要繪制的點數pointPlot = 300# FIR濾波后x軸對應的頻率幅值范圍FIR_X = np.linspace(0, sfreq/2, int(mneFIRFreq.shape[0]/2))# 小波包濾波后x軸對應的頻率幅值范圍WP_X = np.linspace(0, sfreq/2, int(wpFreq.shape[0]/2))# 繪制FIR濾波后的頻譜分布plt.plot(FIR_X[:pointPlot], mneFIRFreq[:pointPlot], c=(1, 0, 0), label='FIR_Filter')# 繪制小波包濾波后的頻譜分布plt.plot(WP_X[:pointPlot], wpFreq[:pointPlot], c=(0, 1, 0), label='WP_FIlter')# 繪制圖例和圖名plt.legend()plt.title(bandFreqs[bandIndex]['name'])plt.show()

注：小波包濾波后數據長度會發生改變，這主要和小波層數以及計算方式有關，超出原始長度的數據可以不用關心，上面的代碼中在進行頻譜分析前，也就是計算FFT前對數據的長度進行了處理，這樣可以保證分析出頻譜數據的一致性。

?

WP和FIR在不同頻段濾波效果對比：

從濾波的結果來看，小波包的濾波效果在低頻段還不錯，可是到了高頻段其效果一般。

https://blog.csdn.net/zhoudapeng01/article/details/108214064（FIR和STFT對比）

對比FIR、STFT以及WP這三種方法，STFT相對效果更好些。

?

Delta頻段（1Hz-3Hz）對應的結果

Theta頻段（4Hz-7Hz）對應的結果

Alpha頻段（8Hz-13Hz）對應的結果

Beta頻段（14Hz-31Hz）對應的結果

Gamma頻段（31Hz-40Hz）對應的結果

對應的資源包下載：https://download.csdn.net/download/zhoudapeng01/12793085

總結

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

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