waveletthresholddenoise.txt

小波变换用于语音信号阈值去噪的一个matlab程序

clear
x=wavread('s5.wav');        %语音信号数据格式转换
s=x(1:16384);
subplot(311);
plot(s);

n=randn(size(s));           %生成白噪声
sn=200*s+n;                 %信号染噪
subplot(312);
plot(sn)
 
[c,l]=wavedec(sn,5,'db7');  %信号五层分解
ca5=appcoef(c,l,'db7',5);
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);
softd1=wthresh(cd1,'s',Z1);
softd2=wthresh(cd2,'s',Z2);
softd3=wthresh(cd3,'s',Z3);
softd4=wthresh(cd4,'s',Z4);
softd5=wthresh(cd5,'s',Z5);
c2=[a5 softd5 softd4 softd3 softd2 softd1];
sn1=waverec(c2,l,'db7');
subplot(313);
plot(sn1)                    
%自动搜索最佳阈值程序

d(1)=length(cd1);   %得到各层小波系数的长度
d(2)=length(cd2); 
d(3)=length(cd3); 
d(4)=length(cd4); 
d(5)=length(cd5); 
k(1)=sqrt(2*log(d(1))/d(1));   %得到各层的初始阈值
k(2)=sqrt(2*log(d(2))/d(2));
k(3)=sqrt(2*log(d(3))/d(3));
k(4)=sqrt(2*log(d(4))/d(4));
k(5)=sqrt(2*log(d(5))/d(5));
for i=1:d(1)
    t2(1)=k(1);      %计算第一层小波系数的最佳阈值
    t1=0;
    for i1=1:d(1)
        g=(cd1(i1)-t2(i)+2*t2(i)/(1+exp(2*cd1(i1)/t2(i))));
        g1=-1+2/(1+(1+exp(2*cd1(i1)/t2(i))))+4*cd1(i1)*exp(2*cd1(i1)/t2(i))/(t2(i)
            *((1+exp(2*cd1(i1)/t2(i))))^2);
        g2=8*cd1(i1)*exp(2*cd1(i1)/t2(i))/((t2(i)^2)*(1+exp(2*cd1(i1)/t2(i))))^2)
           -16*cd1(i1)*exp(4*cd1(i1)/t2(i))/((t2(i)^2)*((1+exp(2*cd1(i1)/t2(i))))^3);
        t1=t1+2*g*g1+2*g2;
    end
    if t1<=0.0000001*t2(i)
       t2(i+1)=t2(i);
       z=t2(i+1);
    else 
        t2(i+1)=t2(i)-t1;
        z=t2(i+1);
    end
        Z1(i)=cd1(i)-z+2*z/(1+exp(2*cd1(i)/z));
end

for i=1:d(2)   %计算第二层小波系数的最佳阈值 
    t2(1)=k(2);
    t1=0;
    for i1=1:d(2)
         g=(cd2(i1)-t2(i)+2*t2(i)/(1+exp(2*cd2(i1)/t2(i))));
        g1=-1+2/(1+(1+exp(2*cd2(i1)/t2(i))))+4*cd2(i1)*exp(2*cd2(i1))/(t2(i)/t2(i)
            *((1+exp(2*cd2(i1)/t2(i))))^2);
        g2=8*cd2(i1)*exp(2*cd2(i1)/t2(i))/((t2(i)^2)*(1+exp(2*cd2(i1)/t2(i))))^2)
           -16*cd2(i1)*exp(4*cd2(i1)/t2(i))/((t2(i)^2)*((1+exp(2*cd2(i1)/t2(i))))^3);
        t1=t1+2*g*g1+2*g2;
    end
    if t1<=0.0000001*t2(i)
       t2(i+1)=t2(i);
       z=t2(i+1);
    else 
        t2(i+1)=t2(i)-t1;
        z=t2(i+1);
    end
        Z2(i)=cd2(i)-z+2*z/(1+exp(2*cd2(i)/z));
end
for i=1:d(3)   %计算第三层小波系数的最佳阈值 
    t2(1)=k(3);
    t1=0;
    for i1=1:d(3)
         g=(cd3(i1)-t2(i)+2*t2(i)/(1+exp(2*cd3(i1)/t2(i))));
        g1=-1+2/(1+(1+exp(2*cd3(i1)/t2(i))))+4*cd3(i1)*exp(2*cd3(i1))/(t2(i)/t2(i)
            *((1+exp(2*cd3(i1)/t2(i))))^2);
        g2=8*cd3(i1)*exp(2*cd3(i1)/t2(i))/((t2(i)^2)*(1+exp(2*cd3(i1)/t2(i))))^2)
           -16*cd3(i1)*exp(4*cd3(i1)/t2(i))/((t2(i)^2)*((1+exp(2*cd3(i1)/t2(i))))^3);
        t1=t1+2*g*g1+2*g2;
    end
    if t1<=0.0000001*t2(i)
       t2(i+1)=t2(i);
       z=t2(i+1);
    else 
        t2(i+1)=t2(i)-t1;
        z=t2(i+1);
    end
        Z3(i)=cd3(i)-z+2*z/(1+exp(2*cd3(i)/z));
end
for i=1:d(4)   %计算第四层小波系数的最佳阈值 
    t2(1)=k(4);
    t1=0;
    for i1=1:d(4)
         g=(cd4(i1)-t2(i)+2*t2(i)/(1+exp(2*cd4(i1)/t2(i))));
        g1=-1+2/(1+(1+exp(2*cd4(i1)/t2(i))))+4*cd4(i1)*exp(2*cd4(i1))/(t2(i)/t2(i)
            *((1+exp(2*cd4(i1)/t2(i))))^2);
        g2=8*cd4(i1)*exp(2*cd4(i1)/t2(i))/((t2(i)^2)*(1+exp(2*cd4(i1)/t2(i))))^2)
           -16*cd4(i1)*exp(4*cd4(i1)/t2(i))/((t2(i)^2)*((1+exp(2*cd4(i1)/t2(i))))^3);
        t1=t1+2*g*g1+2*g2;
    end
    if t1<=0.0000001*t2(i)
       t2(i+1)=t2(i);
       z=t2(i+1);
    else 
        t2(i+1)=t2(i)-t1;
        z=t2(i+1);
    end
        Z4(i)=cd2(i)-z+2*z/(1+exp(2*cd4(i)/z));
end
for i=1:d(5)   %计算第五层小波系数的最佳阈值 
    t2(1)=k(5);
    t1=0;
    for i1=1:d(2)
         g=(cd5(i1)-t2(i)+2*t2(i)/(1+exp(2*cd5(i1)/t2(i))));
        g1=-1+2/(1+(1+exp(2*cd5(i1)/t2(i))))+4*cd5(i1)*exp(2*cd5(i1))/(t2(i)/t2(i)
            *((1+exp(2*cd5(i1)/t2(i))))^2);
        g2=8*cd5(i1)*exp(2*cd5(i1)/t2(i))/((t2(i)^2)*(1+exp(2*cd5(i1)/t2(i))))^2)
           -16*cd5(i1)*exp(4*cd5(i1)/t2(i))/((t2(i)^2)*((1+exp(2*cd5(i1)/t2(i))))^3);
        t1=t1+2*g*g1+2*g2;
    end
    if t1<=0.0000001*t2(i)
       t2(i+1)=t2(i);
       z=t2(i+1);
    else 
        t2(i+1)=t2(i)-t1;
        z=t2(i+1);
    end
        Z5(i)=cd2(i)-z+2*z/(1+exp(2*cd5(i)/z));
end