脉搏波提取程序

1、clear all;close all;%原始信号的导入x=load(D:实验数据libaosenpilaoqianv5.txt);x=x;x=x(2000:25000);fs=1000;figure(1);subplot(211);plot(x);grid;xlabel(时间(ms);ylabel(幅值);title(原始信号);%去除基线漂移k = .7; % cut-off valuefc=0.3/fs;alpha = (1-k*cos(2*pi*fc)-sqrt(2*k*(1-cos(2*pi*fc)-k2*sin(2*pi*fc)2)/(1-k);y = zeros(size(x);

2、for i = 1:size(x,1) y(1,:) = filtfilt(1-alpha,1 -alpha,x(1,:);end x1=x-y;subplot(212);plot(x1);grid on;xlabel(时间(ms);ylabel(幅值);title(去除基线漂移信号);%butter带通滤波figure(2);Wp=0.9 50/500;Ws=0.3 140/500;n,Wn = buttord(Wp,Ws,3,10);b,a = butter(n,Wn);freqz(b,a,512,fs);grid;title(巴特沃斯带通滤波)grid;%cheby去除50hz工频干扰f

3、igure(3);Wp1=30 65/500;Ws1=45 55/500;rp=3;rs=60;n1,Wn1 = cheb2ord(Wp1,Ws1,rp,rs);b1,a1 = cheby2(n1,rs,Wn1,stop);freqz(b1,a1,512,fs);grid;title(带阻切比雪夫型滤波器)grid;%信号之间的对比figure(4);x2=filtfilt(b,a,x1);s21=filtfilt(b1,a1,x2);subplot(313);plot(s21);grid;xlabel(时间（ms）);ylabel(幅值);title(滤波后的信号);subplot(311)

4、;plot(x1);grid;xlabel(时间（ms）);ylabel(幅值);title(原始信号);subplot(312);plot(x2);grid;%频谱图figure(5);F=fft(s21);plot(0:length(s21)-1/(length(s21)*(1/fs),abs(F)*2/length(s21);grid;xlabel(频率/HZ);ylabel(幅值);title(幅频图);xlim(0 100);%R波的提取A=s21;PM=max(A);MM=min(A);G=(PM-MM)*0.3; %峰值大小波动范围不超过最大波形高度的0.3倍 if max(A(

5、1:200)=max(A(1:550)P1(1),t1(1)=max(A(1:200); cnt1=1; elseif max(A(1:200)=max(A(1:550) cnt1=0; endfor i=201:length(A)-200 %因为周期大于600，波峰大于其左右各200个点内的所有值if PM-A(i)G&A(i)=max(A(i-200:i+200) P1(cnt1+1)=A(i); t1(cnt1+1)=i; cnt1=cnt1+1; endendif max(A(length(A)-199:length(A)=max(A(length(A)-550:length(A)P1

6、(cnt1+1),t1(cnt1+1)=max(A(length(A)-199:length(A); cnt1=cnt1+1;endfigure(6);subplot(211);disp(峰值P1)disp(P1);disp(对应时间t1)disp(t1)disp(峰值个数cnt1)disp(cnt1)stem(P1); grid on; xlabel(峰值数);ylabel(峰值);title(R波序列); % % 峰值的特征显示 mp=mean(P1);disp(峰值的均值mp=)disp(mp)vp=var(P1);disp(峰值的方差vp=)disp(vp)sp=std(P1);dis

7、p(峰值的标准差sp=)disp(sp)% %求周期Tj=1;for m=2:length(t1) T(j)=t1(m)-t1(m-1); j=j+1;enddisp (心电信号的周期T=)disp(T)subplot(212);stem(T); grid on; xlabel(心跳数);ylabel(周期值);title(周期序列);mT=mean(T);disp(周期均值mT=)disp(mT)vT=var(T);disp(周期的方差vT=)disp(vT)sT=std(T);disp(周期的标准差sT=)disp(sT)%信号的功率谱估计figure(7);w=hanning(128);nfft=1024;noverlap=64;p,f=spectrum(x,nfft,noverlap,w,fs);plot(f,p);xl