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

1、小波分析實(shí)驗(yàn)：實(shí)驗(yàn)2 二維離散小波變換(Mallat快速算法)實(shí)驗(yàn)?zāi)康模涸诶斫怆x散小波變換原理和Mallat快速算法的基礎(chǔ)上，通過(guò)編程對(duì)圖像進(jìn)行二維離散小波變換，從而加深對(duì)二維小波分解和重構(gòu)的理性和感性認(rèn)識(shí)，并能提高編程能力，為今后的學(xué)習(xí)和工作奠定基礎(chǔ)。實(shí)驗(yàn)工具：計(jì)算機(jī)，matlab6.5 附錄：（1）二維小波分解函數(shù)%二維小波分解函數(shù) function Y=mallatdec2(X,wname,level)%輸入:X      載入的二維圖像像數(shù)值; %     level  小波分

2、解次(級(jí))數(shù)設(shè)定值(如果設(shè)定值超過(guò)最高可分解次數(shù),按最高分解次數(shù)分解)%      wname  小波名字wavelet name%輸出：Y     多極小波分解后的小波系數(shù)矩陣 h,g=wfilters(wname,'d');  %h,g分別為低通和高通濾波器X=double(X);t=1;hh=size(X,2); while t<=level      %先進(jìn)行行小波變換  &#

3、160; for  row=1:hh        Y(row,1:hh)=mdec1(X(row,1:hh),h,g) ;    end    %再進(jìn)行列小波變換    for col=1:hh        temp=mdec1( Y(1:hh,col)',h,g);      

4、0; 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)function y=mdec1(x,h,g)%輸入：x 行數(shù)組%     h為低通濾波器%     g為高通濾波器%輸出： y 進(jìn)行一級(jí)小波分解后的系數(shù)lenx=size(x,2);lenh=size(h,2); rh=h(end:-1:1);rrh=zeros(1,(lenx-lenh),rh

5、;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);  %use para 1r2=dyaddown(ifft(fft(x).*fft(rrg,lenx),1);y=r1,r2; （2）二維小波重構(gòu)函數(shù)%二維小波重構(gòu)函數(shù)function Y=mallatrec2(X,wname,level)%輸入:X  &#

6、160;   載入的小波系數(shù)矩陣; %     level  小波分解次(級(jí))數(shù)設(shè)定值(如果設(shè)定值超過(guò)最高可分解次數(shù),按最高分解次數(shù)分解)%      wname  小波名字wavelet name%輸出：Y     重構(gòu)圖像矩陣 h,g=wfilters(wname,'d');  %h,g分別為重構(gòu)低通濾波器和重構(gòu)高通濾波器 hz=size(X,2);h1=hz/(2(level

7、-1); while h1<=hz     % 對(duì)列變換    for col=1:h1        temp=mrec1(X(1:h1,col)',h,g)'        X(1:h1,col)=temp;    end    %再對(duì)行變換    for row=1:h1 

8、;       temp=mrec1(X(row,1:h1),h,g);         X(row,1:h1)=temp;     end     h1=h1*2;     endY=X;   %內(nèi)部子函數(shù)，對(duì)一行小波系數(shù)進(jìn)行重構(gòu)function y=mrec1(x,h,g)%輸入：x 行數(shù)組% 

9、    h為低通濾波器%     g為高通濾波器%輸出： y 進(jìn)行一級(jí)小波重構(gòu)后值lenx=size(x,2);  r3=dyadup(x(1,1:lenx*0.5),0);   %內(nèi)插零use para 0r4=dyadup(x(1,(lenx*0.5+1):lenx),0);   %use para 0y=ifft(fft(r3,lenx).*fft(h,lenx)+ ifft(fft(r4,lenx).*fft(g,lenx);  

10、60;（3）測(cè)試函數(shù)（主函數(shù)）%測(cè)試函數(shù)（主函數(shù)）clc;clear;X=imread('E:Libin的文檔CourseCourse_wavelet實(shí)驗(yàn)2要求exp2LENA.bmp');%路徑 X=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);

11、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ǔ)上，通過(guò)編程實(shí)現(xiàn)對(duì)一維信號(hào)進(jìn)行連續(xù)小波變換，（實(shí)驗(yàn)中采用的是墨西哥帽小波），從而對(duì)連續(xù)小波變換增加了理性和感性的

12、認(rèn)識(shí)，并能提高編程能力，為今后的學(xué)習(xí)和工作奠定基礎(chǔ)。實(shí)驗(yàn)工具：計(jì)算機(jī)，matlab6.5  程序附錄：(1) 墨西哥帽小波函數(shù),按照(*)式編程% mexh.mfunction Y=mexh(x)if abs(x)<=8Y=exp(-x*x/2)*(1-x2);else     Y=0;End (2) 實(shí)驗(yàn)程序，按照(*)式編程，詳細(xì)過(guò)程請(qǐng)參考“本實(shí)驗(yàn)采取的一些小技巧”% clc;clear;load('data.mat');len=length(dat);lna=70;    &#

13、160;      % （尺度a）的長(zhǎng)度a=zeros(1,lna);wfab=zeros(lna,len);   %小波系數(shù)矩陣mexhab=zeros(1,len);   % 離散化小波系數(shù)矩陣 for s=1:lna           %s 表示尺度      for k=1:len     

14、0;  mexhab(k)=mexh(k/s);    end       for t=1:len    % t 表示位移      wfab(s,t)=(sum(mexhab.*dat)/sqrt(s);   %將積分用求和代替      mexhab=mexh(-1*t/s),mexhab(1:len-1);  %mexhab修改第一項(xiàng)并右移&

15、#160;   endend figure(1);plot(dat);title('原始數(shù)據(jù)圖');figure(2);  %小波系數(shù)譜image(wfab);colormap(pink(128);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,

16、2);%contour(pwfab);%title('瞬態(tài)功率譜等值線'); (3)test函數(shù)。%test 函數(shù)clc;clear;for i=1:200    dat(i)=sin(2*pi*i*0.05);     %正弦波函數(shù)endlen=length(dat);lna=40;wfab=zeros(lna,len);mexhab=zeros(1,len);for s=1:lna         

17、0; %s 表示尺度      for k=1:len        mexhab(k)=mexh(k/s);    end       for t=1:len               % t 表示位移      w

18、fab(s,t)=(sum(mexhab.*dat)/sqrt(s);   %將積分用求和代替      mexhab=mexh(-1*t/s),mexhab(1:len-1);  %mexhab修改第一項(xiàng)并右移    endendfigure(1);plot(dat);title('orignal dat');figure(2);  %小波系數(shù)譜image(wfab);colormap(pink(128);title('正弦波的小波系數(shù)圖')

19、;（4）用fft實(shí)現(xiàn)cwt%按照?qǐng)A周卷積定理，原周卷積和線性卷積的關(guān)系L>=M+N-1%按照?qǐng)A周卷積的定義，相關(guān)和線性卷積的關(guān)系（原始算法和線性卷積的關(guān)系）%注意畫(huà)圖理解clc;clear;t1=cputime; load('data.mat');len=length(dat);lna=70;           % a（尺度）的長(zhǎng)度a=zeros(1,lna);    % a 表示尺度b=zeros(1,len); 

20、60;   % b 表示位移wfab=zeros(lna,len);   %小波系數(shù)矩陣mexhab=zeros(1,2*len-1);   data=zeros(1,len-1),dat;Ydata=fft( data ,4*len);for s=1:lna              for k=1:2*len-1       mexhab(k)=mexh

21、(k-len)/s);        end      temp=ifft( Ydata.*fft( mexhab,4*len ) ,4*len);   wfab(s,:)=real(temp(2*len-1:3*len-2)/sqrt(s); %為什么要取實(shí)部而不是取模，我也不是很清楚，可是有種感覺(jué)endfigure(1);plot(dat);title('原始數(shù)據(jù)圖');figure(2);  %小波系數(shù)譜image(wfab);colormap(pink(128);title('小波系數(shù)譜 ');cputime-t1   4