小波作为一种信号处理的工具在脑波分析中应用很多,常用的有连续小波变换、小波包分析等等。小波涉及的相关介绍和公式推导有很多资料,文章末尾推荐了几个链接。本文主要介绍连续小波变换,小波包分解重构,对应频段能量计算这3种应用在Python中的实现。
数据来源为BCI竞赛公开数据集中的部分数据,剔除了无效数据。有关数据的描述见链接:
1、连续小波变换(主要用于时频域分析)
这里使用连续小波变换进行时频域分析,数据只是示例,代码中的参数在实际应用的时候需要根据实际情况进行调整。代码中有关小波尺度的计算很有意思,这里单独拿出来详细说明下。
一般用小波的尺度来衡量小波的频率,两者之间的转换关系为:
其中Fs为采样频率,wcf为小波的中心频率。
我们以下文代码中的参数为例,当一共想要分析totalscal个频率的时候,第i个频率
对应的是
,带入上面的等式中
该部分对应的代码如下:
import numpy as np
import matplotlib.pyplot as plt
import pywt
import mne
mne.set_log_level(False)
######################################################连续小波变换##########
# totalscal小波的尺度,对应频谱分析结果也就是分析几个(totalscal-1)频谱
def TimeFrequencyCWT(data,fs,totalscal,wavelet='cgau8'):
# 采样数据的时间维度
t = np.arange(data.shape[0])/fs
# 中心频率
wcf = pywt.central_frequency(wavelet=wavelet)
# 计算对应频率的小波尺度
cparam = 2 * wcf * totalscal
scales = cparam/np.arange(totalscal, 1, -1)
# 连续小波变换
[cwtmatr, frequencies] = pywt.cwt(data, scales, wavelet, 1.0/fs)
# 绘图
plt.figure(figsize=(8, 4))
plt.subplot(211)
plt.plot(t, data)
plt.xlabel(u"time(s)")
plt.title(u"Time spectrum")