denoise.m

基于db8小波的

function out=denoise(in)
%基于db8小波的，平移不变小波，层依赖硬阈值去噪。
%分解层数为5层

db8_g=[0.230377813;0.714846571;0.630880768;-0.027983769;-0.187034812;0.030841382;0.032883012;-0.010597402];%尺度滤波器
db8_h=[-0.010597402;-0.032883012;0.030841382;0.187034812;-0.027983769;-0.630880768;0.714846571;-0.230377813];%小波滤波器

Ldata=3979;%待处理的数据长度
Lfilter=8;%滤波器的长度

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%数组声明


%in=zeros(Ldata,1);%声明待处理的数据
out=zeros(Ldata,1);%声明处理结果

filter_g1=zeros(8,1);
filter_h1=zeros(8,1);
filter_g2=zeros(15,1);
filter_h2=zeros(15,1);
filter_g3=zeros(29,1);
filter_h3=zeros(29,1);
filter_g4=zeros(57,1);
filter_h4=zeros(57,1);
filter_g5=zeros(113,1);
filter_h5=zeros(113,1);

W1=zeros(Ldata,1);
V1=zeros(Ldata,1);
W2=zeros(Ldata,1);
V2=zeros(Ldata,1);
W3=zeros(Ldata,1);
V3=zeros(Ldata,1);
W4=zeros(Ldata,1);
V4=zeros(Ldata,1);
W5=zeros(Ldata,1);
V5=zeros(Ldata,1);

W1_dn=zeros(Ldata,1);
W2_dn=zeros(Ldata,1);
W3_dn=zeros(Ldata,1);
W4_dn=zeros(Ldata,1);
W5_dn=zeros(Ldata,1);

h1=zeros(8,1);
g1=zeros(8,1);
h2=zeros(8+15-1,1);
g2=zeros(8+15-1,1);
h3=zeros(8+15-1+29-1,1);
g3=zeros(8+15-1+29-1,1);
h4=zeros(8+15-1+29-1+57-1,1);
g4=zeros(8+15-1+29-1+57-1,1);
h5=zeros(8+15-1+29-1+57-1+113-1,1);
g5=zeros(8+15-1+29-1+57-1+113-1,1);

D1=zeros(Ldata,1);
D2=zeros(Ldata,1);
D3=zeros(Ldata,1);
D4=zeros(Ldata,1);
D5=zeros(Ldata,1);
S5=zeros(Ldata,1);

D1_dn=zeros(Ldata,1);
D2_dn=zeros(Ldata,1);
D3_dn=zeros(Ldata,1);
D4_dn=zeros(Ldata,1);
D5_dn=zeros(Ldata,1);

%分解时的中间变量
temp_h11=zeros(length(filter_h1)-1,1);%声明中间过程数据
temp_h12=zeros(Ldata+length(filter_h1)-1,1);%声明中间过程数据
temp_h13=zeros(Ldata+length(filter_h1)-1+length(filter_h1)-1,1);%声明中间过程数据
temp_g11=zeros(length(filter_g1)-1,1);%声明中间过程数据
temp_g12=zeros(Ldata+length(filter_g1)-1,1);%声明中间过程数据
temp_g13=zeros(Ldata+length(filter_g1)-1+length(filter_g1)-1,1);%声明中间过程数据

temp_h21=zeros(length(filter_h2)-1,1);%声明中间过程数据
temp_h22=zeros(Ldata+length(filter_h2)-1,1);%声明中间过程数据
temp_h23=zeros(Ldata+length(filter_h2)-1+length(filter_h2)-1,1);%声明中间过程数据
temp_g21=zeros(length(filter_g2)-1,1);%声明中间过程数据
temp_g22=zeros(Ldata+length(filter_g2)-1,1);%声明中间过程数据
temp_g23=zeros(Ldata+length(filter_g2)-1+length(filter_g2)-1,1);%声明中间过程数据

temp_h31=zeros(length(filter_h3)-1,1);%声明中间过程数据
temp_h32=zeros(Ldata+length(filter_h3)-1,1);%声明中间过程数据
temp_h33=zeros(Ldata+length(filter_h3)-1+length(filter_h3)-1,1);%声明中间过程数据
temp_g31=zeros(length(filter_g3)-1,1);%声明中间过程数据
temp_g32=zeros(Ldata+length(filter_g3)-1,1);%声明中间过程数据
temp_g33=zeros(Ldata+length(filter_g3)-1+length(filter_g3)-1,1);%声明中间过程数据

temp_h41=zeros(length(filter_h4)-1,1);%声明中间过程数据
temp_h42=zeros(Ldata+length(filter_h4)-1,1);%声明中间过程数据
temp_h43=zeros(Ldata+length(filter_h4)-1+length(filter_h4)-1,1);%声明中间过程数据
temp_g41=zeros(length(filter_g4)-1,1);%声明中间过程数据
temp_g42=zeros(Ldata+length(filter_g4)-1,1);%声明中间过程数据
temp_g43=zeros(Ldata+length(filter_g4)-1+length(filter_g4)-1,1);%声明中间过程数据

temp_h51=zeros(length(filter_h5)-1,1);%声明中间过程数据
temp_h52=zeros(Ldata+length(filter_h5)-1,1);%声明中间过程数据
temp_h53=zeros(Ldata+length(filter_h5)-1+length(filter_h5)-1,1);%声明中间过程数据
temp_g51=zeros(length(filter_g5)-1,1);%声明中间过程数据
temp_g52=zeros(Ldata+length(filter_g5)-1,1);%声明中间过程数据
temp_g53=zeros(Ldata+length(filter_g5)-1+length(filter_g5)-1,1);%声明中间过程数据

%重构时的中间变量
temp_H11=zeros(length(h1)-1,1);%声明中间过程数据
temp_H12=zeros(Ldata+length(h1)-1,1);%声明中间过程数据
temp_H13=zeros(Ldata+length(h1)-1+Ldata+length(h1)-1-1,1);%声明中间过程数据
temp_G11=zeros(length(g1)-1,1);%声明中间过程数据
temp_G12=zeros(Ldata+length(g1)-1,1);%声明中间过程数据
temp_G13=zeros(Ldata+length(g1)-1+Ldata+length(g1)-1-1,1);%声明中间过程数据

temp_H21=zeros(length(h2)-1,1);%声明中间过程数据
temp_H22=zeros(Ldata+length(h2)-1,1);%声明中间过程数据
temp_H23=zeros(Ldata+length(h2)-1+Ldata+length(h2)-1-1,1);%声明中间过程数据
temp_G21=zeros(length(g2)-1,1);%声明中间过程数据
temp_G22=zeros(Ldata+length(g2)-1,1);%声明中间过程数据
temp_G23=zeros(Ldata+length(g2)-1+Ldata+length(g2)-1-1,1);%声明中间过程数据

temp_H31=zeros(length(h3)-1,1);%声明中间过程数据
temp_H32=zeros(Ldata+length(h3)-1,1);%声明中间过程数据
temp_H33=zeros(Ldata+length(h3)-1+Ldata+length(h3)-1-1,1);%声明中间过程数据
temp_G31=zeros(length(g3)-1,1);%声明中间过程数据
temp_G32=zeros(Ldata+length(g3)-1,1);%声明中间过程数据
temp_G33=zeros(Ldata+length(g3)-1+Ldata+length(g3)-1-1,1);%声明中间过程数据

temp_H41=zeros(length(h4)-1,1);%声明中间过程数据
temp_H42=zeros(Ldata+length(h4)-1,1);%声明中间过程数据
temp_H43=zeros(Ldata+length(h4)-1+Ldata+length(h4)-1-1,1);%声明中间过程数据
temp_G41=zeros(length(g4)-1,1);%声明中间过程数据
temp_G42=zeros(Ldata+length(g4)-1,1);%声明中间过程数据
temp_G43=zeros(Ldata+length(g4)-1+Ldata+length(g4)-1-1,1);%声明中间过程数据

temp_H51=zeros(length(h5)-1,1);%声明中间过程数据
temp_H52=zeros(Ldata+length(h5)-1,1);%声明中间过程数据
temp_H53=zeros(Ldata+length(h5)-1+Ldata+length(h5)-1-1,1);%声明中间过程数据
temp_G51=zeros(length(g5)-1,1);%声明中间过程数据
temp_G52=zeros(Ldata+length(g5)-1,1);%声明中间过程数据
temp_G53=zeros(Ldata+length(g5)-1+Ldata+length(g5)-1-1,1);%声明中间过程数据

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%求W1
filter_h1=db8_h;
for jj=1:length(filter_h1)-1
    temp_h11(8-jj)=in(Ldata+1-jj);
end
temp_h12=[temp_h11;in];
temp_h13=conv(filter_h1,temp_h12);
for kk=1:Ldata
    W1(kk)=temp_h13(kk+7);
end

%求V1
filter_g1=db8_g;
for jj=1:length(filter_g1)-1
    temp_g11(8-jj)=in(Ldata+1-jj);
end
temp_g12=[temp_g11;in];
temp_g13=conv(filter_g1,temp_g12);
for kk=1:Ldata
    V1(kk)=temp_g13(kk+7);
end

%求W2
for ii=1:length(filter_h1)
    filter_h2(2*ii-1)=filter_h1(ii);
end
for jj=1:length(filter_h2)-1
    temp_h21(15-jj)=V1(Ldata+1-jj);
end
temp_h22=[temp_h21;V1];
temp_h23=conv(filter_h2,temp_h22);
for kk=1:Ldata
    W2(kk)=temp_h23(kk+14);
end

%求V2
for ii=1:length(filter_g1)
    filter_g2(2*ii-1)=filter_g1(ii);
end
for jj=1:length(filter_g2)-1
    temp_g21(15-jj)=V1(Ldata+1-jj);
end
temp_g22=[temp_g21;V1];
temp_g23=conv(filter_g2,temp_g22);
for kk=1:Ldata
    V2(kk)=temp_g23(kk+14);
end

%求W3
for ii=1:length(filter_h2)
    filter_h3(2*ii-1)=filter_h2(ii);
end
for jj=1:length(filter_h3)-1
    temp_h31(29-jj)=V2(Ldata+1-jj);
end
temp_h32=[temp_h31;V2];
temp_h33=conv(filter_h3,temp_h32);
for kk=1:Ldata
    W3(kk)=temp_h33(kk+28);
end

%求V3
for ii=1:length(filter_g2)
    filter_g3(2*ii-1)=filter_g2(ii);
end
for jj=1:length(filter_g3)-1
    temp_g31(29-jj)=V2(Ldata+1-jj);
end
temp_g32=[temp_g31;V2];
temp_g33=conv(filter_g3,temp_g32);
for kk=1:Ldata
    V3(kk)=temp_g33(kk+28);
end

%求W4
for ii=1:length(filter_h3)
    filter_h4(2*ii-1)=filter_h3(ii);
end
for jj=1:length(filter_h4)-1
    temp_h41(57-jj)=V3(Ldata+1-jj);
end
temp_h42=[temp_h41;V3];
temp_h43=conv(filter_h4,temp_h42);
for kk=1:Ldata
    W4(kk)=temp_h43(kk+56);
end

%求V4
for ii=1:length(filter_g3)
    filter_g4(2*ii-1)=filter_g3(ii);
end
for jj=1:length(filter_g4)-1
    temp_g41(57-jj)=V3(Ldata+1-jj);
end
temp_g42=[temp_g41;V3];
temp_g43=conv(filter_g4,temp_g42);
for kk=1:Ldata
    V4(kk)=temp_g43(kk+56);
end

%求W5
for ii=1:length(filter_h4)
    filter_h5(2*ii-1)=filter_h4(ii);
end
for jj=1:length(filter_h5)-1
    temp_h51(113-jj)=V4(Ldata+1-jj);
end
temp_h52=[temp_h51;V4];
temp_h53=conv(filter_h5,temp_h52);
for kk=1:Ldata
    W5(kk)=temp_h53(kk+112);
end

%求V5
for ii=1:length(filter_g4)
    filter_g5(2*ii-1)=filter_g4(ii);
end
for jj=1:length(filter_g5)-1
    temp_g51(113-jj)=V4(Ldata+1-jj);
end
temp_g52=[temp_g51;V4];
temp_g53=conv(filter_g5,temp_g52);
for kk=1:Ldata
    V5(kk)=temp_g53(kk+112);
end

%求W1_dn
noise1=median(abs(W1))/0.6745;%经验公式
threshold1=noise1*((2*log(3979))^0.5);
for kk=1:Ldata
    if abs(W1(kk))>threshold1
        W1_dn(kk)=W1(kk);
    end
end

%求W2_dn
noise2=median(abs(W2))/0.6745;%经验公式
threshold2=noise2*((2*log(3979))^0.5);
for kk=1:Ldata
    if abs(W2(kk))>threshold2
        W2_dn(kk)=W2(kk);
    end
end

%求W3_dn
noise3=median(abs(W3))/0.6745;%经验公式
threshold3=noise3*((2*log(3979))^0.5);
for kk=1:Ldata
    if abs(W3(kk))>threshold3
        W3_dn(kk)=W3(kk);
    end
end

%求W4_dn
noise4=median(abs(W4))/0.6745;%经验公式
threshold4=noise4*((2*log(3979))^0.5);
for kk=1:Ldata
    if abs(W4(kk))>threshold4
        W4_dn(kk)=W4(kk);
    end
end

%求W5_dn
noise5=median(abs(W5))/0.6745;%经验公式
threshold5=noise5*((2*log(3979))^0.5);
for kk=1:Ldata
    if abs(W5(kk))>threshold5
        W5_dn(kk)=W5(kk);
    end
end

%求D1
h1=filter_h1;
for jj=1:length(h1)-1
    temp_H11(jj)=W1(jj);
end
temp_H12=[W1;temp_H11];
temp_H13=xcorr(temp_H12,h1);
for kk=1:Ldata
    D1(kk)=temp_H13(kk+Ldata+length(h1)-1-1);
end

%求D2
h2=conv(filter_g1,filter_h2);
for jj=1:length(h2)-1
    temp_H21(jj)=W2(jj);
end
temp_H22=[W2;temp_H21];
temp_H23=xcorr(temp_H22,h2);
for kk=1:Ldata
    D2(kk)=temp_H23(kk+Ldata+length(h2)-1-1);
end

%求D3
h3=conv(conv(filter_g1,filter_g2),filter_h3);
for jj=1:length(h3)-1
    temp_H31(jj)=W3(jj);
end
temp_H32=[W3;temp_H31];
temp_H33=xcorr(temp_H32,h3);
for kk=1:Ldata
    D3(kk)=temp_H33(kk+Ldata+length(h3)-1-1);
end

%求D4
h4=conv(conv(conv(filter_g1,filter_g2),filter_g3),filter_h4);
for jj=1:length(h4)-1
    temp_H41(jj)=W4(jj);
end
temp_H42=[W4;temp_H41];
temp_H43=xcorr(temp_H42,h4);
for kk=1:Ldata
    D4(kk)=temp_H43(kk+Ldata+length(h4)-1-1);
end

%求D5
h5=conv(conv(conv(conv(filter_g1,filter_g2),filter_g3),filter_g4),filter_h5);
for jj=1:length(h5)-1
    temp_H51(jj)=W5(jj);
end
temp_H52=[W5;temp_H51];
temp_H53=xcorr(temp_H52,h5);
for kk=1:Ldata
    D5(kk)=temp_H53(kk+Ldata+length(h5)-1-1);
end

%求S5
g5=conv(conv(conv(conv(filter_g1,filter_g2),filter_g3),filter_g4),filter_g5);
for jj=1:length(g5)-1
    temp_G51(jj)=V5(jj);
end
temp_G52=[V5;temp_G51];
temp_G53=xcorr(temp_G52,g5);
for kk=1:Ldata
    S5(kk)=temp_G53(kk+Ldata+length(g5)-1-1);
end

%求D1_dn
h1=filter_h1;
for jj=1:length(h1)-1
    temp_H11(jj)=W1_dn(jj);
end
temp_H12=[W1_dn;temp_H11];
temp_H13=xcorr(temp_H12,h1);
for kk=1:Ldata
    D1_dn(kk)=temp_H13(kk+Ldata+length(h1)-1-1);
end

%求D2_dn
h2=conv(filter_g1,filter_h2);
for jj=1:length(h2)-1
    temp_H21(jj)=W2_dn(jj);
end
temp_H22=[W2_dn;temp_H21];
temp_H23=xcorr(temp_H22,h2);
for kk=1:Ldata
    D2_dn(kk)=temp_H23(kk+Ldata+length(h2)-1-1);
end

%求D3_dn
h3=conv(conv(filter_g1,filter_g2),filter_h3);
for jj=1:length(h3)-1
    temp_H31(jj)=W3_dn(jj);
end
temp_H32=[W3_dn;temp_H31];
temp_H33=xcorr(temp_H32,h3);
for kk=1:Ldata
    D3_dn(kk)=temp_H33(kk+Ldata+length(h3)-1-1);
end

%求D4_dn
h4=conv(conv(conv(filter_g1,filter_g2),filter_g3),filter_h4);
for jj=1:length(h4)-1
    temp_H41(jj)=W4_dn(jj);
end
temp_H42=[W4_dn;temp_H41];
temp_H43=xcorr(temp_H42,h4);
for kk=1:Ldata
    D4_dn(kk)=temp_H43(kk+Ldata+length(h4)-1-1);
end

%求D5_dn
h5=conv(conv(conv(conv(filter_g1,filter_g2),filter_g3),filter_g4),filter_h5);
for jj=1:length(h5)-1
    temp_H51(jj)=W5_dn(jj);
end
temp_H52=[W5_dn;temp_H51];
temp_H53=xcorr(temp_H52,h5);
for kk=1:Ldata
    D5_dn(kk)=temp_H53(kk+Ldata+length(h5)-1-1);
end

%求out
out=(((((D5_dn+S5)/2+D4_dn)/2+D3_dn)/2+D2_dn)/2+D1_dn)/2;

%输入与输出的图形
subplot(2,1,1),plot(in),axis([0,3979,-100,1000]),
subplot(2,1,2),plot(out),axis([0,3979,-100,1000]),