癫痫和正常脑电的AR模型、AR谱及特征比较

课程设计报告
题目：癫痫和正常脑电的AR模型、AR谱及特征比较专业：生物医学工程
班级：09级2班
学号：**********
*名：***
指导教师：田心、刘婷
天津医科大学生物医学工程学院
2012年5 月24 日
（一）开题背景
癫痫是大脑神经元突发性异常放电，导致短暂的大脑功能障碍的一种慢性疾病。

其特征是突然和一过性症状，由于异常放电的神经元在大脑的部位不同，而有多种多样的表现。

2005年国际抗癫痫联盟（ILAE）对癫痫的定义做了修订，其推荐的定义为：癫痫是一种脑部疾患，其特点是持续存在能产生癫痫发作的脑部持久性改变，并出现相应的神经生物学认知心理学以及社会学等方面的后果。

诊断癫痫至少需要一次癫痫发作。

然而，迄今为止，国内学界的共识为反复出现的癫痫发作才可诊断为癫痫，若仅有一次发作就不应诊断为癫痫。

因此，癫痫高发人群正确理解ILAE的癫痫定义，对癫痫的诊断具有非常重要的实际意义。

脑电图是脑神经细胞的电生理活动在大脑皮层或头皮表面表现出的电现象。

临床的常规脑电图检查可为脑损伤、脑血栓、发育异常、新陈代谢和内分泌等疾病的诊断、预后和治疗提供信息。

它也是是脑神经活动的表现，其信息含量丰富。

随着数字处理技术的发展，近年来借助计算机来分析研究脑电图成为一个强有力的手段。

二十年来，在诸如谱分析、波形分析、模型建立等方面都有很大进展。

脑电信号分析[1]已经在脑科学研究中占据了越来越重要的地位，临床医学方面，脑电信号处理可以为某些脑疾病（如癫痫）提供诊断依据和提供有效的治疗手段。

对脑电信号进行预处理，去除噪声和剔除伪差，建立FFT功率谱，AR模型，AR谱等是本次试验的关键。

对脑电信号[1]建立模型的方法有很多种：移动平均（MA）模型、自回归（AR）模型、自回归移动平均（ARMA）模型。

由于AR模型估计的AR谱更适用于脑电信号，因而本文选择的是AR参数模型谱估计[2]。

其中在对AR模型系数的估计方法选择上，我们选用了Marple算法，因为Marple算法使所有的系数都可以自由选择，着眼于调整模型的全部系数使均方误差极小化，能取得更好的结果。

（二）课题目的
1．学习脑电的测试方法
2．在Matlab平台上应用Marple算法建立正常人和癫痫患者脑电的自回归AR模型3．估计AR模型的阶数及各个参数
4．应用模型参数对两类脑电进行频谱估计（AR谱）
（三）课题研究的主要内容
1．测取正常人认知脑区（额部和枕部）的脑电，建立数据文件；对给出的颞叶癫痫患者同一脑区的脑电建立数据文件
2.对以上两类脑电进行预处理，去除噪声和剔除伪差
3.对以上两类脑电建立AR模型
4.估计以上两类脑电的AR谱，并与它们的FFT功率谱进行比较
5.分析比较正常人和癫痫患者在认知脑区脑电的AR模型、AR谱的特征
（四）原理和方法
1. 脑电图及脑电的特点
1.1 脑电图电极安排
为了区分电极和两大脑半球的关系，通常右侧用偶数，左侧用奇数。

从鼻根至枕骨粗隆连一正中矢状线，再从两瞳孔向上、向后与正中矢状线等距的平行线顺延至枕骨粗隆称左右瞳枕线。

1、从枕骨粗隆向上约2cm，左右旁开3cm与左右瞳枕线相交处为左右枕极(9、10)。

2、沿瞳枕线人发际约1cm处为左右额前极(1、 2)。

3、左右外耳道连线与左右瞳枕线相交处为左右中央极(5、6)。

4、左右额前极与中央极之中点处为左右额极(3、4)。

5、左右中央极与枕极之中点处为左右顶极(7、8)。

6、左右中央极与外耳道之中点处为左右颞中极(11、 12)。

7、左右瞳孔与外耳道中点处为左右颞前极(13； 14)。

8、左右乳突上发际内约1em处为左右颞后圾(15、 16)。

9、两耳垂电极的连接方式有两种：
(1)两耳垂连在一起，同时与地线相连，而将无关电极放在头顶中部或鼻根部，或第七颈椎处。

(2)两耳垂分别或共同作为作用电极与头皮电极配对作为脑电团机输入端。

另在头顶处，附加一个电极与地线相连。

目前，国际上已广泛采用10-20系统电极法，其前后方向的测量是以鼻根到枕骨粗隆连成的正中线为准，在此线左右等距的相应部位定出左右前额点(FP1，FP2)、额点(F3，
F4)、中央点(C3．C4)、顶点(P3，，P4)和枕点(O1，O2)，前额点的位置在鼻根上相当于鼻根至枕骨粗隆的10％处，额点在前额点之后相当于鼻根至前额点距离的二倍即鼻根正中线距离20％处，向后中央、顶、枕诸点的间隔均为20％，10-20系统电极的命名即源于此。

1.2 脑电的特点
健康人除个体差异外，在一生不同的年龄阶段，脑电图都各有其特点，但就正常成人脑电图来讲，其波形、波幅、频率和位相等都具有一定的特点。

临床上根据其频率的高低将波形分成以下四种：
β波：频率在13C/S 以上，波幅约为δ波的一半，额部及中央区最明显。

α波：频率在8～ 13C/S ，波幅25～75μV ，以顶枕部最明显，双侧大致同步，重复节律地出现δ波称θ节律。

Φ波：频率为4～7C/S ，波幅20～40μV ，是儿童的正常脑电活动，两侧对称，颞区多见。

δ波：频率为4C/S 以下，δ节律主要在额区，是正常儿童的主要波率，单个的和非局限性的小于20μV 的δ波是正常的，局灶性的δ波则为异常。

δ波和θ波统称为慢波。

在病理状态下如癫痫、脑肿瘤或脑外伤等情况下，大脑活动受到抑制，脑电图出现异常波活动。

癫痫病人的脑电图常见特征波形如下：
棘波：出现方式一般为孤立性或散在性，癫痫病大发作时可以是连续性。

棘波的周期小于 83ms ，其波顶尖锐。

常见于局限性癫痫、癫痫大发作等。

棘慢复合波：棘波和θ波或δ波连接起来的称为棘慢复合波或简称为棘慢波。

其频率有 3Hz 、2Hz 或不规则的。

这种波形是癫痫小发作时特有的。

尖波：又称为锐波，其波形近似于棘波，但周期较长，通常为 83ms 以上，形态为快直上升而缓慢下降的三角波形，波幅可达 200μV 以上，也是一种病理波，多见于癫痫。

尖慢复合波：是由一个尖波和一个慢波组成的复合波，亦称为不全棘慢波或慢棘慢波，会现于局限性癫痫。

2．AR 参数模型
其含义是认为随机信号是有某一激励信号激励一确定性系统所产生的。

只要激励信号的功率和系统的参数已知，随机信号的研究可以转化成对模型参数和性质的研究以及激励信号通过此系统后输入输出的关系。

输入为白噪声的线性系统，输入u （n ）和输出x （n ）的关系可以写作
)()()(0
1
k n u b k n x a n x q
k k p
k k -+--=∑∑== （1）
如果式（1.1）除0b =1以外所有的k b 都等于0，则
)()()(1
n u k n x a n x p
k k +--=∑=
上式为p 阶自回归（AR ）模型。

该模型在n 时刻的输出是n 时刻的输入u （n ）和p 个过去的输出x （n-k ）的加权和。

模型的传递函数为
k
p
k k z a z A z H -=∑+==
111)
(1
)( （2）
3． AR 模型阶次的选择
在AR 模型的参数估计中，模型阶次p 选择是一个关键的问题。

确定AR 模型阶次的方法有很多，较常用的两种是最终预测误差准则（FPE ）和Akaike 信息准则（AIC ）。

（1）最终预测误差准则（FPE ）
)
1()
1()(+-++=k N k N k FPE κ
ρ （3）
（2）Akaike 信息准则（AIC ）
k N k AIC k 2)ln()(+=ρ （4）
式中N 为数据)(n x N 的长度，当阶次k 由1增加时，FPE （k ）和AIC （k ）将在某一个k 处取得极小值，将此时的k 定为最合适的阶次p 。

本文应用的是最终预测误差准则（FPE ）。

4． AR 模型系数的估计—Marple 算法
AR 模型功率谱估计的三种常用算法：L —D 算法、Burg 算法和Marple 算法，：L —D 算法虽然简单，但是分辨率较差；Burg 算法容易出现谱线分裂和伪峰，受初相的影响产生频率偏移，但是它能满足大多数的应用要求，且计算不太复杂；Marple 算法运算量较大，但性能最好，分辨率最高。

Marple 算法又称为不受约束的最小二乘法，它的主要思路是为了摆脱因采用递推运算对确定预测系数的约束，让每一预测系数(模型参数)的确定直接与前、后向预测的总的平方误差最小(最小二乘法)联系起来。

即令总的平方误差最小。

5． AR 谱
即AR 模型估计出的功率谱。

由于AR 模型是一个有理分式，因而估计出的谱要比经典法的谱平滑。

AR 谱分辨率较高，方差反比于数据长度和信噪比。

6．FFT 功率谱
功率谱定义为 2^|)(|1
)(w X N
w p N PER （5）
由于FFT 功率谱是对每一个信号点求得 F 变换，其功率谱是离散的。

经典谱估计都可用 FFT 快速计算。

但谱的分辨率较低（它正比于2pi/N ），而且由于不可避免的窗函数的影响，使得真正的谱P （w ）在窗口主瓣内的功率向边瓣部分“泄露”，降低了分辨率，产生假峰值。

（五）步骤
1. 在计算机上用Matlab 仿真平台显示处理后脑电时序图。

2. 用Matlab 仿真平台建立正常人和癫痫患者脑电的FFT 功率谱，分析比较正常人和癫痫患者脑电的FFT 功率谱。

3. 用Matlab 仿真平台建立正常人和癫痫患者脑电的AR 模型。

4. 用 Matlab 仿真平台估计正常人和癫痫患者脑电的 AR 谱，分析比较正常人和癫痫患者在认知脑区的AR 模型、AR 谱的特征。

（六）结果
6.1在计算机上用Matlab 仿真平台显示处理后脑电时序图
shixutu-normal3-FP1
shixutu-epi1-FP1
shixutu-epi2-FP1
shixutu-epi3-FP1
图1 正常人和癫痫患者的左额区脑电时序图Y/ μV; X/n
shixutu-normal3-01
shixutu-epi1-01
shixutu-epi2-01
shixutu-epi3-01
图2 正常人和癫痫患者的左枕区脑电时序图Y/ μV; X/n
shixutu-normal3-02
shixutu-epi1-02
shixutu-epi2-02
0510152025
shixutu-epi3-02
图3 正常人和癫痫患者的右枕区脑电时序图Y/ μV; X/n
由图1，2，3可以看出，在枕部正常脑电基线较为平稳，相对基线的振幅较小；癫痫脑电基线不平稳，相对基线的振幅较大。

6.2正常人和癫痫患者脑电的FFT 功率谱
4
normal3 fft
4
epil fft
epi2 fft
4
epi3 fft
Frequency (Hz)
图4 正常人和癫痫患者的左额区fft 功率谱
4
normal1 fft C3
4
epil fft C3
4
epi2 fft C3
4
epi3 fft C3
Frequency (Hz)
图5 正常人和癫痫患者的中央区fft 功率谱
4
normal1 fft P4
4
epil fft P4
4
epi2 fft P4
4
epi3 fft P4
Frequency (Hz)
图6 正常人和癫痫患者的左顶点区fft 功率谱
由图4，5，6可以看出正常脑电多为α波，频带范围为10-12Hz ；癫痫脑电多为δ和θ波，频带范围在0-5Hz 。

6.3由AR 模型估计正常人和癫痫患者脑电的AR 谱
0510152025303540
0.511.52
051015
2025303540
0.511.5
24
epil-AR-c3
Frequency (Hz)
图7 正常人和癫痫患者左中央脑电AR 谱
5
10
15
2025
30
35
40
00.511.5
24
normal1-AR-p4
051015
2025303540
5000
10000
epil-AR-p4
Frequency (Hz)
图8 正常人和癫痫患者右顶点脑电AR 谱
0510152025303540
0.511.52x 10
051015
2025303540
5000
10000
epil-AR-f8
Frequency (Hz)
图9 正常人与癫痫患者右颞脑电AR 谱
由图9-16比较可以看出在脑区正常脑电多为α波，频带范围为10-12Hz ；癫痫脑电多为δ和θ波，频带范围在0-5Hz 。

6.4 AR 谱与FFT 功率谱比较
051015
2025303540
01234
4
normal1-fft-p4
Frequency (Hz)A m p l i t u d e
5
10
15
202530
35
40
012
34
4
normal1-AR-p4
Frequency (Hz)
A m p l i t u d e
图12 正常人fft 功率谱与AR 谱比较
由图12可以看出，AR 谱比FFT 功率谱平滑，AR 谱比FFT 功率谱分辨率高
（七）结论
（1）正常脑电基线较为平稳，相对基线的振幅较小；癫痫脑电基线不平稳，相对基线的振幅较大。

（2）正常脑电多为α波，频带范围为10-12Hz ；癫痫脑电多为δ和θ波，频带范围在0-5Hz 。

（3）AR 谱比FFT 功率谱平滑，AR 谱比FFT 功率谱分辨率高
（八）讨论
（1）基于脑电时序图的讨论
正常人额部枕部脑电基线较为平稳，相对于基线的振幅较小，波线较密集；癫痫患者额部枕部脑电基线不平稳，相对于基线的振幅较大，波线较稀疏。

可以推断大脑额部枕部为癫痫患者的发病区。

（2）基于脑电AR 谱的讨论
正常人和癫痫额部枕部脑电的AR 谱及FFT 功率谱有明显差别，正常人额部枕部脑电的功率能量主要集中在 10-12Hz ，即α波段，只有很小一部分能量集中在 0-5Hz ；而癫痫患者额部枕部脑电的功率能量主要集中在0-5Hz ，即δ和θ波段，只有很小一部分能量集中在10-12Hz 。

(3)AR模型阶数的讨论
应用最终预测误差准则（FPE），当阶次k由1增加时，FPE（k）将在某一个k处取得极小值，将此时的k定为最合适的阶次p。

(4)脑电的AR 谱与FFT 功率谱图的对比
FFT 功率谱是每一个信号点求得的傅里叶变换，其功率谱是离散的，从图中我们可以看到它的谱值上下跳动很剧烈，方差性能不好，分辨率低，谱线不平滑。

AR模型估计出的AR谱谱线平滑，分辨率高，方差性能得到改善。

AR谱线近似于FFT功率谱线或者说反映出FFT功率谱的趋势。

一方面说明AR谱性能优于FFT功率谱，另一方面也验证了AR模型的正确性。

（5）癫痫病发病区的讨论
对于除额部枕部外的其他脑区是否为癫痫的发病区，尚未明确。

还须进一步研究。

（九）参考文献
1.李颖洁. 脑电信号分析方法及其应用. 科学出版社. 2009年
2.陈海英. AR模型功率谱估计常用算法的性能比较. 漳州师范学院学报(自然科学
版).2009年第一期(总第63期)
3.生物医学信号处理.天津医科大学生物医学工程系.2009年6月
附录：Matlab程序
第一导联FP1时序图
clc;
clf;
y1=load('D:\normal3.txt');
y2=load('D:\epi1.txt');
y3=load('D:\epi2.txt');
y4=load('D:\epi3.txt');
y1=y1(:,1);
x1=4000/200*(0:3999)/4000; subplot(4,1,1);
plot(x1,y1);
title('shixutu-normal3-FP1'); y2=y2(:,2);
x2=10394/200*(0:10393)/10394; subplot(4,1,2);
plot(x2,y2);
title('shixutu-epi1-FP1');
y3=y3(:,2);
x3=4322/200*(0:4321)/4322; subplot(4,1,3);
plot(x3,y3);
title('shixutu-epi2-FP1’);
y4=y4(:,2);
x4=2782/200*(0:2781)/2782; subplot(4,1,4);
plot(x4,y4);
title('shixutu-epi3-FP1');
第一导联FP1FFT功率谱：
clc;
y1=load('D:\normal3.txt');
y1=y1(:,1);
%滤波器去除工频干扰
n=3;
Wn=[48 52]/100;
[b,a] = butter(n,Wn,'stop'); y1=filter(b,a,y1);
x1=200*(0:1999)/4000;
f1=fft(y1,4000);
f1=abs(f1);
subplot(4,1,1)
plot(x1,f1(1:2000));
axis([0 40 0 5*10^4]);
title('normal3 fft');
y2=load('D:\epi1.txt');
y2=y2(:,2);
x2=200*(0:5196)/10394;
f2=fft(y2);
f2=abs(f2);
subplot(4,1,2)
plot(x2,f2(1:5197));
axis([0 40 0 1.5*10^5]);
title('epil fft');
y3=load('D:\epi2.txt');
y3=y3(:,2);
x3=200*(0:2161)/4322;
f3=fft(y3);
f3=abs(f3);
subplot(4,1,3)
plot(x3,f3(1:2162));
axis([0 40 0 1.5*10^4]);
title('epi2 fft');
y4=load('D:\epi3.txt');
y4=y4(:,2);
x4=200*(0:1391)/2782;
f4=fft(y4);
f4=abs(f4);
subplot(4,1,4)
plot(x4,f4(1:1392));
axis([0 40 0 2*10^4]);
title('epi3 fft');
xlabel('Frequency (Hz)');
第五导联c3的AR模型谱：
clc
p1=90;
y1=load('D:\normal3.txt');
y1=detrend(y1(:,5));
n=3;n=3;
Wn=[48 52]/100;
[b,a] = butter(n,Wn,'stop');
y1=filter(b,a,y1);
fs=200;
[a1,e1] = arburg(y1,p1)
[pxx1,w1] = pmcov(y1,p1); subplot(2,1,1)
plot(w1*fs/2/pi,pxx1) ;
axis([0 40 0 2*10^4]);
title('normal1-AR-c3');
p2=354;
y2=load('D:\epi1.txt');
y2=y2(:,6);
[a2,e2] = armcov(y2,p2)
[pxx2,w2] = pmcov(y2,p2); subplot(2,1,2)
plot(w2*fs/2/pi,pxx2);
axis([0 40 0 2*10^4]);
title('epil-AR-c3');
xlabel('Frequency (Hz)');
正常第八导联（p4）FFT和AR比较：clc;
y1=load('D:\normal3.txt');
y1=y1(:,8);
%滤波器去除工频干扰
n=3;
Wn=[48 52]/100;
[b,a] = butter(n,Wn,'stop');
y1=filter(b,a,y1);
x1=200*(0:1999)/4000;
f1=fft(y1,4000);
f1=abs(f1);
subplot(2,1,1)
plot(x1,f1(1:2000));
title('normal1-fft-p4');
axis([0 40 0 5*10^4]);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
p1=90;
fs=200;
[pxx1,w1] = pmcov(y1,p1); subplot(2,1,2)
plot(w1*fs/2/pi,pxx1) ; axis([0 40 0 5*10^4]); title('normal1-AR-p4'); xlabel('Frequency (Hz)'); ylabel('Amplitude');。