脉搏波提取程序

clearall;closeall;%%原始信号的导入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-offvaluefc=0.3/fs;alpha=(1-k*cos(2*pi*fc)-sqrt(2*k*(1-cos(2*pi*fc))-k^2*sin(2*pi*fc)^2))/(1-k);y=zeros(size(x));fori=1:size(x,1)y(1,:)=filtfilt(1-alpha,[1-alpha],x(1,:));endx1=x-y;subplot(212);plot(x1);gridon;xlabel('时间(ms)');ylabel('幅值');title('去除基线漂移信号');%%butter带通滤波figure(2);Wp=[0.950]/500;Ws=[0.3140]/500;[n,Wn]=buttord(Wp,Ws,3,10);[b,a]=butter(n,Wn);freqz(b,a,512,fs);grid;title('巴特沃斯带通滤波')grid;%%cheby去除50hz工频干扰figure(3);Wp1=[3065]/500;Ws1=[4555]/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);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([0100]);%%R波的提取A=s21;PM=max(A);MM=min(A);G=(PM-MM)*0.3;%峰值大小波动范围不超过最大波形高度的0.3倍ifmax(A(1:200))==max(A(1:550))[P1(1),t1(1)]=max(A(1:200));cnt1=1;elseifmax(A(1:200))~=max(A(1:550))cnt1=0;endfori=201:length(A)-200%因为周期大于600，波峰大于其左右各200个点内的所有值ifPM-A(i)G&&A(i)==max(A(i-200:i+200))P1(cnt1+1)=A(i);t1(cnt1+1)=i;cnt1=cnt1+1;endendifmax(A(length(A)-199:length(A)))==max(A(length(A)-550:length(A)))[P1(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);gridon;xlabel('峰值数');ylabel('峰值');title('R波序列');%%峰值的特征显示mp=mean(P1);disp('峰值的均值mp=')disp(mp)vp=var(P1);disp('峰值的方差vp=')disp(vp)sp=std(P1);disp('峰值的标准差sp=')disp(sp)%%求周期Tj=1;form=2:length(t1)T(j)=t1(m)-t1(m-1);j=j+1;enddisp('心电信号的周期T=')disp(T)subplot(212);stem(T);gridon;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);xlabel('频率(Hz)');ylabel('功率谱(dB)');title('信号的功率谱图');xlim([050]);[z,f1]=max(p);disp('功率谱最大值z=')disp(z)disp('最大值对应频率f1=')disp(f(f1))