心电特征提取程序

clc;
tic
x=load （‘E:\心电脉搏实验资料\20110412睡意状态\xuwei_v5_0.txt‘）;
x=x‘;  %转为行向量
x=x（660001:720000）;
d=x;
sfreq=1000;%采样 
s=length（d）;%信号长度
s2=1./sfreq;%采样周期
s3=round（s./s2）;%有多少周期信号，若是小数区里最近的整数
N=length（d）;
r=1:1:N-1;
C=[1r+1];% 序列 x对应的采样点数
fs=1000;%采样频率
figure（1）
subplot（211）plot（d）;grid;
xlabel（‘时间/ms‘）;
ylabel（‘幅值/mv‘）;
title（‘原始信号‘）;
axis tight;
%低通去除基线漂移
k = .7; % cut-off value
fc=0.8/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（d））;
for i = 1:size（d1）  %可以去掉
   y（i:） = filtfilt（1-alpha[1 -alpha]d（i:））;
end
x=d-y;
subplot（212）;
plot（x）;grid;
xlabel（‘时间/ms‘）;
ylabel（‘幅值/mv‘）;
title（‘去除基线漂移后的信号‘）;
axis tight;
% bandpass滤波
% figure（3）;
[ba]=butter（2[0.5 100]/200‘bandpass‘）;
y=filter（bax）;
% 椭圆带阻滤波器
figure（4）;
wp1=48;
wp2=52;
ws1=49;
ws2=51;
wp=[wp1 wp2];
ws=[ws1 ws2];
rp=0.1;
rs=50;
fs=1000;
[nwn]=ellipord（2*wp/fs2*ws/fsrprs）;
[numden]=ellip（nrprswn‘stop‘）;
[hw]=freqz（numden）;
% 椭圆带阻滤波器滤除50Hz工频干扰
h=filter（numdeny）;
%利用小波分解进行第一次降噪
[cl]=wavedec（h5‘coif5‘）;  %小波分解
a5=wrcoef（‘a‘cl‘coif5‘5）;%重构第五层逼近信号
d5=wrcoef（‘d‘cl‘coif5‘5）;%重构第1--5层细节信号
d4=wrcoef（‘d‘cl‘coif5‘4）;
d3=wrcoef（‘d‘cl‘coif5‘3）;
d2=wrcoef（‘d‘cl‘coif5‘2）;
d1=wrcoef（‘d‘cl‘coif5‘1）;
a2=wrcoef（‘a‘cl‘coif5‘2）;
%去除躁声，重构信号
y=a5+d5+d4;
%进行二次降噪
wavelet=‘sym5‘;
level=5;
alpha=1.5;
sorh=‘s‘;
[cl]=wavedec（ylevelwavelet）;
[thrnkeep]=wdcbm（clalpha）;
[xdcxdlxdperf0perfl2]=wdencmp（‘lvd‘ clwaveletlevelthrsorh）;%重构信号
xx=xd;
for n=13:length（xx）-11
   ex（n） = （xx（n）^2） - （（xx（n-10）） * （xx（n+10）））;%采用Teager能量算子的方法来确定ECG信号R波的位置
end
subplot（311）plot（ex‘r‘）;grid;
axis tight
axis xy
%*************************************************
m=bartlett（7）;
mm=m（2:6）;
eey=conv（exmm）;
subplot（312）;plot（eey‘b‘）;grid;
axis tight
for i=1:length（eey）-1
%     %eey1（i）=eey（i+1）-eey（i）;
%     %eeym=max（eey）-min（eey）;
%     %eey2（i）=abs（eey1（i）/eeym）;
    eey1（i）=abs（eey（i））;
end
subplot（313）;
plot（eey1‘g‘）;grid;
axis tight
c=10;
sm0=0;
for i=1:length（eey1）
    sm0=sm0+eey1（i）;
end
sm1=1/length（eey1）;
thr2=c*sm1*sm0;
%%%%%%
eey2=eey1;
for i=1:length（eey2）
   if eey2（i）>=thr2 
      reg1（i）=1;
   else
      reg1（i）=0;
   end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k3=0;
k4=0;
k1=round（0.5*sfreq）;
for n=1:round（length（eey2）/k1）
  for i1=（（n-1）*k1+1）:n*k1
    for k2=1:10
    if i1<=（length（reg1）-11） & reg1（i1）==1 & reg1（i1+1）==0 & reg1（i1+1+k2）==1
       reg1（:[i1:i1+1+k2]）=1;
    end
    end
  end
end
for n=1:round（length（eey2）/k1）
  for i2=（（n-1）*k1+1）:n*k1
      if i2<=length（reg1）-1 & reg1（i2）==0 & reg1（i2+1）==1
         k3=k3+1;
         qrs1（k3）=i2;
      end
      if i2<=（length（reg1）） & reg1（i2）==1 & reg1（i2+1）==0
        k4=k4+1;
        qrs2（k4）=i2;
      end
  end
end
c=zeros（1k3）;b=zeros（1k3）;mr=[];
 for m1=1:k3 
       [c（m1）b（m1）]=max（xd（qrs1（m1）:qrs2（m1）