清醒与睡意状态脉搏信号分析.doc

数字信号处理综合训练说 明 书题目：清醒与睡意状态脉搏信号分析学院：电气工程与信息工程学院 班级：电子信息科学与技术（2）班 姓名：吴永凯 学号：09260218 指导老师：张爱华 2012年5月31日目录1 摘要.2 2 脉搏信号处理的内容要求.3 3信号预处理.3 4 脉搏信号的时域分析.4 5 脉搏信号频域分析.5 6程序清单.5 7 结果分析.24 8心得体会，致谢.25 9参考文献.26 1 摘要 睡意是伴随着体力疲劳或精神疲劳 而出现的一种状态。睡意会大大降低人的警觉度，决策能力和记忆力。在许多危险性操作和对精度要求比较高的作业环境，如车辆驾驶，航空航天，内外科手术，高空风险作业等，如果操作人员跌入睡意状态，后果将不堪设想。因此，对睡意状态的准确检测和实时估计具有重要的科学意义和迫切的现实需求。 生理信号可以很好的反映出身体各子系统的生理状态和病理变化，广泛地应用于各种临床病症的检测和诊断之中。通过生理记录仪采集脉搏、心电等信号，采用数字滤波、特征提取、分类等方法，对受试者在清醒状态和睡意状态下的脉搏信号进行分析，探讨睡意对脉搏心电信号的影响，以期为睡意检测和评估奠定基础。 2 脉搏信号处理的内容要求 （1）设计滤波器，实现对脉搏信号的噪声抑制。 （2）时域分析：波形特征检测。 （3）功率谱分析：对消噪后的信号进行功率谱分析。要求计算信号的功率谱，峰值频率和多频段能量。 （4）信号特征分析。 3信号预处理人体的脉搏信号是一种微弱的信号，信噪比比较的低，正常人的脉象信号在020Hz的频率范内，且大约99%的能量分布在010Hz，在检测和采集时，由于受到仪器，人体方面的影响，所采集的信号中通常包括干扰信号,如图2.1中的脉搏信号，其中最主要的干扰信号有：（1）基线漂移、人体呼吸等低频干扰，频率小于1Hz；（2）工频干扰，是固定频率的干扰，频率为50Hz；（3）其他一些干扰信号，主要是高频干扰。所以，为了对信号做准确的分析，在分析处理之前必须做一些必要的预处理。针对信号中存在噪声的特点，基线漂移和呼吸等低频干扰在1Hz 以下，而脉搏信号主要在低频范围。在设计滤波器前，介绍一下常用数字滤波器知识：数字滤波器可分为无限冲激响应数字滤波器（IIR）和有限冲激响应数字滤波器（FIR），其具体设计过程不再详述，就一些性能区别作简介。从性能上说，IIR数字滤波器可用较低的阶数获得高的选择性，幅频特性好，但缺点是相位不是严格线性的。相反，FIR滤波器却可以得到严格的线性相位，但FIR只能用较高的阶数达到高的选择性；对于同样的滤波器设计指标，FIR滤波器所要求的阶数可以比IIR滤波器高5到10倍，成本较高，信号延时也较大。综合我们选用IIR滤波器，对常用一些IIR滤波器性能做比较。利用Butterworth滤波器、Chebysheve I型滤波器、Che-bysheve型滤波器、椭圆滤波器都可以进行低通、高通、带通、带阻滤波器的设计，但是各有特点。其中Butterworth滤波器通带内的幅频响应曲线能得到最大限度的平滑，但牺牲了截止频率的坡度；Chebysheve I型滤波器通带内等波纹，阻带内单调；Chebysheve型滤波器通带内单调，然而阻带内等波纹；椭圆滤波器阻带和通带内都是等波纹的，但下降的坡度更大，而且可以以更低的阶数实现和其他两类滤波器一样的性能指标。设计滤波器时要根据实际情况综合考虑截止特性和相位失真的要求等要求。截止特性好的，相位失真就严重，两者不可兼得。Matlab信号工具箱提供直接设计IIR数字滤波器的函数，直接调用这些函数就可以很方便地设计满足我们要求的滤波器。本设计中设计的低通滤波器的参数为：Wp = 30/(fs/2);Ws = 90/(fs/2);Rp = 3; Rs = 60。 选用椭圆低通滤波器滤除其他一些干扰信号等引起的干扰。跟前面一样设计零相位滤波器去噪。 本设计中设计的高通滤波器的参数为Wp1= 0.9/(fs/2);Ws1= 0.1/(fs/2);，Rp=0.5，Rs=40。 50Hz固定工频在频域中是一个点，设计零相位切比雪夫I型滤波器滤除工频干扰。跟前面一样设计零相位滤波器去噪。 本设计中设计的带阻滤波器的参数为：Wp2=45,55/(fs/2);Ws2=49,51/(fs/2);Rp=1;Rs=40; 4 脉搏信号的时域分析脉搏信号的时域特征提取 脉搏信号的特征提取是脉搏信号分析中的关键，能不能准确的提取出来脉搏信号的特征关系到能不能准确的分析脉搏信号，其时域中提取的主要特征有脉搏的峰值、谷值、周期等。而在这其中最重要的是峰值的提取，准确的提取出来峰值，可为进一步分析其周期等其他参数打下基础。 提取信号的峰值基本思想是基于阈值的方法，统计分析采集到得脉搏信号发现，虽然每个脉搏波的波峰值大小并不相同，但他们总在一个范围内波动，波动的范围基本上不超过最大波形高度的0.3倍，为了更为可靠的检测波峰，以0.4倍为参考。因此可以认为波峰点是在每个脉搏周期中波形的最大值附近，大于其邻域内所有点的点。由于所用数据采样频率为1000HZ，通过对大量脉搏波测量数据分析，把邻域长度定义为400个数据点可以准确识别波峰点。 提取出波峰点的位置后，其他时域的特征周期等可相应而出。 峰值提取程序流程如下： 5频域分析脉搏信号频域特征得提取脉搏信号在频域中主要进行功率谱分析，功率谱反应了随机信号各频率成分功率能量分布。为获取清醒与睡意状态脉搏信号能量变化情况，对脉搏信号进行功率谱分析。本设计中用MATLAB中的函数periodogram来实现的饿，其调用形式如下：Pxx,f = periodogram(k,window,n,fs)；其中z表示为滤波后的脉搏信号；window表示为窗函数；n表示为采样点数；fs表示为采样频率；Pxx表示为输出功率谱的估计值，其大小用分贝表示；f表示为对应的频率点。得到功率谱计算最大功率及所对应的频率。具体程序见程序清单6 程序清单清醒状态脉搏信号处理clc;clear;%载入脉搏信号；x=load(E:15guanpeigang_maibo_zao.txt);x=x(1000:19500);N=length(x);x=x+rand(1,N);figure(1);plot(x);title(清醒状态脉搏信号);grid on;%低通滤波器滤波；fs=1000;Wp = 30/(fs/2);Ws = 90/(fs/2);Rp = 3; Rs = 60;n,Wp = ellipord(Wp,Ws,Rp,Rs);b,a = ellip(n,Rp,Rs,Wp);figure(2);freqz(b,a,512,1000); title( Elliptic Lowpass Filter)y=filter(b,a,x);figure(3);plot(y);title(低通滤波后脉搏信号);grid on;%高通滤波器滤波；Wp1= 0.9/(fs/2);Ws1= 0.1/(fs/2); n,Wn=cheb1ord(Wp1,Ws1,0.5,40);b,a = cheby1(n,0.5,Wn,high);figure(4);freqz(b,a,512,1000); title( Chebyshev Type I Highpass Filter)m=filtfilt(b,a,y);figure(5);plot(m);title(高通滤波后脉搏信号);grid on;%带阻滤波器滤波；Wp2=45,55/(fs/2);Ws2=49,51/(fs/2);Rp=1;Rs=40;n,Wn=buttord(Wp2,Ws2,Rp,Rs);b,a=butter(n,Wn,stop);figure(6);freqz(b,a,512,1000); title(Type Butter bandstop Filter)k=filtfilt(b,a,m);figure(7);plot(k);title(带阻滤波后脉搏信号);grid on;%峰值检测;PM=max(k);%取整个信号最大值;figure(8);plot(k);G=(PM-min(k)*0.4;%设置阈值;cnt=0;%设置计数器;for i=1:(length(k)-100) if PM-k(i)G&k(i)=max(k(i-50:i+50) P1(cnt+1)=k(i); T1(cnt+1)=i; cnt=cnt+1; endendstem(T1,P1,r);xlabel(时间（ms);ylabel(幅度);title(峰值检测点);for i=1:length(T1)-1 T3(i)=T1(i+1)-T1(i);endP3=mean(P1); %峰值的平均值;P4=var(P1); %求峰值的方差;Q3=std(P1); %求峰值的标准差;Q4=mean(T3); %求平均周期；%谷值检测;PN=min(k);%取整个信号最小值;figure(9);plot(k);G=(max(k)-PN)*0.4;%设置阈值;cnt=0;%设置计数器;for i=200:(length(k)-200) if k(i)-PNG&k(i)=min(k(i-170:i+170) P2(cnt+1)=k(i); T2(cnt+1)=i; cnt=cnt+1; endendstem(T2,P2,g);xlabel(时间（ms);ylabel(幅度);title(波谷检测点);for i=1:length(T2)-1 T3(i)=T2(i+1)-T2(i);endP5=mean(P2); %谷值的平均值;P6=var(P2); %求谷值的方差;Q5=std(P2); %求谷值的标准差;Q6=mean(T3); %求平均周期；t=1:1:length(k);figure(10);plot(t,k,g,T1,P1,bo,T2,P2,r*);xlabel(时间（ms);ylabel(幅度);title(峰谷检测点);grid on;%FFT变换；s=fft(k,1024);figure(11);f=1:length(s)/2;plot(f,abs(s(1:512);xlabel(频率/(Hz);ylabel(幅度);title(FFT(幅频);grid on;%功率谱；n=length(k);fs=1000;window=hamming(length(k);Pxx,f = periodogram(k,window,n,fs);figure(12);plot(f,Pxx);xlabel(频率(Hz);ylabel(功率谱(dB);title(功率谱);grid on;z2,f2=max(Pxx);%求功率谱最大值时,对应频率值;L=f(f2);峰值的平均值P3 = 4.2554峰值的方差P4= 0.8343峰值的标准差Q3=0.91342平均周期Q4=589.1333谷值的平均值P5=-2.1051谷值的方差P6= 0.23826谷值的标准差Q=50.48812平均周期Q6=769.7667功率谱最大值z2=19.2066功率谱最大值时,对应频率L= 1.0810睡意状态脉搏信号处理clc;clear;%载入脉搏信号；x=load(E:15guanpeigang_maibo_wu.txt);x=x(1000:19500);N=length(x);x=x+rand(1,N);figure(1);plot(x);title(睡意状态脉搏信号);grid on;%低通滤波器滤波；fs=1000;Wp = 30/(fs/2);Ws = 90/(fs/2);Rp = 3; Rs = 60;n,Wp = ellipord(Wp,Ws,Rp,Rs);b,a = ellip(n,Rp,Rs,Wp);figure(2);freqz(b,a,512,1000); title( Elliptic Lowpass Filter)y=filter(b,a,x);figure(3);plot(y);title(低通滤波后脉搏信号);grid on;%高通滤波器滤波；Wp1= 0.9/(fs/2);Ws1= 0.1/(fs/2); n,Wn=cheb1ord(Wp1,Ws1,0.5,40);b,a = cheby1(n,0.5,Wn,high);figure(4);freqz(b,a,512,1000); title( Chebyshev Type I Highpass Filter)m=filtfilt(b,a,y);figure(5);plot(m);title(高通滤波后脉搏信号);grid on;%带阻滤波器滤波；Wp2=45,55/(fs/2);Ws2=49,51/(fs/2);Rp=1;Rs=40;n,Wn=buttord(Wp2,Ws2,Rp,Rs);b,a=butter(n,Wn,stop);figure(6);freqz(b,a,512,1000); title(Type Butter bandstop Filter)k=filtfilt(b,a,m);figure(7);plot(k);title(带阻滤波后脉搏信号);grid on;%峰值检测;PM=max(k);%取整个信号最大值;figure(8);plot(k);G=(PM-min(k)*0.4;%设置阈值;cnt=0;%设置计数器;for i=1:(length(k)-100) if PM-k(i)G&k(i)=max(k(i-50:i+50) P1(cnt+1)=k(i); T1(cnt+1)=i; cnt=cnt+1; endendstem(T1,P1,r);xlabel(时间（ms);ylabel(幅度);title(峰值检测点);for i=1:length(T1)-1 T3(i)=T1(i+1)-T1(i);endP3=mean(P1); %峰值的平均值;P4=var(P1); %求峰值的方差;Q3=std(P1); %求峰值的标准差;Q4=mean(T3); %求平均周期；%谷值检测;PN=min(k);%取整个信号最小值;figure(9);plot(k);G=(max(k)-PN)*0.4;%设置阈值;cnt=0;%设置计数器;for i=200:(length(k)-200) if k(i)-PNG&k(i)=min(k(i-170:i+170) P2(cnt+1)=k(i); T2(cnt+1)=i; cnt=cnt+1; endendstem(T2,P2,g);xlabel(时间（ms);ylabel(幅度);title(波谷检测点);for i=1:length(T2)-1 T3(i)=T2(i+1)-T2(i);endP5=mean(P2); %谷值的平均值;P6=var(P2); %求谷值的方差;Q5=std(P2); %求谷值的标准差;Q6=mean(T3); %求平均周期；t=1:1:length(k);figure(10);plot(t,k,g,T1,P1,bo,T2,P2,r*);xlabel(时间（ms);ylabel(幅度);title(峰谷检测点);grid on;%FFT变换；s=fft(k,1024);figure(11);f=1:length(s)/2;plot(f,abs(s(1:512);xlabel(频率/(Hz);ylabel(幅度);title(FFT(幅频);grid on;%功率谱；n=length(k);fs=1000;window=hamming(length(k);Pxx,f = periodogram(k,window,n,fs);figure(12);plot(f,Pxx);xlabel(频率(Hz);ylabel(功率谱(dB);title(功率谱);grid on;z2,f2=max(Pxx);%求功率谱最大值时,对应频率值;L=f(f2);峰值的平均值P3= 6.5493峰值的方差P4= 2.9414峰值的标准差Q3= 1.7150平均周期Q4=973.7222谷值的平均值P5= -2.3926谷值的方差P6=0.29069谷值的标准差Q5= 0.53916平均周期Q6= 993.5000功率谱最大值z2=34.7902功率谱最大值时,对应频率L= 1.0270 7结果分析提取出脉搏的特征信号，通过分析这些不同的特征信息，找出清醒与睡意变化比较显著的信息，作为清醒与睡意的测评依据。提取出来清醒与睡意的基本特征参数如下表所示：睡意状态峰值平均值(v)峰值的方差峰值的标准差平均周期(ms)4.25540.83430.91342589.1333谷值平均值(v)谷值的方差谷值的标准差平均周期(ms)-2.10510.2382650.48812769.7667功率谱最大值（dB)功率谱最大值对应频率(Hz)19.20661.0810清醒状态峰值平均值(v)峰值的方差峰值的标准差平均周期(ms)6.54932.94141.7150973.7222谷值平均值(v)谷值的方差谷值的标准差平均周期(ms)-2.39260.290690.53916993.5000功率谱最大值（dB)功率谱最大值对应频率(Hz)34.79021.0270对比清醒与睡意的基本参数，可以看出睡意状态的峰值平均值、峰值的标准差、最大功率和最大功率所对应的频率比清醒状态的明显变大，而平均周期却明显变小，故我们以期用这些参数作为清