基于Matlab的小波分解、去噪與重構

1、現(xiàn)代信號處理大作業(yè)基于Matlab的小波分解、去噪與重構目錄一 作業(yè)內(nèi)容及要求31.1 作業(yè)內(nèi)容31.2 作業(yè)要求3二 系統(tǒng)原理32.1 小波變換原理32.2 閾值去噪原理3三 系統(tǒng)分析及設計53.1 圖像分解53.2 高頻去噪53.3 圖像重構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 重構函數(shù)124.5.1 二維重構函數(shù)124.5.2 一維重構函數(shù)13五 結果分析及檢驗145.1 結果分析145.2 結果檢驗16六 心得體會1

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

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

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

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

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

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

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

9、列卷積具體步驟是先將去噪后的子圖分別對行進行，再將低頻分量和垂直分量進行二元上采樣、低通濾波并合成；將水平分量和對角線分量進行二元上采樣、高通濾波并合成，就得到低頻、高頻（也稱為平均、細節(jié)）兩部分系數(shù)；再對這兩部分系數(shù)進行小波合成，步驟與第一級重構一致，最后通過wkeep選取我們所需的原始圖像的大小，得到原始圖像。該流程圖如下圖2所示。2Lo_RHi_RHi_RLo_R列卷積列卷積列卷積CACHCDCV22222Lo_RLo_R行卷積行卷積WkeepX圖2 重構流程圖四 程序編寫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');% 導入和顯示原始圖像%load woman;photo=X; %X含有被裝載的信號figure;colormap(map); %設置顏色映射到矩陣圖image(photo);axis off; %關閉原始信號的坐標軸title('原始信號');% 用自定義的一維和二維函數(shù)分解圖像% 得到低頻部分（cA），水平

12、部分（cH），垂直部分（cV），對角線部分（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;% 對高頻部分（cH,cV,cD）進行去噪處理% 得到新的高

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

14、#39;type'分別取'h','v','d'代表三個方向%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'); %對新的小波分解結構(nc，s)進行重構new_cHH=wkeep(new_cH,L,'c'); %提取重構的圖像，'c'代表中心部分，但因為L為完整長度，所以用'l'或'r'也可以，下同%用小波函數(shù)sym5對cV進行1層小波分解C2,S2=wavedec2(cV,1,'sym5');%求取信號去噪的默認閾值、軟閾值、并保留低頻部分thr2,sorh2,keepapp2=ddencmp('den','wv',cV);%對三個方向高頻系數(shù)進行軟閾值處理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'); %對新的小波分解結構(nc，s)進行重構new_cVV=wkeep(new_cV,L,'c'); %提取重構的圖像%用小波函數(shù)sym5對cD進行1層小波分解C3,S3=wavedec2(cD,1,'sym5

17、');%求取信號去噪的默認閾值、軟閾值、并保留低頻部分thr3,sorh3,keepapp3=ddencmp('den','wv',cD);%對三個方向高頻系數(shù)進行軟閾值處理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、%對新的小波分解結構(nc，s)進行重構new_cDD=wkeep(new_cD,L,'c'); %提取重構的圖像%subplot(221);colormap(map);image(new_cHH);title('去噪后的水平分量');subplot(222);colormap(map);image(new_cVV);title('去噪后的垂直分量');subplot(223);colormap(map);image(new_cDD);title('去噪后的對角線分量'); % cA和去噪后的子圖（new_cH,new_cV,ne

19、w_cD）進行重構 % reconstruction為重構之后的圖形%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() 對輸入的r*c維矩陣 x 進行二維小波分解，輸出四個分解系數(shù)子矩陣LL,HL,LH,HH% 輸入?yún)?shù)：x 輸入矩陣，為r*c維矩陣。% 輸出參數(shù)：LL,HL,LH,HH 是分解系數(shù)矩陣的四個相等大小的子矩陣,大小均為 r/2 * c/2

20、維% LL：低頻部分分解系數(shù)； HL：垂直方向分解系數(shù)；% LH：水平方向分解系數(shù)； HH：對角線方向分解系數(shù)。 row,col=size(x); % 讀取輸入矩陣的大小for j=1:row % 首先對輸入矩陣的每一行序列進行一維離散小波分解 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 % 再對輸入矩陣的每一列序列進行一維離散小波分解 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('分解后的對角分量');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) 對輸入序列x進行一維離散小波分解，輸出分解序列cA,cD% 輸入?yún)?shù)：x輸入序列；% lpd低通濾波器；% hpd高通濾波器；% dim小波分解級數(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); % 通過下采樣求出本層分解后的細節(jié)部分系數(shù) n2=length(dnh); cA=dnl; % 下采樣后的平均部分系數(shù)進入下一層分解 cD=cD,dnh; % 將本層分解所得的細節(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ù)累加到對應序號的c矩陣數(shù)中endend4.4 采樣函數(shù)4.4.1 下采樣函數(shù)function y=downspl(x)% 函數(shù) Y=DOWMSPL(X) 對輸入序列進行下抽樣，輸出序列 Y。% 下抽樣是對輸入序列取其偶數(shù)位，舍棄奇數(shù)位。例如 x=x1,x2,x3,x4,x5，則 y=x2,x4. N=length(x); % 讀取輸入序列長度M=floor(N/2); % 輸出序列的長度是輸入序列長度的一半（帶小數(shù)時取整數(shù)部分）i=1:M;y(i)=x(2*i);4.4.2 上采樣函數(shù)function y=upspl(x)% 函數(shù) Y = UPSPL(X) 對輸入的一維序列x進行上抽樣，

26、即對序列x每個元素之間% 插零，例如 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ù)位按次序等于相應位置的輸入序列元素 if mod(i,2) y(i)=x(i+1)/2); else y(i)=0; endend4.5 重構函數(shù)4.5.1 二維重構函數(shù)function y=myidwt2(LL,HL,LH,HH,Lo_R,Hi_R,map)% 函數(shù) MYIDWT2() 對輸入的子矩陣序列進行逆小波變換，

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

28、1,cd1,Lo_R,Hi_R);% 重構序列 yt1(:,k)=tmp1; % 將重構序列存入待輸出矩陣 yt 的相應列，此時 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、細節(jié)系數(shù)序列cd2 cd2=tmp_mat(j,col/2+1:col); tmp2=myidwt(ca2,cd2,Lo_R,Hi_R); % 重構序列 yt(j,:)=tmp2; % 將重構序列存入

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