小波平移不变量法去噪的快速算法(附程序)

编程实现平移不变量去噪的快速算法，结合例子验证该方法的去噪性能[程序说明]：共包含三个M函数文件，分别是：shift_left函数，实现信号序列向左循环平移一位；shift_right函数，实现信号序列向右循环平移一位，这两个函数在程序实现和功能上都很简单，在此不再赘述；TI_Denoise函数，是该算法的主函数，用快速算法实现含噪信号的TI去噪。以下是该函数的帮助文档，比较清楚地说明了各参数的意义和该函数的用法：%FastTI_Denoisingof1-dsignalwithwaveletthresholding.%Usage%y=TI_Denoise(signal,wavename,L)%Inputs%signal1-dnoisysignal,length(signal)=2^J.Jmustbeanpositiveinteger.%wavenamenameofwavelet%LLow-Frequencycutoffforshrinkage.L=J.%Outputs%ythesignalafterbeingdenoised另外一个文件main调用了TI_Denoise函数，分别用不同种类的小波，对系统自带的几个含噪信号进行了平移不变量去噪，并将去噪前后的信号显示出来。其中的signal可以选用系统本身带有的几个含噪信号之一，小波wavename和分解级数L也可任意选择（L=J），以此对各种条件下的TI去噪性能进行比较。具体的M程序文件见附件。[算法流程]：TI小波去噪的功能基本上全是在TI_Denoise函数中实现的，主要流程如下：1.快速TI前向小波分解原理与课本上一致，不再赘述；这里调用了shift_left函数，来对各分辨级上的低频信号进行平移；做小波变换时，直接调用了dwt，并将延拓方式设为周期型’per’，保证了分解后系数的总数目不变；每一级dwt后的系数均直接存入TI表中，这里将每一级的低频系数也存入了TI表中，置于第一列，到下一级时再将其更新，直至最后一级的低频系数存入，则不再发生变化。2.软阈值法去噪先利用最高分辨级的小波系数对噪声的标准差做出估计，得到，再采用全局阈值对各小波系数进行软阈值法处理，这里的n为信号的长度。3.信号的重构类似TI小波分解的逆过程，在各个分辨级上，分别从TI表中取出低频系数和小波系数，用周期延拓的idwt进行重构，并用到了序列的反向平移shift_right，用每一级重构出的低频系数更新TI表的第一列，直至得到最终的重构信号。按照课本上的步骤，整个算法的思路很清晰，也并未涉及到比较复杂的运算，核心部分的操作只有两个：一个是平移，可以通过简单的循环平移来实现；另一个是正交小波变换与反变换，也可以通过直接调用dwt和idwt来实现。唯一有些麻烦的是，该算法涉及的系数比较多，对应关系比较繁琐，非常容易出错。[实验结果]：一、采用不同的小波进行去噪1.对含噪的方块信号（长度为1024），取分解级数为710，分别用Haar、Sym2、Sym8小波对其进行TI小波去噪，得到的结果如下：可见，对于这样的方块信号，用Haar小波来进行去噪，重构的效果是最好的；即便是消失矩为2的连续小波sym2，其重构的效果也要略优于消失矩为8的更加光滑的sym8。2.对含噪的Bumps信号（长度为1024），取分解级数为710，分别用Haar、Sym2、Sym8小波对其进行TI小波去噪，得到的结果如下：对于Bumps这样连续性不太好的信号，用Haar小波来进行去噪，重构的效果要稍好一点，但用各种小波来重构的差别已经很不明显。3.对含噪的Doppler信号（长度为1024），取分解级数为710，分别用Haar、Sym2、Sym8小波对其进行TI小波去噪，得到的结果如下：可见，对于Doppler这样的连续信号，用光滑性最好的Sym8小波来进行去噪，重构的效果是最好的；用Haar重构出的信号，在连续性上则要略差一些。总的来说，尽管利用TI方法对信号去噪的效果仍与所选择的小波有关，但是采用不同小波所带来的性能差异已经很不明显；与其它的小波去噪方法，比如传统的小波阈值法相比，去除了伪吉布斯现象，表现出了很好的去噪效果。二、分解级数对去噪效果的影响最后，分析了一下分解级数L对去噪效果的影响。在TI小波阈值去噪法中，并未指定分解级数L的大小，只要其满足L=J即可，其中2^J=n，n为信号长度。这里我们选用Blocks信号和Doppler信号为例，分别用Haar和Sym8小波对其进行TI阈值去噪，并取不同的分解级数L=1,4,7。这两个信号长为1024，J=10。从下图的实验结果可以看出，分解级数越多，其去噪的效果也越好，这是很好理解的，与理论分析一致。当然，在实际工程应用中，考虑到分解级数多所带来的高复杂度和长的延迟，应该选用L适中的情况，以达到去噪效果、计算复杂度以及实时性的折中和综合最优。[有关程序附件]共有四个程序，分别为：main.m，shift_left.m，shift_right.m，TI_Denoise.mmain.m-------------------------------------------------------------------------------%loadnoisysignalswhichcanbenoisbloc,noisbump,ornoisdoppetc.loadnoisbloc;signal=noisbloc;subplot(4,1,1);plot(signal);axistight;title('Noisysignal:Blocks');%FastTI_Denoisingoftheloadedsignalwavename='haar';L=1;y1=TI_Denoise(signal,wavename,L);subplot(4,1,2);plot(y1);axistight;text(450,12,'Haar,L=1');wavename='haar';L=4;y2=TI_Denoise(signal,wavename,L);subplot(4,1,3);plot(y2);axistight;text(450,12,'Haar,L=4');wavename='haar';L=7;y3=TI_Denoise(signal,wavename,L);subplot(4,1,4);plot(y3);axistight;text(450,12,'Haar,L=7');-------------------------------------------------------------------------------shift_left.m-------------------------------------------------------------------------------functiony=shift_left(x);y=[x(2:length(x))x(1)];-------------------------------------------------------------------------------shift_right.m-------------------------------------------------------------------------------functiony=shift_right(x);n=length(x);y=[x(n)x(1:n-1)];-------------------------------------------------------------------------------TI_Denoise.m-------------------------------------------------------------------------------functiony=TI_Denoise(signal,wavename,L)%FastTI_Denoisingof1-dsignalwithwaveletthresholding.%Usage%y=TI_Denoise(signal,wavename,L)%Inputs%signal1-dNoisysignal,length(Noisy)=2^J.Jmustbeanpositiveinteger.%wavenamenameofwavelet%LLow-Frequencycutoffforshrinkage.LJ%Outputs%ythesignalafterbeingdenoised%%%FastTIforwardwavelettransformn=length(signal);%InitializeTItabletable=zeros(n,L+1);table(:,1)=signal';%Differentdecomposescalesfrom1toLfori=1:L,%Sizeofthelow_frequencycoefficientsorwaveletcoefficientsatdifferentscalessize=n/2^i;%thelow_frequencycoefficientsofthepriorscalelow_coef=table(:,1);forj=0:2:2^i-2,%Getthelow_frequencycoefficientsbelta(J+1-i,j/2)sig0=low_coef(size*j+1:size*(j+2));%DWTtogetthelow_frequencycoefficientsbelta(J-i,j)andwaveletcoefficientsalpha(J-i,j)[a1,b1]=dwt(sig0,wavename,'mode','per');%StorethecoefficientsintoTItabletable((j*size+1):(j+1)*size,1)=a1;table((j*size+1):(j+1)*size,i+1)=b1;%Shift1onbelta(J+1-i,j/2)sig1=shift_left(sig0);%DWTtogetthelow_frequencycoefficientsbelta(J-i,j+1)andwaveletcoefficientsalpha(J-i,j+1)[a2,b2]=dwt(sig1,wavename,'mode','per');%StorethecoefficientsintoTItabletable(((j+1)*size+1):(j+2)*size,1)=a2;table(((j+1)*size+1):(j+2)*size,i+1)=b2;endend%%%Softthresholddenoising%estimatethestandarddeviationofadditivenoise.d=table(:,2);sigma=median(abs(d))/0.6745;%Setthethresholdthresh=sqrt(2*log(n))*sigma;%Softthresholdfori=2:L+1,forj=1:n,ifabs(table(j,i))thresh,table(j,i)=0;elseiftable(j,i)=thresh,table(j,i)=table(j,i)-thresh;elsetable(j,i)=table(j,i)+thresh;endendend%%%Reconstructthesignal%DifferentreconstructionscalesfromLto1fori=L:-1:1,size=n/2^i;forj=0:2:2^i-2,%getthelow_frequencycoefficientsbelta(J-i,j)andwa