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

编程实现平移不变量去噪的快速算法，结合例子验证该方法的去噪性能程序说明：共包含三个M 函数文件，分别是：shift_left 函数，实现信号序列向左循环平移一位；shift_right 函数，实现信号序列向右循环平移一位，这两个函数在程序实现和功能上都很简单，在此不再赘述；TI_Denoise 函数，是该算法的主函数，用快速算法实现含噪信号的TI 去噪。以下是该函数的帮助文档，比较清楚地说明了各参数的意义和该函数的用法：% Fast TI_Denoising of 1-d signal with wavelet thresholding.% Usage% y=TI_Denoise(signal,wavename,L)% Inputs% signal 1-d noisy signal, length(signal)= 2J. J must be an positive integer.% wavename name of wavelet% L Low-Frequency cutoff for shrinkage. L = J.% Outputs% y the signal after being denoised另外一个文件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 即可，其中2J=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-%load noisy signals which can be noisbloc, noisbump, or noisdopp etc.load noisbloc;signal=noisbloc;subplot(4,1,1);plot(signal);axis tight;title(Noisy signal: Blocks);%Fast TI_Denoising of the loaded signalwavename=haar;L=1;y1=TI_Denoise(signal,wavename,L);subplot(4,1,2);plot(y1);axis tight;text(450,12,Haar,L=1);wavename=haar;L=4;y2=TI_Denoise(signal,wavename,L);subplot(4,1,3);plot(y2);axis tight;text(450,12,Haar,L=4);wavename=haar;L=7;y3=TI_Denoise(signal,wavename,L);subplot(4,1,4);plot(y3);axis tight;text(450,12,Haar,L=7);-shift_left.m-function y=shift_left(x);y=x(2:length(x) x(1);-shift_right.m-function y=shift_right(x);n=length(x);y= x(n) x(1:n-1);-TI_Denoise.m-function y=TI_Denoise(signal,wavename,L)% Fast TI_Denoising of 1-d signal with wavelet thresholding.% Usage % y=TI_Denoise(signal,wavename,L)% Inputs% signal 1-d Noisy signal, length(Noisy)= 2J. J must be an positive integer.% wavename name of wavelet % L Low-Frequency cutoff for shrinkage. L J% Outputs % y the signal after being denoised%Fast TI forward wavelet transformn=length(signal);%Initialize TI tabletable=zeros(n,L+1);table(:,1)=signal;%Different decompose scales from 1 to Lfor i=1:L, %Size of the low_frequency coefficients or wavelet coefficients at different scales size=n/2i; %the low_frequency coefficients of the prior scale low_coef=table(:,1); for j=0:2:2i-2, %Get the low_frequency coefficients belta(J+1-i,j/2) sig0=low_coef(size*j+1:size*(j+2); %DWT to get the low_frequency coefficients belta(J-i,j) and wavelet coefficients alpha(J-i,j) a1,b1=dwt(sig0,wavename,mode,per); %Store the coefficients into TI table table(j*size+1):(j+1)*size,1)=a1; table(j*size+1):(j+1)*size,i+1)=b1; %Shift 1 on belta(J+1-i,j/2) sig1=shift_left(sig0); %DWT to get the low_frequency coefficients belta(J-i,j+1) and wavelet coefficients alpha(J-i,j+1) a2,b2=dwt(sig1,wavename,mode,per); %Store the coefficients into TI table table(j+1)*size+1):(j+2)*size,1)=a2; table(j+1)*size+1):(j+2)*size,i+1)=b2; endend%Soft threshold denoising%estimate the standard deviation of additive noise.d=table(:,2);sigma=median(abs(d)/0.6745;%Set the thresholdthresh=sqrt(2*log(n)*sigma;%Soft thresholdfor i=2:L+1, for j=1:n, if abs(table(j,i)=thresh, table(j,i)=table(j,i)-thresh; else table(j,i)=table(j,i)+thresh; end endend%Reconstruct the signal%Different reconstruction scales from L to 1for i=L:-1:1, size=n/2i; for j=0:2:2i-2, %get the low_frequency coefficients belta(J-i,j) and wavelet coefficients alpha(J-i,j) a1=table(j*size+1):(j+1)*size,1); b1=table(j*size+1):(j+1)*size,i+1); %IDWT sig0=idwt(a1,b1,wavename,mode,per); %get the low_frequency coefficients belta(J-i,j+1) and wavelet coefficients alpha(J-i,j+1) a2=table(j+1)*size+1):(j+2)*size,1); b2=table(j+1)*size+1):(j+2)*size,i+1); %IDWT sig1=idwt(a2,b2,waven