脉搏信号处理y105

数字信号处理综合训练说明书题目：脉搏信号间期序列的时域、频域分析学院：电气工程与信息工程学院班级：电子信息科学与技术（1）班姓名：马艳霞学号：082601372011年7月15日8目录1脉搏信号处理的基本思路…………………………….12信号预处理……………………………………………..23脉搏信号的时域分析…………………………………..54脉搏信号频域分析……………………………………..85程序清单………………………………………………..96心得体会………………………………………………..237参考文献………………………………………………..2495程序清单%加载采样数据；y=[y10y11y12y13y14y15];t=1:6000;%时域范围；figure(1);subplot(2,1,1);y=-y;plot(t,y);gridon;title('脉搏信号');xlabel('时间/ms');ylabel('相对幅度');%截取整数倍周期data1=min(y(1:500));data2=min(y(5000:5500));fori=1:500ify(i)==data1m1=i;endendforj=5000:5500ify(j)==data2m2=j;endend10mb=m2-m1;y=y(m1:m2);t=1:mb+1;subplot(2,1,2);plot(t,y);title('脉搏信号');xlabel('时间/ms');ylabel('相对幅度');gridon;原始的脉搏信号和经过整数倍提取后的信号如下图所示图1%带通滤波器,虑除工频.基线漂移和肌电干扰；figure(2);fs=1500;%采样频率；wp=[0.9,30]/(fs/2);%设置通带截止频率;ws=[0.1,200]/(fs/2);%设置阻带截止频率;rp=0.5;%通带波纹系数;11rs=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('幅度');gridon;y1=filtfilt(num,den,y);%脉搏信号通过零相位滤波器；subplot(2,1,2);plot(y1);title('脉搏信号通过带通滤波器');xlabel('时间/ms');ylabel('幅度');gridon;脉搏信号通过零相位带通滤波器后，高频信号和基线漂移信号有所衰减，信号变得清晰。图212%小波包消噪；%进行一维小波消噪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('幅度');gridon;%小波包消噪，固定阈值nn=length(y1);thr=sqrt(2*log(nn*log(nn)/log(2)));13y2=wpdencmp(y1,'s',4,'db4','sure',thr,1);subplot(2,1,2);plot(y2);title('小波包消噪后的信号y2');xlabel('时间/ms');ylabel('幅度');gridon;图3小波包消除信号中夹杂的突变信号，使得脉搏信号能准确反映人体的生理状况。经过消噪后的脉搏信号如图3所示，比未经处理的信号更清晰，准确。figure(4);s1=fft(y,fs);%原始脉搏信号傅里叶变换；subplot(2,1,1);plot(abs(s1));%幅频特性；axis([0,400,0,40]);%定坐标范围；title('原脉搏信号幅频响应');xlabel('频率/Hz');ylabel('幅度');gridon;14s2=fft(y2,fs);subplot(2,1,2);plot(abs(s2));axis([0,400,0,15]);%定坐标范围；title('消噪后脉搏信号幅频特性');xlabel('频率/Hz');ylabel('幅度');gridon;图4原始的脉搏信号经过带通滤波器以及小波消噪后，主要频率集中在0—40HZ之间，极低频段的基线漂移信号被滤除，高频段的工频干扰和肌电干扰也被滤除。%时域分析；%求1s中出现的最大值，取其0.6倍作为阈值，计算y的平均幅值pRv13=y1;v13=v13-mean(v13);%去基频直流分量;v13=[v13zeros(1,length(y)-length(v13))];15m=ceil(length(v13)/1000)-1;forj=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;forq=1:z[s0(q+1),t(q)]=max(v13(1+500*(q-1):500*q));s8(q)=max(s0(q+1));ifs8(q)=thr1rr=rr+1;%峰值点的个数；endendPrate=20*rr%脉率；%标出y主波的位置；forl=1:length(v13)%范围；ifv13(l)thr1%大于阈值时，v13=v14;v14(l)=v13(l);16elsev14(l)=0;%小于阈值时v14=0;endendz1=round(length(v13)/1000)-1;forj=1:z1[s(j+1),t(j)]=max(v14(1+1000*(j-1):1000*j));s9(j)=max(s(j+1));%[s4,t]=max(x3(2000:2890));forl=1+1000*(j-1):1000*jifv14(l)==s9(j)v14(l)=v14(l);elsev14(l)=0;endendendR1=0;forl=1:length(v14)ifv14(l)0R1=R1+1;%峰值个数；end17endR1=R1-1;figure(5);subplot(2,1,1);plot(v13);%去基频直流分量后的脉搏波;xlabel('时间/ms');ylabel('幅值');gridon;subplot(2,1,2);plot(v14);title('峰值检测点');xlabel('时间/ms');ylabel('幅值');gridon;图5虽然每个脉搏波的波峰值大小并不相同，但他们总在一个范围内波动，波动的范围基本上不超过最大波形高度的0.3倍，为了更为可靠的检测波峰，以0.6倍为参考。因此可以认为波峰点是在每个脉搏周期中波形的最大值附近，大于其邻域内所有点的点。18%提取一个周期波形figure(6);data1=min(y2(750:1500));data2=min(y2(1500:2500));fori=750:1500ify2(i)==data1N11=i;%单周期起点；endendforj=1500:2500ify2(j)==data2N21=j;%单周期终点；endendNc=N21-N11;%单周期长度；y4=y2(N11:N21);subplot(2,1,1);plot(y4);%绘制单周期波形；gridon;title('小波包消噪后的信号y2的单周期波形');xlabel('时间/ms');ylabel('幅值');s4=fft(y4,fs);%傅里叶变换；subplot(2,1,2);19plot(abs(s2));axis([0,400,0,10]);%定坐标范围；title('单周期脉搏信号幅频特性');xlabel('频率/Hz');ylabel('幅度');gridon;图6上图为脉搏信号的单周期时域波形和频域波形图%求R点位置,计算PULSE间期MEAN;R=zeros(1,rr);u3=0;forr=1:(length(v14)-2)v15(r)=v14(r+1)-v14(r);v15(r+1)=v14(r+2)-v14(r+1);ifv15(r)0&&v15(r+1)0u3=u3+1;20R(u3)=r+1;ifu3=2RR(u3-1)=R(u3)-R(u3-1);%间期序列；endendendMEAN=sum(RR)/(R1)%脉搏信号间期平均值；%计算PRR间期总体标准差差SDHR、心率变异指数CV;fori1=1:length(RR)HR(i1)=(RR(i1)-MEAN).^2;endSDRR=sqrt(sum(HR)/(length(HR)-1))%脉搏信号间期标准差；fori1=1:length(RR)-1HR1(i1)=HR(i1+1)-HR(i1);HR2(i1)=HR1(i1).^2;HR3(i1)=mean(HR1);endRMSSD=sqrt(sum(HR2)/(length(HR)-1))%相邻PR间期差的均方根；SDSD=sqrt(sum((HR1-HR3).^2)/(length(HR)-1))%全部PR间期的标准差；CV=SDRR/MEAN%心率变异系数；21根据公式计算的脉搏信号间期序列的平均值，标准差，均方根，差值标准差以及心率的变异系数。%频域分析（功率谱估计）;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('滤波前的功率谱图：PSD—Welch方法');gridon;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');22title('滤波后的功率谱图：PSD—Welch方法');gridon;LF=min(pxx);%低频功率；HF=max(pxx);%高频功率；SE=LF/HF%功率比；图7脉搏信号在频域中主要分析其功率谱，功率谱反映的能量随频率的分布情况，Welch法采用加窗交叠求功率谱，可