前言
在这里我将例举一种脑信号的功率谱估计,可分为如下几个模块:
- 产生脑信号原始数据
- 对原始脑信号进行FFT变换(查看频谱图)
- 对原始信号的预处理阶段(低通滤波——巴特沃斯滤波器)
- 对滤波信号的小波阈值去噪(软阈值小波去噪)
- 对去噪信号进行FFT变换(查看频谱图)
- 四个阶段的波过带通滤波器绘制
一. 产生脑信号原始数据
通过data = rand(1,9999)产生电信号原始数据,设置采样时间间隔(t=1), 取时间1-3000秒;
原始数据:
对原始信号进行傅里叶变换绘制脑电波信号的频谱图
频谱图:
二. 巴特沃斯滤波器
在对原始信号进行预处理滤波处理,我们这里选用的是巴特沃斯滤波器。
2.1 设置Fs、fp、fs、Ap、As
2.2 计算归一化角频率Wp=fp/(Fs/2)和Ws=fs/(Fs/2);
2.3 计算阶数和截止频率 [N,Wc]=buttord(Wp,Ws,Ap,As)
2.4 计算H(z)分子、分母多项式系数 [b,a]=butter(N,Wc,'low')
2.5 计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率) [H,F]=freqz(b,a,500,Fs)
2.6 绘制频率幅度和频率相位图
2.7 绘制巴特沃斯滤波后的频谱图
滤波器图:
滤波后的时域+频谱图:
三. 小波去噪
这里我们选用软阈值小波去噪,借助wdencmp函数完成
[c,s]=wavedec2(Q,2,'db5'); [thr,sorh,keepapp] = ddencmp('den','wv',Q); [xc,cxc,lxc,perf0,perfl2]=wdencmp('gbl',Q,'sym4',2,thr,sorh,keepapp);
四.功率谱估计
在进行功率谱估计分为四个阶段:
4.1 第一段波过带通滤波器------------delta
4.2 第二段波过带通滤波器------------theta
4.3 第三段波过带通滤波器------------alpha
4.4 第四段波过带通滤波器------------beta
使用核心函数:
f=design(fdesign.bandpass(fs(1),fp(1),fp(2),fs(2),as,rp,as,fo),'butter') 进行功率普估计。
第一段波过带滤波器图:
第二段波过带通滤波器:
第三段波过带通滤波器:
第四段波过带通滤波器: