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

博客對應(yīng)的代碼和數(shù)據(jù)

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 # 設(shè)置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)# 頻譜由低到高的對應(yīng)關(guān)系，這里需要注意小波變換的頻帶排列默認(rèn)并不是順序排列，所以這里需要使用’freq‘排序。freqTree = [node.path for node in wp.get_level(maxlevel, 'freq')]# 計算maxlevel最小頻段的帶寬，采樣頻率的一半freqBand = (sfreq/2) / (2 ** maxlevel)bandResult = []#######################根據(jù)實際情況計算頻譜對應(yīng)關(guān)系，這里要注意系數(shù)的順序for iter_freq in band:# 構(gòu)造空的小波包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個頻段是否在要分析的范圍內(nèi)if (iter_freq['fmin'] <= bandMin and iter_freq['fmax'] >= bandMax):# 給新構(gòu)造的小波包參數(shù)賦值# print('freq',bandMin, bandMax,'fmin',iter_freq['fmin'],'fmax',iter_freq['fmax'])new_wp[freqTree[i]] = wp[freqTree[i]].data# 計算對應(yīng)頻率的數(shù)據(jù)bandResult.append(new_wp.reconstruct(update=True))return bandResult########################################小波包變換-重構(gòu)造分析不同頻段的特征(注意maxlevel，如果太小可能會導(dǎo)致部分頻段分析不到)######################### # 定義WP函數(shù) # epochsData:epochs的數(shù)據(jù)（mumpy格式） # sfreq:采樣頻率 # wavelet:小波類型 # maxlevel:小波層數(shù) # 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格式的數(shù)據(jù)epochs = mne.read_epochs(r'F:\BaiduNetdiskDownload\BCICompetition\BCICIV_2a_gdf\Train\Fif\A02T_epo.fif')# 繪圖驗證結(jié)果plt.figure(figsize=(15, 10))# 獲取采樣頻率sfreq = epochs.info['sfreq']# 想要分析的目標(biāo)頻帶bandIndex = 0# 想要分析的channelchannelIndex = 0# 想要分析的epochepochIndex = 0# 繪制原始數(shù)據(jù)plt.plot(epochs.get_data()[epochIndex][channelIndex], label='Raw')# 計算FIR濾波后的數(shù)據(jù)并繪圖（注意這里要使用copy方法，否則會改變原始數(shù)據(jù)）firFilter = epochs.copy().filter(bandFreqs[bandIndex]['fmin'], bandFreqs[bandIndex]['fmax'])plt.plot(firFilter.get_data()[epochIndex][channelIndex], c=(1, 0, 0), label='FIR_Filter')# 計算小波包濾波后的數(shù)據(jù)并繪圖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濾波后的數(shù)據(jù)進(jìn)行FFT變換mneFIRFreq = np.abs(fft.fft(firFilter.get_data()[epochIndex][channelIndex]))# 對小波包濾波后的數(shù)據(jù)進(jìn)行FFT變換，需要注意小波包變換后數(shù)據(jù)的點數(shù)可能會發(fā)生變化，這里截取數(shù)據(jù)保持一致性pointNum = epochs.get_data()[epochIndex][channelIndex].shape[0]wpFreq = np.abs(fft.fft(wpFilter[bandIndex][epochIndex][channelIndex][:pointNum]))# 想要繪制的點數(shù)pointPlot = 300# FIR濾波后x軸對應(yīng)的頻率幅值范圍FIR_X = np.linspace(0, sfreq/2, int(mneFIRFreq.shape[0]/2))# 小波包濾波后x軸對應(yīng)的頻率幅值范圍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()

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

?

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

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

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

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

?

Delta頻段（1Hz-3Hz）對應(yīng)的結(jié)果

Theta頻段（4Hz-7Hz）對應(yīng)的結(jié)果

Alpha頻段（8Hz-13Hz）對應(yīng)的結(jié)果

Beta頻段（14Hz-31Hz）對應(yīng)的結(jié)果

Gamma頻段（31Hz-40Hz）對應(yīng)的結(jié)果

對應(yīng)的資源包下載：https://download.csdn.net/download/zhoudapeng01/12793085

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。