心电信号的滤波和特征值提取

clc;clear all;
sig=load（‘ecg101_1.mat‘）;%信号
sig=struct2cell（sig）;
isig=cell2mat（sig）;
t=isig;
t=t‘;

%除低频
[cl]=wavedec（t5‘db5‘）;
d1=wrcoef（‘d‘cl‘db5‘1）;
d2=wrcoef（‘d‘cl‘db5‘2）;
d3=wrcoef（‘d‘cl‘db5‘3）;
d4=wrcoef（‘d‘cl‘db5‘4）;
d5=wrcoef（‘d‘cl‘db5‘5）;
a5=wrcoef（‘a‘cl‘db5‘5）;
x=d1+d2+d3+d4+d5;

%除高频干扰
[cl]=wavedec（x5‘db5‘）;
sigma=wnoisest（cl1）;
alpha=1;
thr=wbmpen（clsigmaalpha）;
t1=wdencmp（‘gbl‘cl‘db5‘5thr‘s‘1）;

%压缩
[cl]=wavedec（t15‘db5‘）;
alpha=1.1;m=l（1）;
[thrnkeep]=wdcbm（clalpham）;
t1=wdencmp（‘lvd‘cl‘db5‘5thr‘h‘）;

%%分层取特征值
[Lo_DHi_D]=orthfilt（‘tbsp‘）;   %低通
Lo_D=（1/sqrt（2））*Lo_D;
Hi_D=（1/sqrt（2））*Hi_D;
N=length（t1）;
Lo_D=[Lo_Dzeros（1N-length（Lo_D））];%补零（圆周卷积，且增大分辨率便于观察）
Hi_D=[Hi_Dzeros（1N-length（Hi_D））];%补零（圆周卷积，且增大分辨率便于观察）
%五层分解
a1=ifft（fft（t1）.*fft（Lo_D））;%一层低频
d1=ifft（fft（t1）.*fft（Hi_D））;%高频
Lo_D1=dyadup（Lo_D‘c‘）;
Hi_D1=dyadup（Hi_D‘c‘）;
Lo_D1=Lo_D1（1:N）;
Hi_D1=Hi_D1（1:N）;
a2=ifft（fft（a1）.*fft（Lo_D1））;%二层低频
d2=ifft（fft（a1）.*fft（Hi_D1））;%高频
Lo_D2=dyadup（Lo_D1‘c‘）;
Hi_D2=dyadup（Hi_D1‘c‘）;
Lo_D2=Lo_D2（1:N）;
Hi_D2=Hi_D2（1:N）;
a3=ifft（fft（a2）.*fft（Lo_D2））;%三层低频
d3=ifft（fft（a2）.*fft（Hi_D2））;%高频
Lo_D3=dyadup（Lo_D2‘c‘）;
Hi_D3=dyadup（Hi_D2‘c‘）;
Lo_D3=Lo_D3（1:N）;
Hi_D3=Hi_D3（1:N）;
a4=ifft（fft（a3）.*fft（Lo_D3））;%四层低频
d4=ifft（fft（a3）.*fft（Hi_D3））;%高频
Lo_D4=dyadup（Lo_D3‘c‘）;
Hi_D4=dyadup（Hi_D3‘c‘）;
Lo_D4=Lo_D4（1:N）;
Hi_D4=Hi_D4（1:N）;
a5=ifft（fft（a4）.*fft（Lo_D4））;%五层低频
d5=ifft（fft（a4）.*fft（Hi_D4））;%高频

%****************正负极大值对***********************
swd=d3;
points=length（d3）;
ddw=zeros（size（swd））;
pddw=ddw;
nddw=ddw;
%小波系数的大于0的点
posw=swd.*（swd>0）;
%斜率大于0
pdw=（（posw（1:points-1）-posw（2:points））<0）;
%正极大值点
pddw（2:points-1）=（（pdw（1:points-2）-pdw（2:points-1））>0）;
%小波系数小于0的点
negw=swd.*（swd<0）;
ndw=（（negw（1:points-1）-negw（2:points））>0）;
%负极大值点
nddw（2:points-1）=（（ndw（1:points-2）-ndw（2:points-1））>0）;
%或运算
ddw=pddw|nddw;
ddw（1）=1;
ddw（points）=1;
%求出极值点的值其他点置0
wpeak=ddw.*swd;
wpeak（1）=wpeak（1）+1e-10;
wpeak（points）=wpeak（points）+1e-10;
%画出各尺度下极值点
figure（6）
plot（wpeak）;axis tight;grid on;
ylabel（strcat（‘j= ‘num2str（i）））
interva2=zeros（1points）;
intervaqs=zeros（1points）;
Mj4=wpeak;
%画出尺度3极值点
figure（7）
plot （Mj4）;
title（‘尺度3下小波系数的模极大值点‘）;
posi=Mj4.*（Mj4>0）;
%求正极大值的平均
thposi=（max（posi（1:round（points/4）））+max（posi（round（points/4）:2*round（points/4）））+max（posi（2*round（points/4）:3*round（points/4）））+max（posi（3*round（points/4）:4*round（points/4））））/4;
posi=（posi>thposi/3）;
nega=Mj4.*（Mj4<0）;
%求负极大值的平均
thnega=（min（nega（1:round（points/4）））+min（nega（round（points/4）:2*round（points/4）））+min（nega（2*round（points/4）:3*round（points/4）））+min（nega（3*round（points/4）:4*round（points/4））））/4;
nega=-1*（nega%找出非0点
interva=posi+nega;
loca=find（interva）;
for i=1:length（loca）-1
if abs（loca（i）-loca（i+1））<80
diff（i）=interva（loca（i））-interva（loca（i+1））;
else
diff（i）=0;
end
end
%找出极值对
loca2=find（diff==-2）