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）