waveletthreshod.m

小波阈值去噪

clear;clc;
%-----------------------------------------------------------------
%从程序的运行结果来看，文献1，3，4的去噪效果比较好
%其中文献4对高信噪比的的情况不是很好，在高信噪比时，软硬阈值的效果最好
%-----------------------------------------------------------------------
%测试数据的选取
% fun='blocks';
% snr=5;
% jN=5;  %分解的层数
% N=13;  %数据长度为2^N
% [x,s]=wnoise(4,N,sqrt(snr));
% ps=sum(x.^2)/length(x);
% sigma_noise=abs(sqrt(ps/(10^(snr/10))));
% noise=sigma_noise*randn(1,length(x));%noise噪声的方差是(sigma_noise.^2)
% s=x+noise;
% figure,
% subplot(211);plot(x);
% subplot(212);plot(s);
% subplot(211);plot(x);title('纯净信号x');
% subplot(212);plot(s);title('混合信号x');

%调幅信号,纯净信号
tic;
% fs=5e+6;          %采样率50M          
% ts=1/fs;             
% fc=10.7e+6;       %载频10.7M
% t0=2;             %数据长度N=t0*fs;
% t=[0:ts:t0];      %模拟信号的数字化
% m=sinc(100*t);    %消息信号
% c=cos(2*pi*fs.*t);%载波信号
% x=m.*c;
% N=t0*fs;
% toc;
%t0=.2;						    % signal duration
%ts=0.001;					    % sampling interval
%fc=250;						% carrier frequency
%snr=20;						% SNR in dB (logarithmic)
%fs=1/ts; 					    % sampling frequency
%t=[-t0/2:ts:t0/2]; 		    % time vector
%m=sinc(100*t); 				% the message signal
%c=cos(2*pi*fc.*t);				% the carrier signal
%x=m.*c;						% the DSB-AM modulated signal
%N=t0*fs;
%---------------------------------------------------------------------
%-------------------------BPSK信号,数字调制-----------------------------------
codes=6;                                %码元个数，即输入调制信号的长度
sigma=1;                                %调制信号的幅度
fs=600e3;                               %采样率600KHz
fb=1e3;                                 %波特率1KHz，fb表示对输入调制信号的采样率
fc=100e3;                               %载频100KHz
Modulate=2;                             %为选择调制方式
N0=fs/fb;                               %一个码元周期内的采样点数，一个输入信号长度内的采样点数
N=N0*codes;                             %总的采样点数（已调信号的长度）
p0=pi*rand(1,1);                        %产生初始相位
symbols=randint(1,codes,[0,1]); 	      %产生基带码元
x_B = ones(N0,1)*symbols;
x_BB = x_B(:)';                         %根据波特率要求产生码元
signal_base = x_BB;                     %产生基带信号
signal=sigma*dmod(symbols,fc,fb,[fs p0],'psk',Modulate); 
                                        %产生psk调制信号，p0是载频的初始相位
x=signal;
%-------------------------加入指定强度的噪声---------------------------------
snr=5;
ps=sum(x.^2)/N;
sigma_noise=abs(sqrt(ps/(10^(snr/10))));
nn=randn(1,N);
enn=sum(nn)/N;                          %随机数nn的均值
nn=nn-enn;                              %使nn均值为0 
noise=sigma_noise*nn;
s=x+noise;
wname='db7';
jN=6;  %分解的层数
[c,l]=wavedec(s,jN,wname);
snrs=20*log10(norm(x)/norm(s-x));
mmses=mmse(s-x);

%高频分量的索引
first = cumsum(l)+1;
first1=first;
first = first(end-2:-1:1);
ld   = l(end-1:-1:2);
last = first+ld-1;
%--------------------------------------------------------------------------
%软阈值与硬阈值
cxdsoft=c;
cxdhard=c;
for j=1:jN                 %j是分解尺度
    flk = first(j):last(j);%flk是di在c中的索引
    thr(j)=sqrt(log(length(flk)))/log(j+1);
    for k=0:(length(flk)-1)%k是位移尺度
        djk=c(first(j)+k);%为了简化程序
        djk2=c(first(j)+k);
        absdjk=abs(djk);
        thr1=thr(j);
        %阈值处理
        if absdjk<thr1
           djk=0;
       else
           djk=sign(djk)*(absdjk-thr1);   
       end  
       if abs(djk2)<thr1
           djk2=0;
       else
           djk=djk2;   
       end 
            cxdhard(first(j)+k)=djk2;
            cxdsoft(first(j)+k)=djk;
   end             
end
snewsoft=waverec(cxdsoft,l,wname);
snrysoft=20*log10(norm(x)/norm(snewsoft-x));
mmsesoft=mmse(snewsoft-x);
snewhard=waverec(cxdhard,l,wname);
snryhard=20*log10(norm(x)/norm(snewhard-x));
mmsehard=mmse(snewhard-x);
%--------------------------------------------------------------------------
%阈值处方法一
%文献【基于新阈值函数及最优尺度的小波去噪研究】
cxd1=c;
for j=1:jN                 %j是分解尺度
    flk = first(j):last(j);%flk是di在c中的索引
    thr(j)=sqrt(log(length(flk)))/log(j+1);
    for k=0:(length(flk)-1)%k是位移尺度
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        %阈值处理
        if absdjk<thr1
           djk=djk.^23/(exp(22/abs(thr1))*thr1.^22);
       else
           djk=sign(djk)*(absdjk-thr1+thr1/exp(22/absdjk));    
       end         
            cxd1(first(j)+k)=djk;
   end             
end
%信号重构
snew1=waverec(cxd1,l,wname);
figure,
subplot(211);plot(x);title('纯净信号x');
subplot(212);plot(snew1);title('采用方法一去噪后信号');

%信噪比，均方差分析
snry1=20*log10(norm(x)/norm(snew1-x));
mmsey1=mmse(snew1-x);
%--------------------------------------------------------------------------
%阈值处方法二
%文献【几种基于小波阈值去噪的改进方法】作者朱艳琴
%新方法一
cxd21=c;
a=0.5;
for j=1:jN
    flk = first(j):last(j);
    thr(j)=sqrt(2*log((j+1)/j))*median(c(flk))/0.6745;
    for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=0;
        else
            djk=sign(djk)*(abs(djk)-a*thr1);
        end
        cxd21(first(j)+k)=djk;
    end
end
%新方法二
cxd22=c;
a=0.5;
for j=1:jN
    flk = first(j):last(j);
      for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=0;
        else
            djk=sign(djk)*sqrt(abs(djk).^2-a*thr1.^2);
        end
        cxd22(first(j)+k)=djk;
    end
end
%新方法三
cxd23=c;
for j=1:jN
    flk = first(j):last(j);
    for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=0;
        else
            djk=0.5*djk+0.5*sign(djk)*(abs(djk)-a*thr1);
        end
        cxd23(first(j)+k)=djk;
    end
end
%新方法四
cxd24=c;
for j=1:jN
    flk = first(j):last(j);
    thr22(j)=sqrt(2*log(length(flk)))/(log(exp(1)+j-1))*std(c(flk));
      for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        thr2=thr22(j);
        if absdjk<thr2
            djk=0;
        elseif absdjk<thr1;
            djk=sign(djk)*thr*(absdjk-thr2)/(thr-thr2);
        else
            djk=djk;
        end
        cxd24(first(j)+k)=djk;
    end
end
%信号重构
snew21=waverec(cxd21,l,wname);
snew22=waverec(cxd22,l,wname);
snew23=waverec(cxd23,l,wname);
snew24=waverec(cxd24,l,wname);
%信噪比及均方差分析
snry21=20*log10(norm(x)/norm(snew21-x));
snry22=20*log10(norm(x)/norm(snew22-x));
snry23=20*log10(norm(x)/norm(snew23-x));
snry24=20*log10(norm(x)/norm(snew24-x));
mmsey21=mmse(snew21-x);
mmsey22=mmse(snew22-x);
mmsey23=mmse(snew23-x);
mmsey24=mmse(snew24-x);
figure,
subplot(211);plot(x);title('纯净信号x');
subplot(212);plot(snew24);title('采用方法二去噪后的信号');

%--------------------------------------------------------------------------
%阈值处方法三
%文献【小波阈值去噪的一种改进方案】
cxd3=c;
a=0.5;
for j=1:jN
    flk = first(j):last(j);
    thr(j)=sqrt(2*log(2.^N))/log(j+1);
    for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=0;
        else
            djk=sign(djk)*(abs(djk)-a*thr1+2*thr1*sqrt(a*(1-1)/exp(absdjk)));
        end
        cxd3(first(j)+k)=djk;
    end
end
%信号重构
snew3=waverec(cxd3,l,wname);
figure,
subplot(211);plot(x);title('纯净信号x');
subplot(212);plot(snew3);title('采用方法三去噪后信号');
%信噪比，均方差分析
snry3=20*log10(norm(x)/norm(snew3-x));
mmsey3=mmse(snew3-x);
%--------------------------------------------------------------------------
%阈值处方法四
%文献【基于一种新的阈值函数的小波域信号去噪】作者张维强
cxd4=c;
a=0.5;
for j=1:jN
    flk = first(j):last(j);
    thr(j)=sqrt(2*log(2.^N))/log(j+1)*median(flk)/0.6745;
    for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=0;
        else
            djk=sign(djk)*(abs(djk)-thr1/exp((absdjk-thr1)/8));
        end
        cxd4(first(j)+k)=djk;
    end
end
snew4=waverec(cxd4,l,wname);
%信噪比，均方差分析
snry4=20*log10(norm(x)/norm(snew4-x));
mmsey4=mmse(snew4-x);
figure,
subplot(211);plot(x);title('纯净信号x');
subplot(212);plot(snew4);title('采用方法四去噪后信号');
%--------------------------------------------------------------------------
%阈值处方法五
%文献【小波阈值去噪函数的改进方法分析】作者刘卫东 
n=0.1;
a=0.5;
b=1;
m=0.9;
aa=6;
%新方法一
cxd51=c;
a=0.5;
for j=1:jN
    flk = first(j):last(j);
    thr(j)=sqrt(2*log((j+1)/j))*median(c(flk))/0.6745;
    for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=0;
        elseif djk>thr1
            djk=djk-thr1/exp((djk-thr1)/n);
        else
             djk=djk+thr1/exp((-djk-thr1)/n);
        end
        cxd51(first(j)+k)=djk;
    end
end
%新方法二
cxd52=c;
for j=1:jN
    flk = first(j):last(j);
      for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=0;
        else
            djk=djk-a*thr1+2*a*thr1/(1+exp(djk));
        end
        cxd52(first(j)+k)=djk;
    end
end
%新方法三
cxd53=c;
for j=1:jN
    flk = first(j):last(j);
    for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr1=thr(j);
        if absdjk<thr1
            djk=djk.^(2*b+1)/((2*b+1)*djk.^(2*b));
        elseif djk>thr1
            djk=djk-thr1+thr1/(2*b+1);
        else
            djk=djk+thr1-thr1/(2*b+1);
        end
        cxd53(first(j)+k)=djk;
    end
end
%新方法四
cxd54=c;
for j=1:jN
    flk = first(j):last(j);
        for k=0:(length(flk)-1)
            djk=c(first(j)+k);%为了简化程序
            absdjk=abs(djk);
            thr1=thr(j);
            bb=(2*aa+1+m*sqrt((2*aa+1).^2-pi.^2))/(pi*thr1.^(2*n+1));
            kk=(1+(bb*thr1.^(2*n+1)).^2)/((2*a+1)*bb*thr1.^(2*n));
            if absdjk<thr1
                djk=kk*atan(bb*djk.^(2*n+1));
            elseif djk>thr1
                djk=djk-thr1+kk*atan(bb*djk.^(2*aa+1));
            else
                djk=djk+thr1-kk*atan(bb*djk.^(2*aa+1));
            end
        cxd54(first(j)+k)=djk;
    end
end
%信号重构
snew51=waverec(cxd51,l,wname);
snew52=waverec(cxd52,l,wname);
snew53=waverec(cxd53,l,wname);
snew54=waverec(cxd54,l,wname);
%信噪比及均方差分析
snry51=20*log10(norm(x)/norm(snew51-x));
snry52=20*log10(norm(x)/norm(snew52-x));
snry53=20*log10(norm(x)/norm(snew53-x));
snry54=20*log10(norm(x)/norm(snew54-x));
mmsey51=mmse(snew51-x);
mmsey52=mmse(snew52-x);
mmsey53=mmse(snew53-x);
mmsey54=mmse(snew54-x);
%--------------------------------------------------------------------------
%阈值处方法六
%文献【一种改进的小波域阈值去噪算法】作者付炜 
m=40;
cxd6=c;
for j=1:jN
    flk = first(j):last(j);
    thr1=thselect(s,'heursure');
    for k=0:(length(flk)-1)
        djk=c(first(j)+k);%为了简化程序
        absdjk=abs(djk);
        thr2=thr1/3;
        if absdjk>thr1
            djk=djk'*(absdjk.^m-thr1.^m)/absdjk.^m;
        elseif djk<thr2
            djk=0;
        else
             djk=djk-thr1+2*thr1/(1+exp(2*djk-thr1));
        end
        cxd6(first(j)+k)=djk;
    end
end
%信号重构
snew6=waverec(cxd6,l,wname);
%信噪比，均方差分析
snry6=20*log10(norm(x)/norm(snew6-x));
mmsey6=mmse(snew6-x);

snr1=[snrs,   snryhard, snrysoft, 0;
      snry1,  0,        0,        0;
      snry21, snry22,   snry23,   snry24;
      snry3,  snry4,    snry6,    0;
      snry51, snry52,   snry53,   snry54];

  mmse=[mmses,  mmsehard, mmsesoft, 0;
        mmsey1, 0,        0,         0;
        mmsey21,mmsey22,  mmsey23,   mmsey24;
        mmsey3, mmsey4 ,  mmsey6 ,   0;
        mmsey51,mmsey52,  mmsey53,   mmsey54];