脉搏波提取程序

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)

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

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);

disp('峰值的标准差sp=')

disp(sp)

% %求周期T

j=1;

for m=2:length(t1)

T(j)=t1(m)-t1(m-1);

j=j+1;

end

disp ('心电信号的周期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);