基于Matlab的小波分解、去噪與重構(gòu)

1、現(xiàn)代信號(hào)處理大作業(yè)基于Matlab的小波分解、去噪與重構(gòu)目錄一 作業(yè)內(nèi)容及要求31.1 作業(yè)內(nèi)容31.2 作業(yè)要求3二 系統(tǒng)原理32.1 小波變換原理32.2 閾值去噪原理3三 系統(tǒng)分析及設(shè)計(jì)53.1 圖像分解53.2 高頻去噪53.3 圖像重構(gòu)6四 程序編寫74.1 main函數(shù)74.2 分解函數(shù)94.2.1 二維分解函數(shù)94.2.2 一維分解函數(shù)104.3 卷積函數(shù)104.4 采樣函數(shù)114.4.1 下采樣函數(shù)114.4.2 上采樣函數(shù)114.5 重構(gòu)函數(shù)124.5.1 二維重構(gòu)函數(shù)124.5.2 一維重構(gòu)函數(shù)13五 結(jié)果分析及檢驗(yàn)145.1 結(jié)果分析145.2 結(jié)果檢驗(yàn)16六 心得體會(huì)1

2、8參考文獻(xiàn)19一 作業(yè)內(nèi)容及要求1.1 作業(yè)內(nèi)容用小波對(duì)圖像進(jìn)行濾波分解、去噪，然后重構(gòu)。1.2 作業(yè)要求用小波對(duì)圖像進(jìn)行濾波分解、去噪，然后重構(gòu)。具體要求：(1) 被處理圖像可選擇：woman, wbarb, wgatlin, detfingr, tire.；(2) 可以選擇db等正交小波、或雙正交小波（或用幾種小波）；(3) 用選用小波的分解濾波器通過定義的卷積函數(shù)conv_my( )對(duì)圖像二維數(shù)組進(jìn)行小波分解，并進(jìn)行下采樣，獲取CA、CV、CD、CH等分解子圖；(4) 對(duì)高頻信號(hào)子圖進(jìn)行去噪處理，可以采用軟閾值、硬閾值等方法；(5) 用選用小波的綜合濾波器對(duì)去噪的子圖進(jìn)行圖像重構(gòu)。二 系

3、統(tǒng)原理2.1 小波變換原理小波變換的一級(jí)分解過程是，先將信號(hào)與低通濾波器卷積再下采樣可以得到低頻部分的小波分解系數(shù)再將信號(hào)與高通濾波器卷積后下采樣得到高頻部分的小波分解系數(shù)；而多級(jí)分解則是對(duì)上一級(jí)分解得到的低頻系數(shù)再進(jìn)行小波分解，是一個(gè)遞歸過程。二維小波分解重構(gòu)可以用一系列的一維小波分解重構(gòu)來實(shí)現(xiàn)。重構(gòu)則是分解的逆過程，對(duì)低頻系數(shù)、高頻系數(shù)分別進(jìn)行上抽樣和低通、高通濾波處理。要注意重構(gòu)時(shí)同一級(jí)的低頻、高頻系數(shù)的個(gè)數(shù)必須相等。2.2 閾值去噪原理圖像去噪的方法是：（1）圖像的小波分解。選擇合適的小波函數(shù)以及合適的分解層數(shù)對(duì)圖像進(jìn)行分解。（2）對(duì)分解后的高頻系數(shù)進(jìn)行閾值處理。對(duì)分解的每一層，選擇合

4、適的閾值對(duì)該層的水平、垂直和對(duì)角三個(gè)方向的高頻系數(shù)進(jìn)行閾值處理。（3）重構(gòu)圖像。根據(jù)小波分解的低頻系數(shù)和經(jīng)閾值量化處理后的高頻系數(shù)進(jìn)行圖像重構(gòu)。本設(shè)計(jì)采用軟閾值去噪，其原理為：當(dāng)小波系數(shù)的絕對(duì)值小于給定的閾值時(shí)，令其為0，大于閾值時(shí)，令其都減去閾值，即：小波閾值在去噪過程中起到?jīng)Q定性的作用。如果太小，那么施加閾值以后的小波系數(shù)中將包含過多的噪聲分量，達(dá)不到去噪的目的；反之，如果太大，那么將去除一部分信號(hào)的分量，從而使由小波系數(shù)重建后的信號(hào)產(chǎn)生過大的失真。MATLAB中實(shí)現(xiàn)閾值獲取的函數(shù)有ddenmp、thselect、wbmpen和wdcbm，本次設(shè)計(jì)中采用ddenmp函數(shù)進(jìn)行閾值獲取。調(diào)用格

5、式為THR,SORH,KEEPAPP=ddencmp('den','wv',X)，函數(shù)ddencmp用于獲取信號(hào)在消噪或壓縮過程中的默認(rèn)閾值。輸入?yún)?shù)X為一維或二維信號(hào)；'den'表示進(jìn)行去噪； 'wv'表示選擇小波。返回值THR是返回的閾值；SORH是軟閾值或硬閾值選擇參數(shù)；KEEPAPP表示保存低頻信號(hào)。三 系統(tǒng)分析及設(shè)計(jì)3.1 圖像分解圖像分解在程序中分為兩部分：一維分解以及二維分解。（1）一維分解先將信號(hào)與低通濾波器卷積再下采樣可以得到低頻部分的小波分解系數(shù)，再將信號(hào)與高通濾波器卷積后下采樣得到高頻部分的小波分解系數(shù)（2）

6、對(duì)輸入的信號(hào)，也就是前面load圖像之后得到的X矩陣，先對(duì)每一行進(jìn)行一維分解，在對(duì)分解后得到的矩陣的每一列進(jìn)行一維分解，最后得到一個(gè)矩陣X2=CACHCVCD；可以將CA、CH、 CV、 CD依次輸出得到分解后的低頻、垂直、水平、對(duì)角分量。原始信號(hào)分別進(jìn)行低通、高通濾波，再分別對(duì)列進(jìn)行二元下抽樣，就得到低頻、高頻（也稱為平均、細(xì)節(jié)）兩部分系數(shù)；再對(duì)這兩部分系數(shù)進(jìn)行小波分解，步驟與第一級(jí)分解一致，最后得到cA：低頻分量，cH：水平分量，cV：垂直分量，cD：對(duì)角線分量。流程圖如下圖1所示。列卷積XLo_DHi_D22Lo_DHi_DHi_DLo_D2222CACHCDCV行卷積行卷積列卷積列卷積

7、列卷積圖1 分解流程圖3.2 高頻去噪本次設(shè)計(jì)采用軟閾值去噪方法，將對(duì)分解出的cH、cV、cD進(jìn)行閾值去噪，選用ddencmp函數(shù)獲取去噪過程中的默認(rèn)軟閾值。其調(diào)用格式為：THR,SORH,KEEPAPP=ddencmp(IN1,IN2,X)。IN1='den'，IN2='wv'，X為一維或二維矩陣信號(hào)。在獲取默認(rèn)閾值后，對(duì)該閾值進(jìn)行軟閾值處理，調(diào)用到函數(shù)wthcoef2進(jìn)行處理，其調(diào)用過程中關(guān)鍵部分是 tmp = (abs(x)-t); tmp = (tmp+abs(tmp)/2; y = sign(x).*tmp;其中t代表前面獲得的默認(rèn)閾值，x代表小波系數(shù)

8、的大小。得到的y即為軟閾值處理后的小波系數(shù)。3.3 圖像重構(gòu)重構(gòu)的過程與分解過程類似，同樣分為兩部分：一維重構(gòu)和二維重構(gòu)。(1)一維重構(gòu)：先對(duì)平均部分系數(shù)進(jìn)行上采樣，再與低通濾波器卷積。然后對(duì)細(xì)節(jié)部分選取重構(gòu)所需的細(xì)節(jié)部分，長度與本層平均部分系數(shù)相同，再對(duì)其上采樣后進(jìn)行高通濾波器卷積。將兩個(gè)得到的結(jié)果相加得到新的平均部分系數(shù)。重復(fù)以上操作知道細(xì)節(jié)部分長度小于平均部分長度。(2)二維重構(gòu)：根據(jù)重構(gòu)流程圖可知，先對(duì)cA和cH進(jìn)行重構(gòu)，再對(duì)cV和cD重構(gòu)，在該程序中所以在該程序中X2=CACHCVCD，所以先對(duì)行左右兩半部分進(jìn)行重構(gòu)，重構(gòu)后得到的矩陣在進(jìn)行列上下兩半部分進(jìn)行重構(gòu)，最終得到的重構(gòu)圖像。

9、列卷積具體步驟是先將去噪后的子圖分別對(duì)行進(jìn)行，再將低頻分量和垂直分量進(jìn)行二元上采樣、低通濾波并合成；將水平分量和對(duì)角線分量進(jìn)行二元上采樣、高通濾波并合成，就得到低頻、高頻（也稱為平均、細(xì)節(jié)）兩部分系數(shù)；再對(duì)這兩部分系數(shù)進(jìn)行小波合成，步驟與第一級(jí)重構(gòu)一致，最后通過wkeep選取我們所需的原始圖像的大小，得到原始圖像。該流程圖如下圖2所示。2Lo_RHi_RHi_RLo_R列卷積列卷積列卷積CACHCDCV22222Lo_RLo_R行卷積行卷積WkeepX圖2 重構(gòu)流程圖四 程序編寫4.1 main函數(shù)clear all;% db4小波分解成低通和高通濾波器%figure;Lo_D,Hi_D,Lo

10、_R,Hi_R=wfilters('db4');subplot(221);stem(Lo_D);grid on;title('分解低通濾波器');xlabel('系數(shù)');ylabel('Lo_D');subplot(222);stem(Hi_D);grid on;title('分解高通濾波器');xlabel('系數(shù)');ylabel('Hi_D');subplot(223);stem(Lo_R);grid on;title('綜合低通濾波器');xlabel(&#

11、39;系數(shù)');ylabel('Lo_R');subplot(224);stem(Hi_R);grid on;title('綜合高通濾波器');xlabel('系數(shù)');ylabel('Hi_R');% 導(dǎo)入和顯示原始圖像%load woman;photo=X; %X含有被裝載的信號(hào)figure;colormap(map); %設(shè)置顏色映射到矩陣圖image(photo);axis off; %關(guān)閉原始信號(hào)的坐標(biāo)軸title('原始信號(hào)');% 用自定義的一維和二維函數(shù)分解圖像% 得到低頻部分（cA），水平

12、部分（cH），垂直部分（cV），對(duì)角線部分（cD）子圖%cA,cH,cV,cD,x1=mydwt2(photo,Lo_D,Hi_D,map); % 顯示經(jīng)過行分解和下采樣之后的圖像%n3,n4=size(x1);L0=x1(1:n3,1:n4/2);H0=x1(1:n3,n4/2+1:n4);figure;colormap(map);image(L0);title('行分解后的低頻分量');figure;colormap(map);image(H0);title('行分解后的高頻分量');figure;% 對(duì)高頻部分（cH,cV,cD）進(jìn)行去噪處理% 得到新的高

13、頻子圖（new_cH,new_cV,new_cD）%用小波函數(shù)sym5對(duì)cH進(jìn)行1層小波分解C1,S1=wavedec2(cH,1,'sym5'); %二維小波多層分解，這里取1，即為1層，輸出C1為各層分解系數(shù)，S1為各層分解系數(shù)長度，下同L=size(x1);%求取信號(hào)去噪的默認(rèn)閾值、軟閾值、并保留低頻部分thr1,sorh1,keepapp1=ddencmp('den','wv',cH);%對(duì)三個(gè)方向高頻系數(shù)進(jìn)行軟閾值處理%wthcoef2的調(diào)用格式為NC =wthcoef2('type',C,S,N,T,SORH)%其中&

14、#39;type'分別取'h','v','d'代表三個(gè)方向%N為尺度向量，T為閾值向量，SORH取'h'和's'分別表示使用硬閾值或軟閾值，下同nc1=wthcoef2('h',C1,S1,1,thr1,'s');nc1=wthcoef2('v',C1,S1,1,thr1,'s');nc1=wthcoef2('d',C1,S1,1,thr1,'s');new_cH=waverec2(nc1,S1,'sym

15、5'); %對(duì)新的小波分解結(jié)構(gòu)(nc，s)進(jìn)行重構(gòu)new_cHH=wkeep(new_cH,L,'c'); %提取重構(gòu)的圖像，'c'代表中心部分，但因?yàn)長為完整長度，所以用'l'或'r'也可以，下同%用小波函數(shù)sym5對(duì)cV進(jìn)行1層小波分解C2,S2=wavedec2(cV,1,'sym5');%求取信號(hào)去噪的默認(rèn)閾值、軟閾值、并保留低頻部分thr2,sorh2,keepapp2=ddencmp('den','wv',cV);%對(duì)三個(gè)方向高頻系數(shù)進(jìn)行軟閾值處理nc2=wthc

16、oef2('h',C2,S2,1,thr2,'s');nc2=wthcoef2('v',C2,S2,1,thr2,'s');nc2=wthcoef2('d',C2,S2,1,thr2,'s');new_cV=waverec2(nc2,S2,'sym5'); %對(duì)新的小波分解結(jié)構(gòu)(nc，s)進(jìn)行重構(gòu)new_cVV=wkeep(new_cV,L,'c'); %提取重構(gòu)的圖像%用小波函數(shù)sym5對(duì)cD進(jìn)行1層小波分解C3,S3=wavedec2(cD,1,'sym5

17、');%求取信號(hào)去噪的默認(rèn)閾值、軟閾值、并保留低頻部分thr3,sorh3,keepapp3=ddencmp('den','wv',cD);%對(duì)三個(gè)方向高頻系數(shù)進(jìn)行軟閾值處理nc3=wthcoef2('h',C3,S3,1,thr3,'s');nc3=wthcoef2('v',C3,S3,1,thr3,'s');nc3=wthcoef2('d',C3,S3,1,thr3,'s');new_cD=waverec2(nc3,S3,'sym5');

18、%對(duì)新的小波分解結(jié)構(gòu)(nc，s)進(jìn)行重構(gòu)new_cDD=wkeep(new_cD,L,'c'); %提取重構(gòu)的圖像%subplot(221);colormap(map);image(new_cHH);title('去噪后的水平分量');subplot(222);colormap(map);image(new_cVV);title('去噪后的垂直分量');subplot(223);colormap(map);image(new_cDD);title('去噪后的對(duì)角線分量'); % cA和去噪后的子圖（new_cH,new_cV,ne

19、w_cD）進(jìn)行重構(gòu) % reconstruction為重構(gòu)之后的圖形%reconstruction = myidwt2(cA,new_cH,new_cV,new_cD,Lo_R,Hi_R,map);%4.2 分解函數(shù)4.2.1 二維分解函數(shù)function LL,HL,LH,HH,x1=mydwt2(x,Lo_D,Hi_D,map)% 函數(shù) MYDWT2() 對(duì)輸入的r*c維矩陣 x 進(jìn)行二維小波分解，輸出四個(gè)分解系數(shù)子矩陣LL,HL,LH,HH% 輸入?yún)?shù)：x 輸入矩陣，為r*c維矩陣。% 輸出參數(shù)：LL,HL,LH,HH 是分解系數(shù)矩陣的四個(gè)相等大小的子矩陣,大小均為 r/2 * c/2

20、維% LL：低頻部分分解系數(shù)； HL：垂直方向分解系數(shù)；% LH：水平方向分解系數(shù)； HH：對(duì)角線方向分解系數(shù)。 row,col=size(x); % 讀取輸入矩陣的大小for j=1:row % 首先對(duì)輸入矩陣的每一行序列進(jìn)行一維離散小波分解 tmp1=x(j,:); ca1,cd1,n1,n2=mydwt(tmp1,Lo_D,Hi_D,1); a1=n1+n2; x1(j,:)=ca1,cd1; % 將分解系數(shù)序列再存入矩陣x1中，得到L|HendL=size(x1); %取x1的長度 for k=1:a1 % 再對(duì)輸入矩陣的每一列序列進(jìn)行一維離散小波分解 tmp2=x1(:,k); ca

21、2,cd2,n1,n2=mydwt(tmp2,Lo_D,Hi_D,1); x2(:,k)=ca2,cd2; % 將分解所得系數(shù)存入矩陣x2中，得到LL,HL;LH,HHendrow,col=size(x2); LL=x2(1:row/2,1:col/2); % LL是矩陣x2的左上角部分LH=x2(row/2+1:row,1:col/2); % LH是矩陣x2的左下角部分HL=x2(1:row/2,col/2+1:col); % HL是矩陣x2的右上角部分HH=x2(row/2+1:row,col/2+1:col); % HH是矩陣x2的右下角部分figure;subplot(221);col

22、ormap(map);image(LL);title('分解后的低頻分量');subplot(222);colormap(map);image(LH);title('分解后的垂直分量');subplot(223);colormap(map);image(HL);title('分解后的水平分量');subplot(224);colormap(map);image(HH);title('分解后的對(duì)角分量');4.2.2 一維分解函數(shù)function cA,cD,n1,n2 = mydwt(x,Lo_D,Hi_D,dim)% 函數(shù) cA

23、,cD,n1,n2=MYDWT(X,LPD,HPD,DIM) 對(duì)輸入序列x進(jìn)行一維離散小波分解，輸出分解序列cA,cD% 輸入?yún)?shù)：x輸入序列；% lpd低通濾波器；% hpd高通濾波器；% dim小波分解級(jí)數(shù)。% 輸出參數(shù)：cA低頻部分的小波分解系數(shù)；% cD高頻部分的小波分解系數(shù)。cA=x; % 初始化cA，cDcD=;for i=1:dim cvl=myconv1(cA,Lo_D); % 低通濾波 dnl=downspl(cvl); % 通過下采樣求出平均部分的分解系數(shù) n1=length(dnl); cvh=myconv1(cA,Hi_D); % 高通濾波 dnh=downspl(cv

24、h); % 通過下采樣求出本層分解后的細(xì)節(jié)部分系數(shù) n2=length(dnh); cA=dnl; % 下采樣后的平均部分系數(shù)進(jìn)入下一層分解 cD=cD,dnh; % 將本層分解所得的細(xì)節(jié)部分系數(shù)存入序列cDend4.3 卷積函數(shù)function c=myconv1(a,b);m=length(a); %獲取a矩陣的長度n=length(b);%獲取b矩陣的長度for i=1:m+n-1 c(i)=0;%定義長度為m+n-1的全為0的矩陣cendfor i=1:m for j=1:n x=a(i)*b(j);%將a矩陣中的數(shù)乘以b矩陣中的數(shù) c(i+j-1)=c(i+j-1)+x;%將上一行得

25、到的數(shù)累加到對(duì)應(yīng)序號(hào)的c矩陣數(shù)中endend4.4 采樣函數(shù)4.4.1 下采樣函數(shù)function y=downspl(x)% 函數(shù) Y=DOWMSPL(X) 對(duì)輸入序列進(jìn)行下抽樣，輸出序列 Y。% 下抽樣是對(duì)輸入序列取其偶數(shù)位，舍棄奇數(shù)位。例如 x=x1,x2,x3,x4,x5，則 y=x2,x4. N=length(x); % 讀取輸入序列長度M=floor(N/2); % 輸出序列的長度是輸入序列長度的一半（帶小數(shù)時(shí)取整數(shù)部分）i=1:M;y(i)=x(2*i);4.4.2 上采樣函數(shù)function y=upspl(x)% 函數(shù) Y = UPSPL(X) 對(duì)輸入的一維序列x進(jìn)行上抽樣，

26、即對(duì)序列x每個(gè)元素之間% 插零，例如 x=x1,x2,x3,x4,上抽樣后為 y=x1,0,x2,0,x3,0,x4;N=length(x); % 讀取輸入序列長度M=2*N-1; % 輸出序列的長度是輸入序列長度的2倍再減一for i=1:M % 輸出序列的偶數(shù)位為0，奇數(shù)位按次序等于相應(yīng)位置的輸入序列元素 if mod(i,2) y(i)=x(i+1)/2); else y(i)=0; endend4.5 重構(gòu)函數(shù)4.5.1 二維重構(gòu)函數(shù)function y=myidwt2(LL,HL,LH,HH,Lo_R,Hi_R,map)% 函數(shù) MYIDWT2() 對(duì)輸入的子矩陣序列進(jìn)行逆小波變換，

27、重構(gòu)出矩陣 y% 輸入?yún)?shù)：LL,HL,LH,HH 是四個(gè)大小均為 r*c 維的矩陣% 輸出參數(shù)：y 是一個(gè)大小為 2r*2c 維的矩陣tmp_mat=LL,HL;LH,HH; % 將輸入的四個(gè)矩陣組合為一個(gè)矩陣flag=0; %循環(huán)標(biāo)志位for i=1:2 if flag=1 row,col=size(yt); % 求出組合矩陣的行列數(shù) for k=1:col % 首先對(duì)組合矩陣tmp_mat的每一列，分開成上下兩半 ca1=yt(1:row/2,k); % 分開的兩部分分別作為平均系數(shù)序列ca1、細(xì)節(jié)系數(shù)序列cd1 cd1=yt(row/2+1:row,k); tmp1=myidwt(ca

28、1,cd1,Lo_R,Hi_R);% 重構(gòu)序列 yt1(:,k)=tmp1; % 將重構(gòu)序列存入待輸出矩陣 yt 的相應(yīng)列，此時(shí) y=L|H end break; end if flag=0 row,col=size(tmp_mat); for j=1:row % 將輸出矩陣 y 的每一行，分開成左右兩半 ca2=tmp_mat(j,1:col/2); % 分開的兩部分分別作為平均系數(shù)序列ca2、細(xì)節(jié)系數(shù)序列cd2 cd2=tmp_mat(j,col/2+1:col); tmp2=myidwt(ca2,cd2,Lo_R,Hi_R); % 重構(gòu)序列 yt(j,:)=tmp2; % 將重構(gòu)序列存入

29、待輸出矩陣 yt 的相應(yīng)行，得到最終的輸出矩陣 y=yt end flag=1; endendy=yt1;L=256,256; %取數(shù)據(jù)原始大小figure;y=wkeep(yt1,L,'c') colormap(map);image(y);title('重構(gòu)圖像');4.5.2 一維重構(gòu)函數(shù)function y=myidwt(cA,cD,Lo_R,Hi_R)% 函數(shù) MYIDWT() 對(duì)輸入的小波分解系數(shù)進(jìn)行逆離散小波變換，重構(gòu)出信號(hào)序列 y% 輸入?yún)?shù)：cA 平均部分的小波分解系數(shù)；% cD 細(xì)節(jié)部分的小波分解系數(shù)；% Lo_R、Hi_R 重構(gòu)所用的低通、高通濾波器。 lca=length(cA); % 求出平均、細(xì)節(jié)部分分解系數(shù)的長度lcd=length(cD); while (lcd)>=(lca) % 每一層重構(gòu)中，cA 和 cD 的長度要相等，故每層重構(gòu)后， % 若lcd小于lca，則重構(gòu)停止