脉搏信号处理y105.doc

数字信号处理综合训练说 明 书题目：脉搏信号间期序列的时域、频域分析学院：电气工程与信息工程学院班级：电子信息科学与技术（1）班姓名：马艳霞学号：08260137 2011年7月15日目录1 脉搏信号处理的基本思路.12 信号预处理.23 脉搏信号的时域分析.54 脉搏信号频域分析.85 程序清单.96 心得体会.237 参考文献.24 5 程序清单%加载采样数据；y=y10 y11 y12 y13 y14 y15;t=1:6000; %时域范围；figure(1);subplot(2,1,1);y=-y; plot(t,y); grid on;title(脉搏信号); xlabel(时间/ms); ylabel(相对幅度);%截取整数倍周期data1=min(y(1:500);data2=min(y(5000:5500);for i=1:500 if y(i)=data1 m1=i; endendfor j=5000:5500 if y(j)=data2 m2=j; endendmb=m2-m1; y=y(m1:m2); t=1:mb+1;subplot(2,1,2); plot(t,y); title(脉搏信号); xlabel(时间/ms); ylabel(相对幅度); grid on;原始的脉搏信号和经过整数倍提取后的信号如下图所示 图1%带通滤波器,虑除工频.基线漂移和肌电干扰；figure(2);fs=1500;%采样频率；wp=0.9,30/(fs/2);%设置通带截止频率;ws=0.1,200/(fs/2);%设置阻带截止频率;rp=0.5;%通带波纹系数;rs=40;%阻带波纹系数;N,wc=buttord(wp,ws,rp,rs);num,den=butter(N,wc);%滤波器分子分母系数向量；H,W=freqz(num,den);subplot(2,1,1);plot(fs*W/(2*pi),abs(H);title(butterworth带通滤波器幅频响应);xlabel(频率/Hz); ylabel(幅度); grid on;y1=filtfilt(num,den,y);%脉搏信号通过零相位滤波器；subplot(2,1,2); plot(y1);title(脉搏信号通过带通滤波器);xlabel(时间/ms); ylabel(幅度); grid on;脉搏信号通过零相位带通滤波器后，高频信号和基线漂移信号有所衰减，信号变得清晰。 图2%小波包消噪；%进行一维小波消噪figure(3);z=wden(y1,heursure,s,mln,3,sym8);%XD,CXD,LXD=wden(X,tptr,sorh,scal,n,wavename)使用小波系数阈值，返回输入信号X除噪后的信号XD，输出参数CXD,LXD表示XD的小波分解结构。 %输入参数中，tptr同thselect（）函数；sorh为 s或h表示软硬阈值；n表示在n层上的小波分解；wavename指定小波名称；scal定义阈值调整比例： %one不设定比例； %sln使用的基于第单层系数噪声估计，设置比例； %mln用噪声层的层相关估计，调整比例。 subplot(2,1,1);plot(z); title(小波消噪后的信号);xlabel(时间/ms); ylabel(幅度); grid on;%小波包消噪，固定阈值nn=length(y1);thr=sqrt(2*log(nn*log(nn)/log(2);y2=wpdencmp(y1,s,4,db4,sure,thr,1);subplot(2,1,2); plot(y2);title(小波包消噪后的信号y2);xlabel(时间/ms); ylabel(幅度); grid on; 图3小波包消除信号中夹杂的突变信号，使得脉搏信号能准确反映人体的生理状况。经过消噪后的脉搏信号如图3所示，比未经处理的信号更清晰，准确。figure(4);s1=fft(y,fs);%原始脉搏信号傅里叶变换；subplot(2,1,1);plot(abs(s1);%幅频特性； axis(0,400,0,40);%定坐标范围；title(原脉搏信号幅频响应);xlabel(频率/Hz); ylabel(幅度); grid on;s2=fft(y2,fs);subplot(2,1,2); plot(abs(s2);axis(0,400,0,15);%定坐标范围；title(消噪后脉搏信号幅频特性);xlabel(频率/Hz); ylabel(幅度); grid on; 图4原始的脉搏信号经过带通滤波器以及小波消噪后，主要频率集中在040HZ之间，极低频段的基线漂移信号被滤除，高频段的工频干扰和肌电干扰也被滤除。%时域分析；%求1s中出现的最大值，取其0.6倍作为阈值，计算y的平均幅值pRv13=y1; v13=v13-mean(v13); %去基频直流分量;v13=v13 zeros(1,length(y)-length(v13);m=ceil(length(v13)/1000)-1;for j=1:m s(j+1),t(j)=max(v13(1+1000*(j-1):1000*j); s5(j)=max(s(j+1); %s4,t=max(x3(2000:2890);endpR=sum(s5)/m;%平均幅值； thr1=0.6*pR;%阈值；%求0.5s中出现的最大值，幅值超过阈值的认为是R波，计算脉率率Prate；rr=0; z=round(length(v13)/500)-1;for q=1:z s0(q+1),t(q)=max(v13(1+500*(q-1):500*q); s8(q)=max(s0(q+1);if s8(q)=thr1 rr=rr+1;%峰值点的个数；endendPrate=20*rr %脉率；%标出y主波的位置；for l=1:length(v13)%范围； if v13(l)thr1 %大于阈值时，v13=v14; v14(l)=v13(l); else v14(l)=0;%小于阈值时v14=0; endendz1=round(length(v13)/1000)-1;for j=1:z1s(j+1),t(j)=max(v14(1+1000*(j-1):1000*j);s9(j)=max(s(j+1);%s4,t=max(x3(2000:2890);for l=1+1000*(j-1):1000*j if v14(l)=s9(j) v14(l)=v14(l); else v14(l)=0; endendendR1=0;for l=1:length(v14) if v14(l)0 R1=R1+1; %峰值个数； endend R1=R1-1;figure(5);subplot(2,1,1);plot(v13); %去基频直流分量后的脉搏波;xlabel(时间/ms); ylabel(幅值); grid on;subplot(2,1,2); plot(v14); title(峰值检测点); xlabel(时间/ms); ylabel(幅值); grid on; 图5虽然每个脉搏波的波峰值大小并不相同，但他们总在一个范围内波动，波动的范围基本上不超过最大波形高度的0.3倍，为了更为可靠的检测波峰，以0.6倍为参考。因此可以认为波峰点是在每个脉搏周期中波形的最大值附近，大于其邻域内所有点的点。%提取一个周期波形figure(6);data1=min(y2(750:1500);data2=min(y2(1500:2500);for i=750:1500 if y2(i)=data1 N11=i;%单周期起点； endendfor j=1500:2500 if y2(j)=data2 N21=j;%单周期终点； endendNc=N21-N11;%单周期长度；y4=y2(N11:N21);subplot(2,1,1);plot(y4);%绘制单周期波形； grid on;title(小波包消噪后的信号y2的单周期波形);xlabel(时间/ms); ylabel(幅值);s4=fft(y4,fs);%傅里叶变换；subplot(2,1,2);plot(abs(s2); axis(0,400,0,10);%定坐标范围；title(单周期脉搏信号幅频特性);xlabel(频率/Hz); ylabel(幅度);grid on; 图6上图为脉搏信号的单周期时域波形和频域波形图%求R点位置,计算PULSE间期MEAN;R=zeros(1,rr);u3=0;for r=1:(length(v14)-2) v15(r)=v14(r+1)-v14(r); v15(r+1)=v14(r+2)-v14(r+1);if v15(r)0&v15(r+1)=2 RR(u3-1)=R(u3)-R(u3-1); %间期序列； endendendMEAN=sum(RR)/(R1) %脉搏信号间期平均值；%计算PRR间期总体标准差差SDHR、心率变异指数CV;for i1=1:length(RR) HR(i1)=(RR(i1)-MEAN).2; end SDRR=sqrt(sum(HR)/(length(HR)-1) %脉搏信号间期标准差； for i1=1:length(RR)-1 HR1(i1)=HR(i1+1)-HR(i1); HR2(i1)=HR1(i1).2; HR3(i1)=mean(HR1); end RMSSD=sqrt(sum(HR2)/(length(HR)-1)%相邻PR间期差的均方根； SDSD=sqrt(sum(HR1-HR3).2)/(length(HR)-1)%全部PR间期的标准差；CV=SDRR/MEAN %心率变异系数；根据公式计算的脉搏信号间期序列的平均值，标准差，均方根，差值标准差以及心率的变异系数。 %频域分析（功率谱估计）;figure(7);y1=y1-mean(y1); %去直流,去趋势；y1=detrend(y1); m1=mb;nfft=256;%设定FFT算法的长度;n1=0:m1-1; t=n1/fs; %数据长度、时间序列window=hanning(256); %选用的窗口noverlap=128; %分段序列重叠的采样点数（长度）dflag=none; %不做趋势处理subplot(2,1,1);pxx=psd(y,nfft,fs,window,noverlap,dflag); %功率谱估计f=(0:nfft/2)*fs/nfft; %求得对应的频率向量plot(f,10*log10(pxx); %绘制功率谱xlabel(频率/Hz); ylabel(功率谱/dB);title(滤波前的功率谱图：PSDWelch方法); grid on;subplot(2,1,2);pxx=psd(y1,nfft,fs,window,noverlap,dflag); %功率谱估计f=(0:nfft/2)*fs/nfft; %求得对应的频率向量plot(f,10*log10(pxx); %绘制功率谱xlabel(频率/Hz);ylabel(功率谱/dB);title(滤波后的功率谱图：PSDWelch方法); grid on;LF=min(pxx);%低频功率；HF=max(pxx);%高频功率；SE=LF/HF %功率比； 图7 脉搏信号在频域中主要分析其功率谱，功率谱反映的能量随频率的分布情况，Welch法采用加窗交叠求功率谱，可以有效减小方差和偏差，一般情况下能接近一致估计的要求，因此采用Welch法。其基本原理是：对数据分段时，使每一段有部分重叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均。对信号加不同的窗函数，谱估计的质量是不同的，相对于矩形窗，Hanning窗的主瓣包括更多的能量，因而使功率谱的的主瓣较窄，分辨率较高域中主要进行功率谱分析，功率谱反应了随机信号各频率成分功率能量分布。因此采用Hanning窗。6 心得体会在本学期的最后两周，学院安排我们进行信号处理综合训练课程设计。在前三周信号检测课设之后学习信号的处理，本次课设内容为脉搏信号的时域和频域分析。利用Matlab软件对数字信号进行处理分析。大约用一周的时间熟悉Matlab,并制定初步方案。利用图书