小波分析實(shí)驗(yàn)：二維離散小波變換

小波分析實(shí)驗(yàn)：實(shí)驗(yàn)2二維離散小波變換（Mallat快速算法）實(shí)驗(yàn)?zāi)康模涸诶斫怆x散小波變換原理和Mallat快速算法的基礎(chǔ)上，通過編程對(duì)圖像進(jìn)行二維離散小波變換，從而加深對(duì)二維小波分解和重構(gòu)的理性和感性認(rèn)識(shí)，并能提高編程能力，為今后的學(xué)習(xí)和工作奠定基礎(chǔ)。實(shí)驗(yàn)工具：計(jì)算機(jī)，matlab6.5分解算法:重構(gòu)算法:c^1=Vh{n-2k)c￡十Vg(刃-24)d￡+分解算法寫成矩陣的形式，（h,g的長(zhǎng)度為4）〃。[0]寸]W]

啪1]哉[。]&[]]h[2]國(guó)3]00「小[°]。[0]寸]W]

啪1]哉[。]&[]]h[2]國(guó)3]00「小[°]100用0]用1網(wǎng)[2]43]°勺一】［1］000??■.也2]段3]0000…覘]用1］」匠加-1］.我-1]?(1)3g[0]g[l]s[2]g[3]0...00-「勺.】［。］-00卻0]卻1]卻2]g[3]：.0*[1]..????????g[2]g[3]0000-g[0]g[l]_..q.LMT]_修1]?(2NJ?[0]h[A]如⑵h[3]0...00■■^T記尹二00勿⑼用1忡[2]勿[3]■-0????????????????。產(chǎn)勺[1].../[2]h\3]0000…勿[0]知1]_5-1]_■g[0]gmg⑵g[3]0???00一0.0.g].g[l]g⑵??g[3].???.0...卻2].g[3].0??00.0.???g[0].g[l]_4二

于是Mallat分解公式為矩陣變換…TOC\o"1-5"\h\zE=(3)vD產(chǎn)QD§_\⑷重構(gòu)算法寫成矩陣變換："h=g+S(5)3MATi?lJ-l峰值信噪比計(jì)算公式：哭…，專皿gMATi?lJ-l峰值信噪比計(jì)算公式：哭…，專皿g｛扁J，｛￡)分別表示原始圖像和重建圖像，本^^的一些好巧…(1)分■編程時(shí)用如下思想：（h,g的長(zhǎng)度為4）v京?！稠?]...京。〕頃1]...W-i]四0]

叩1]4皿-1]知0]切1]妃］徹3]00瑚］那］切2］風(fēng)3］.Ml].切2].雄]??00g⑼如］磯2]亦]00g[0]知］g[2]g[3].叩].g⑵.眼]??000.0???00將數(shù)據(jù)％_】和低通（高通）濾波器進(jìn)行添零到數(shù)據(jù)長(zhǎng)度m,再2抽樣，相加得到g)Stt袒，小波系數(shù)A后一半系數(shù)“內(nèi)插零“高通重構(gòu)「慮波器7重構(gòu)結(jié)果編程時(shí)用如下思想：Vii[O]抑]地］?_iDT0的..0000000????????000ii[O]抑]地］?_iDT0的..0000000????????000??心0..0心..?如]/(!]四../[0]_0叩]加）］

抑］00]00函)]000???g[3]00史］000…敗3]明「77+10-&-1[0]-撲]...—■"]-c"[l].十_^[0]-如11..孚伍-1］_寸_1以-1]_1]務(wù)尺度分解圖像50100150200250圖《1）三皴小波分解圖350100150200250E1C2[原始囹像和重構(gòu)圖像+J峰值信噪比：psnr=N52.8923dK附錄：二維小波分解函數(shù)%二維小波分解函數(shù)functionY=mallatdec2(X,wname,level)%輸入:X載入的二維圖像像數(shù)值；%level小波分解次(級(jí))數(shù)設(shè)定值(如果設(shè)定值超過最高可分解次數(shù)，按最高分解次數(shù)分解)%wname小波名字waveletname%輸出：Y多極小波分解后的小波系數(shù)矩陣[h,g]=wfilters(wname,'d');%h,g分別為低通和高通濾波器X=double(X);t=1;hh=size(X,2);whilet<=level%先進(jìn)行行小波變換forrow=1:hhY(row,1:hh)=mdec1(X(row,1:hh),h,g);end%再進(jìn)行列小波變換forcol=1:hhtemp=mdec1(Y(1:hh,col)',h,g);Y(1:hh,col)=temp';endt=t+1;hh=hh/2;X=Y;end%內(nèi)部子函數(shù)，對(duì)一行(row)矢量進(jìn)行一次小波變換，利用fft實(shí)現(xiàn)functiony=mdec1(x,h,g)%輸入：x行數(shù)組%h為低通濾波器%g為高通濾波器%輸出：y進(jìn)行一級(jí)小波分解后的系數(shù)lenx=size(x,2);rh=h(end：-1:1);rrh=[zeros(1,(lenx-lenh)),rh];rrh=circshift(rrh',1)';rg=g(end：-1:1);rrg=[zeros(1,(lenx-lenh)),rg];rrg=circshift(rrg',1)';r1=dyaddown(ifft(fft(x).*fft(rrh,lenx)),1);%usepara1r2=dyaddown(ifft(fft(x).*fft(rrg,lenx)),1);y=[r1,r2];二維小波重構(gòu)函數(shù)%二維小波重構(gòu)函數(shù)functionY=mallatrec2(X,wname,level)%輸入:X載入的小波系數(shù)矩陣；%level小波分解次(級(jí))數(shù)設(shè)定值(如果設(shè)定值超過最高可分解次數(shù)，按最高分解次數(shù)分解)%wname小波名字waveletname%輸出：Y重構(gòu)圖像矩陣[h,g]=wfilters(wname,'d');%h,g分別為重構(gòu)低通濾波器和重構(gòu)高通濾波器hz=size(X,2);h1=hz/(2“(levelT));whileh1<=hz%對(duì)列變換forcol=1:h1temp=mrec1(X(1:h1,col)',h,g)'X(1:h1,col)=temp;end%再對(duì)行變換forrow=1:h1temp=mrec1(X(row,1:h1),h,g);X(row,1:h1)=temp;endh1=h1*2;endY=X;%內(nèi)部子函數(shù)，對(duì)一行小波系數(shù)進(jìn)行重構(gòu)functiony=mrec1(x,h,g)%輸入：x行數(shù)組%h為低通濾波器%g為高通濾波器%輸出：y進(jìn)行一級(jí)小波重構(gòu)后值lenx=size(x,2);%內(nèi)插零usepara0r3=dyadup(x(1,1:lenx*0.5),0);r4=dyadup(x(1,(lenx*0.5+1):lenx),0);%usepara0y=ifft(fft(r3,lenx).*fft(h,lenx))+ifft(fft(r4,lenx).*fft(g,lenx));測(cè)試函數(shù)(主函數(shù))%測(cè)試函數(shù)(主函數(shù))clc;clear;X=imread('E:\Libin的文檔\Course\Course_wavelet、實(shí)驗(yàn)2要求\exp2\LENA.bmp');%%內(nèi)插零usepara0X=double(X);A=mallatdec2(X,'sym2',3);image(abs(A));colormap(gray(255));title('多尺度分解圖像')；Y=mallatrec2(A,'sym2',3);Y=real(Y);figure(2);subplot(1,2,1);image(X);colormap(gray(255));title('原始圖像')；subplot(1,2,2);image(Y);colormap(gray(255));title('重構(gòu)圖像')；csize=size(X);sr=csize(1);sc=csize(2);mse=sum(sum((Y-X)."2,1))/(sr*sc);psnr=10*log(255*255/mse)/log(10)小波分析實(shí)驗(yàn)：實(shí)驗(yàn)1連續(xù)小波變換實(shí)驗(yàn)?zāi)康模涸诶斫膺B續(xù)小波變換原理的基礎(chǔ)上，通過編程實(shí)現(xiàn)對(duì)一維信號(hào)進(jìn)行連續(xù)小波變換，(實(shí)驗(yàn)中采用的是墨西哥帽小波)，從而對(duì)連續(xù)小波變換增加了理性和感性的認(rèn)識(shí)，并能提高編程能力，為今后的學(xué)習(xí)和工作奠定基礎(chǔ)。實(shí)驗(yàn)工具：計(jì)算機(jī)，matlab6.5實(shí)驗(yàn)原理：~?維連續(xù)小波變換公式：“_Lg*t—bJF*E)二m戶JV(OW*(~y)UTOC\o"1-5"\h\z9/_A=1白尸Jr⑦w(—)出當(dāng)小波函教-03"編程時(shí)先離散化，在通過求和代替積分。。(燈1)由_]_})甲出展)=▽f/(r)x|f/|2V/(—)mF?以根據(jù)秒_1x,ArA?—b二“I2W(.^—k根據(jù)秒本實(shí)驗(yàn)采取的一些小技巧：襯（1）因?yàn)閷?shí)驗(yàn)中的數(shù)據(jù)通過load（data.mat'）（先將數(shù)據(jù)拷到work文件下）命令后，數(shù)據(jù)存在了dat矩陣中，所以就是霰（對(duì)；由于數(shù)捶樣間隔AZ為0.03（常量），所以而以把這個(gè)系數(shù)忽略，實(shí)際上本程序中-1k-b用的是wr（。0）=|s2￡尹修）w/（—L）（2）本實(shí)驗(yàn)中墨西哥帽函數(shù)為：呻維）={"—/）#:、0。otherwise（3）如果按照（*）式計(jì)算小波系數(shù)矩陣，則算法時(shí)間復(fù)雜度為口（疽），在通的電腦上要運(yùn)行5分鐘左右才能出正確結(jié)果。通過研究（**）,可以發(fā)對(duì)于同一尺度"，網(wǎng)*）與例上匚）有如下關(guān)系：網(wǎng)二二aaa二項(xiàng)開始與物H）相同，只是第一項(xiàng)不同，（在程序中有說明），所以a把算法時(shí)間復(fù)雜度降為Q（W）,程序執(zhí)行時(shí)間在5秒鐘內(nèi)。V（4）程序的正確性可由cwt來驗(yàn)證°Q（?）本實(shí)驗(yàn)還編了個(gè)test程序來說明存在哪些分量，即對(duì)正弦波函數(shù)樣后進(jìn)行連續(xù)小波變換，分析圖象結(jié)果。/

實(shí)驗(yàn)結(jié)果：3(1)實(shí)驗(yàn)數(shù)據(jù)圖像和小波系數(shù)圖像“

②test數(shù)據(jù)圖像和小波系數(shù)圖像"正弦波的小波系數(shù)圖204060801001201401S0180200程序附錄：⑴墨西哥帽小波函數(shù)，按照(***)式編程%mexh.mfunctionY=mexh(x)ifabs(x)<=8Y=exp(-x*x/2)*(1-x"2);elseY=0;End實(shí)驗(yàn)程序，按照(**)式編程，詳細(xì)過程請(qǐng)參考“本實(shí)驗(yàn)采取的一些小技巧”%clc;clear;load('data.mat');len=length(dat);lna=70;a=zeros(1,lna);wfab=zeros(lna,len);mexhab=zeros(1,len);%（尺度a）的長(zhǎng)度%小波系數(shù)矩陣%離散化小波系數(shù)矩陣fors=1:lna%s表示尺度fork=1:lenmexhab(k)=mexh(k/s);endfort=1:len%t表示位移%將積分用求和代%mexhab修改第一wfab(s,t)=(sum(mexhab.*dat))/sqrt(s);替mexhab=[mexh(-1*t/s),mexhab(1:len-1)];項(xiàng)并右移endendfigure(1);plot(dat);title('原始數(shù)據(jù)圖’)；figure(2);%小波系數(shù)譜%將積分用求和代%mexhab修改第一title('小波系數(shù)圖’)；%surf(wfab);%title('小波系數(shù)譜網(wǎng)格圖')；%pwfab二wfab.*wfab;%%瞬態(tài)功率譜%figure(3);%subplot(1,2,1);%surf(pwfab);%title('瞬態(tài)功率譜網(wǎng)格圖’)；%subplot(1,2,2);%contour(pwfab);%title('瞬態(tài)功率譜等值線')；test函數(shù)。%test函數(shù)clc;clear;fori=1:200dat(i)=sin(2*pi*i*0.05);%正弦波函數(shù)endlen=length(dat);lna=40;wfab=zeros(lna,len);mexhab=zeros(1,len);fors=1:lna%s表示尺度fork=1:lenmexhab(k)=mexh(k/s);

endfort=1:len%t表示位移wfab(s,t)=(sum(mexhab.*dat))/sqrt(s);%將積分用求和代替mexhab=[mexh(-1*t/s),mexhab(1:len-1)];%mexhab修改第一項(xiàng)并右移endendfigure(1);plot(dat);title('orignaldat');figure(2);%小波系數(shù)譜image(wfab);colormap(pink(128));title('正弦波的小波系數(shù)圖’)；用fft實(shí)現(xiàn)cwt%按照?qǐng)A周卷積定理，原周卷積和線性卷積的關(guān)系L>=M+N-1%按照?qǐng)A周卷積的定義，相關(guān)和線性卷積的關(guān)系(原始算法和線性卷積的關(guān)系)%注意畫圖理解clc;clear;t1二cputime;load('data.mat');len=length(dat);lna=70;%a(尺度)的長(zhǎng)度a=zeros(1,lna);%a表示尺度b=zeros(1,len);%b表示位移wfab=zeros(lna,len);%小波系數(shù)矩陣mexhab二zeros(1,2*len-1);data=[zeros(1,len-1),dat];Ydata二fft(data,4*len);fors=1:lnafork=1:2*len-1mexhab(k)=mexh((k-len)/s);enda=zeros(1,lna);%a表示尺度temp二ifft(Ydata.*fft(mexhab,4*len),4*len);wfab(s,:)=real(temp(2*len-1:3*len-2))/sqrt(s);%為什么要取實(shí)部而不是取模，我也不是很清楚