脉搏波提取程序.docx

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 value fc=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)); 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工频干扰 figure(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); 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(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; end for i=201:length(A)-200 %因为周期大于600，波峰大于其左右各200个点内的所有值 if PM-A(i)GA(i)==max(A(i-200:i+200)) P1(cnt1+1)=A(i); t1(cnt1+1)=i; cnt1=cnt1+1; end end if max(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; end figure(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=) di