小波平移不變量法去噪的快速算法(附程序)小波平移不變量法去噪的快速算法(附程序)小波平移不變量法去噪的快速算法(附程序)xxx公司小波平移不變量法去噪的快速算法(附程序)文件編號(hào)：文件日期：修訂次數(shù)：第1.0次更改批準(zhǔn)審核制定方案設(shè)計(jì)，管理制度編程實(shí)現(xiàn)平移不變量去噪的快速算法，結(jié)合例子驗(yàn)證該方法的去噪性能[程序說明]：共包含三個(gè)M函數(shù)文件，分別是：shift_left函數(shù)，實(shí)現(xiàn)信號(hào)序列向左循環(huán)平移一位；shift_right函數(shù)，實(shí)現(xiàn)信號(hào)序列向右循環(huán)平移一位，這兩個(gè)函數(shù)在程序?qū)崿F(xiàn)和功能上都很簡單，在此不再贅述；TI_Denoise函數(shù)，是該算法的主函數(shù)，用快速算法實(shí)現(xiàn)含噪信號(hào)的TI去噪。以下是該函數(shù)的幫助文檔，比較清楚地說明了各參數(shù)的意義和該函數(shù)的用法：%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另外一個(gè)文件main調(diào)用了TI_Denoise函數(shù)，分別用不同種類的小波，對(duì)系統(tǒng)自帶的幾個(gè)含噪信號(hào)進(jìn)行了平移不變量去噪，并將去噪前后的信號(hào)顯示出來。其中的signal可以選用系統(tǒng)本身帶有的幾個(gè)含噪信號(hào)之一，小波wavename和分解級(jí)數(shù)L也可任意選擇（L<=J），以此對(duì)各種條件下的TI去噪性能進(jìn)行比較。具體的M程序文件見附件。[算法流程]：TI小波去噪的功能基本上全是在TI_Denoise函數(shù)中實(shí)現(xiàn)的，主要流程如下：1.快速TI前向小波分解原理與課本上一致，不再贅述；這里調(diào)用了shift_left函數(shù)，來對(duì)各分辨級(jí)上的低頻信號(hào)進(jìn)行平移；做小波變換時(shí)，直接調(diào)用了dwt，并將延拓方式設(shè)為周期型’per’，保證了分解后系數(shù)的總數(shù)目不變；每一級(jí)dwt后的系數(shù)均直接存入TI表中，這里將每一級(jí)的低頻系數(shù)也存入了TI表中，置于第一列，到下一級(jí)時(shí)再將其更新，直至最后一級(jí)的低頻系數(shù)存入，則不再發(fā)生變化。2.軟閾值法去噪先利用最高分辨級(jí)的小波系數(shù)對(duì)噪聲的標(biāo)準(zhǔn)差做出估計(jì)，得到，再采用全局閾值對(duì)各小波系數(shù)進(jìn)行軟閾值法處理，這里的n為信號(hào)的長度。3.信號(hào)的重構(gòu)類似TI小波分解的逆過程，在各個(gè)分辨級(jí)上，分別從TI表中取出低頻系數(shù)和小波系數(shù)，用周期延拓的idwt進(jìn)行重構(gòu)，并用到了序列的反向平移shift_right，用每一級(jí)重構(gòu)出的低頻系數(shù)更新TI表的第一列，直至得到最終的重構(gòu)信號(hào)。按照課本上的步驟，整個(gè)算法的思路很清晰，也并未涉及到比較復(fù)雜的運(yùn)算，核心部分的操作只有兩個(gè)：一個(gè)是平移，可以通過簡單的循環(huán)平移來實(shí)現(xiàn)；另一個(gè)是正交小波變換與反變換，也可以通過直接調(diào)用dwt和idwt來實(shí)現(xiàn)。唯一有些麻煩的是，該算法涉及的系數(shù)比較多，對(duì)應(yīng)關(guān)系比較繁瑣，非常容易出錯(cuò)。[實(shí)驗(yàn)結(jié)果]：一、采用不同的小波進(jìn)行去噪1.對(duì)含噪的方塊信號(hào)（長度為1024），取分解級(jí)數(shù)為7<10，分別用Haar、Sym2、Sym8小波對(duì)其進(jìn)行TI小波去噪，得到的結(jié)果如下：可見，對(duì)于這樣的方塊信號(hào)，用Haar小波來進(jìn)行去噪，重構(gòu)的效果是最好的；即便是消失矩為2的連續(xù)小波sym2，其重構(gòu)的效果也要略優(yōu)于消失矩為8的更加光滑的sym8。2.對(duì)含噪的Bumps信號(hào)（長度為1024），取分解級(jí)數(shù)為7<10，分別用Haar、Sym2、Sym8小波對(duì)其進(jìn)行TI小波去噪，得到的結(jié)果如下：對(duì)于Bumps這樣連續(xù)性不太好的信號(hào)，用Haar小波來進(jìn)行去噪，重構(gòu)的效果要稍好一點(diǎn)，但用各種小波來重構(gòu)的差別已經(jīng)很不明顯。3.對(duì)含噪的Doppler信號(hào)（長度為1024），取分解級(jí)數(shù)為7<10，分別用Haar、Sym2、Sym8小波對(duì)其進(jìn)行TI小波去噪，得到的結(jié)果如下：可見，對(duì)于Doppler這樣的連續(xù)信號(hào)，用光滑性最好的Sym8小波來進(jìn)行去噪，重構(gòu)的效果是最好的；用Haar重構(gòu)出的信號(hào)，在連續(xù)性上則要略差一些?？偟膩碚f，盡管利用TI方法對(duì)信號(hào)去噪的效果仍與所選擇的小波有關(guān)，但是采用不同小波所帶來的性能差異已經(jīng)很不明顯；與其它的小波去噪方法，比如傳統(tǒng)的小波閾值法相比，去除了偽吉布斯現(xiàn)象，表現(xiàn)出了很好的去噪效果。二、分解級(jí)數(shù)對(duì)去噪效果的影響最后，分析了一下分解級(jí)數(shù)L對(duì)去噪效果的影響。在TI小波閾值去噪法中，并未指定分解級(jí)數(shù)L的大小，只要其滿足L<=J即可，其中2^J=n，n為信號(hào)長度。這里我們選用Blocks信號(hào)和Doppler信號(hào)為例，分別用Haar和Sym8小波對(duì)其進(jìn)行TI閾值去噪，并取不同的分解級(jí)數(shù)L=1,4,7。這兩個(gè)信號(hào)長為1024，J=10。從下圖的實(shí)驗(yàn)結(jié)果可以看出，分解級(jí)數(shù)越多，其去噪的效果也越好，這是很好理解的，與理論分析一致。當(dāng)然，在實(shí)際工程應(yīng)用中，考慮到分解級(jí)數(shù)多所帶來的高復(fù)雜度和長的延遲，應(yīng)該選用L適中的情況，以達(dá)到去噪效果、計(jì)算復(fù)雜度以及實(shí)時(shí)性的折中和綜合最優(yōu)。[有關(guān)程序附件]共有四個(gè)程序，分別為：，，，-------------------------------------------------------------------------------%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');--------------------------------------------------------------------------------------------------------------------------------------------------------------functiony=shift_left(x);y=[x(2:length(x))x(1)];--------------------------------------------------------------------------------------------------------------------------------------------------------------functiony=shift_right(x);n=length(x);y=[x(n)x(1:n-1)];--------------------------------------------------------------------------------------------------------------------------------------------------------------functiony=TI_Denoise(signal,wavename,L)%FastTI_Denoisingof1-dsignalwithwaveletthresholding.%Usage%y=TI_Denoise(signal,wavename,L)%Inputs%signal 1-dNoisysignal,length(Noisy)=2^J.Jmustbeanpositiveinteger.%wavenamenameofwavelet%L Low-Frequencycutoffforshrinkage.L<<J%Outputs%y thesignalafterbeingdenoised%%%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))/;%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)andwaveletcoefficientsalpha(J-i,j)a1=table((j*size+1):(j+1)*size,1);b1=table((j*size+1):(j+1)*size,i+1);%IDWTsig0=idwt(a1,