小波去噪和曲线拟合去基线-无力脉6.5.m

matlab仿真小波多分辨率去噪和曲线拟合去噪

clear all;
close all;
levelw=10;
levelc=7;
level=levelw;
points=4096;
subplot(211);
load('C:\MATLAB6p5\work\脉诊\2004.8.4脉象数据\wuli.txt');
x=wuli';
x1=x(11:(points+10))/10000;
plot(x1),grid on,axis tight;


% 小波去噪和基漂
[c,l]=wavedec(x1,level,'dmey');%对f进行10尺度小波分解
ca10=appcoef(c,l,'dmey',10);%提取第10尺度上的低频小波分解系数，即逼近分量
cd10=detcoef(c,l,10);    %提取各尺度上的高频小波分解系数，即细节分量
cd9=detcoef(c,l,9);
cd8=detcoef(c,l,8);    
cd7=detcoef(c,l,7);
cd6=detcoef(c,l,6);   
cd5=detcoef(c,l,5);
cd4=detcoef(c,l,4);
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
caa10=zeros(1,length(ca10));
cdd10=zeros(1,length(cd10));%把第10尺度上的高频小波系数置零
cdd9=zeros(1,length(cd9));
cdd4=zeros(1,length(cd4));%把第4尺度上的高频小波系数置零
cdd3=zeros(1,length(cd3));
cdd2=zeros(1,length(cd2));
cdd1=zeros(1,length(cd1));
c1=[caa10,cd10,cd9,cd8,cd7,cd6,cd5,cd4,cdd3,cdd2,cdd1];%新的小波分解结构
f1=waverec(c1,l,'dmey');%直接去除干扰尺度上对应的小波分量后的重建脉搏信号
figure;
subplot(211);plot(x1);title('原始脉搏信号');grid on;
subplot(212);plot(f1);title('去噪后脉搏信号');grid on;
% 小波去噪和基漂


% 曲线拟合去基漂
%  信号的小波变换，按级给出概貌和细节的波形
level=levelc;
[swa,swd] = swt(f1,level,'dmey');

%求小波变换的模极大值及其位置
ddw=zeros(size(swd));
%pddw=ddw;
nddw=ddw;

negw=swd.*(swd<0);%负数不变，其余置0
ndw=((negw(:,1:points-1)-negw(:,2:points))>0);%大于后面的数的位置置1
nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);%极小值点置1
%ddw=pddw|nddw;
ddw=nddw;
ddw(:,1)=1;
ddw(:,points)=1;
wpeak=ddw.*swd;
wpeak(:,1)=wpeak(:,1)+1e-10;
wpeak(:,points)=wpeak(:,points)+1e-10;
wwpeak=wpeak(level,:);
wwpeak=wwpeak.*(wwpeak<-7);
% wwpeak(4041)=0;

figure;
subplot(level+1,1,1);
plot(f1);grid on;
for i=1:level
    subplot(level+1,1,i+1);
    plot(wpeak(i,:)); axis tight;grid on;
    ylabel(strcat('j=   ',num2str(i)));
end

figure;
subplot(211);
plot(f1);
axis tight;
grid on;
subplot(212);
plot(wwpeak);
axis tight;
grid on;
%哈哈哈哈


peaktmpn=zeros(1,1);
notmpn=zeros(1,1);
count=0;
for i=2:points-60
    if wwpeak(i)<0
        count=count+1;
        notmpn(count)=i;
    end
end
baseline=zeros(1,points);
per=12;
for i=1:per
    baseline(notmpn(i))=f1(notmpn(i));
    peaktmpn(i)=f1(notmpn(i));
end
for j=1:per-1
    for i=notmpn(j)+1:notmpn(j+1)-1
        baseline(i)=(peaktmpn(j+1)-peaktmpn(j))/(notmpn(j+1)-notmpn(j)-2)*(i-notmpn(j)-1)+peaktmpn(j);
    end
end
figure;
subplot(311);plot(f1(notmpn(1):notmpn(per)));grid on;title('含噪脉搏信号');
subplot(312);plot(baseline(notmpn(1):notmpn(per)));grid on;title('基线');
f2=f1(notmpn(1):notmpn(per))-baseline(notmpn(1):notmpn(per));
subplot(313);plot(f2);grid on;title('去基线后脉波信号');
% 曲线拟合去基漂