小波实验报告双树复小波变换

1、相位：Zdc (j, n) = arctan2.代码实现。我采用 Matlab 函数编程实现。具体程序见shift_test_2D.m , drawcirc.m , setfig.m ,dtwavexfm2.m , dtwaveifm2.m,waveifm2.m,wavexfm2.mSkelMap.m。设- h , 和一；h¥h分别是双正交对偶尺度函数与对偶小波， h1 n 是相应的低通滤波器和高通滤波器，即它们满足_实部：h(t) 处nFt)八2、n看t)二.2"nho(n) h(2t -n)h(n) h(2t -n)K(n )"h(2t-n)*h(t)二 2、

2、hi(n)%h(2t -n)n双树复小波变换可以通过离散小波变换 产生虚部。虚部:DWTg(t)"办n*g(t)二云n冷(t) = 3ng(t)=云n实现：一个 DWT产生实部，另一个go(n) g(2t- n)gi(n) g(2t- n)g°(n)"g(2t- n)gi(n)"g(2t- n)一、题目：双树复小波变换二、目的：双树复小波和实小波变换的比较三、算法及其实现： 提取阶梯型边界点1.算法。' c(t)八 h(t)L g(t)幅值:|dc(j,n)|二dh(j, n)2 dg(j, n) d g(j，n)<dh(j,n)四、实现工

3、具：Matlab五、程序代码：(1) shift_test_2D.m :% shift_test_2D.m% M-file to perform a 4-level wavelet transform on a circle using Q-shift% dual wavelet tree and DWT, and to compare shift in varia nee properties.% Nick Ki ngsbury, Cambridge Un iversity. May 2002.clear all close all% Draw a circular disc.x = rou

4、nd(drawcirc(64,1,0,0,256) - 0.5) * 200); setfig(1);colormap(gray(256)image(mi n(max(x+128,1),256);set(gca,'positio n',0.1 0.25 .25 .5);axis('off);axis('image');% draw(xx);title(' In put (256 x 256)','Fo ntSize',14); draw now% Do 4 levels of CWT.YI,Yh = dtwavexfm2(

5、x,4,' near_sym_b','qshift_b');% Loop to recon struct output from coefs at each level in turn.% Starts with the fin est level.titl = '1st''2 nd''3rd''4th''Low'yy = zeros(size(x) .* 2 3);yt1 = 1:size(x,1); yt2 = 1:size(x,2);for mlev = 1:5,mask =

6、zeros(6,5);mask(:,mlev) = 1;z = dtwaveifm2(YI*mask(1,5),Yh,' near_sym_b','qshift_b',mask); figure;draw(z);draw nowyy(yt1,yt2)=乙yt2 = yt2 + size(x,2)/2;end% disp('Press a key .')% pause% Now do same with DWT.% Do 4 levels of Real DWT usi ng 'antonini' (9,7)-tap filters

7、.Yl,Yh = wavexfm2(x,4,'a nto nini');yt1 = 1:size(x,1) + size(x,1); yt2 = 1:size(x,2);for mlev = 1:5,mask = zeros(3,5);mask(:,mlev) = 1;z = waveifm2(Yl*mask(1,5),Yh,'a nto nini ',mask);figure;draw(z);draw nowyy(yt1,yt2)=乙yt2 = yt2 + size(x,2)/2;endfigure;setfig(gcf);colormap(gray(256)

8、image(mi n(max(yy+128,1),256);set(gca,'positio n',0.1 0.1 .8 .8);axis('off);axis('image');hold onplot(128* 1;1*1:4 0;6+1,128*0;4*1 1 1 1 2;2+1,'-k');hold offtitle('Comp onents of recon structed "disc" images','F on tSize',14); text(-0.01*size(yy,

9、2),0.25*size(yy,1),'DT CWT','horiz','r'); text(0.02*size(yy,2),1.02*size(yy,1),'wavelets:','horiz','r','vert','t'); text(-0.01*size(yy,2),0.75*size(yy,1),'DWT','horiz','r');endfor k=1:4, text(k*128-63,size(yy,1)*

10、1.02,spri ntf('level %d',k),'Fo ntSize',14,'horiz','c','vert','t');text(5*128+1,size(yy,1)*1.02,'level 4 scali ng fn.','Fo ntSize',14,'horiz','c','vert','t');draw now% print -deps circrecq.epsdisp('Pr

11、ess a key to see perfect recon structi on property .')pause% Accumulate the images from lowba nd upwards to show perfect recon structi on.sy = size(x,2)/2;for mlev = 4:-1:1,yt2 = 1:sy + (mlev-1)*sy;yy(:,yt2) = yy(:,yt2) + yy(:,yt2+sy);endfigure;setfig(gcf);colormap(gray(256)image(mi n(max(yy+128

12、,1),256);set(gca,'positio n',0.1 0.1 .8 .8);axis('off);axis('image');title('Accumulated recon structio ns from each level of DT CWT ','Fo ntSize',14); text(size(yy,2)*0.5,size(yy,1)*1.02,'Accumulatedreconstructions from each level of DWT','Fo ntSize

13、9;,14,'hor','c','vert','t');draw nowreturn(2) function p = drawcirc(r,w,dx,dy,N)% function p = drawcirc(r,w,dx,dy,N)% Gen erate an image of size N*N pels, containing a circle% radius r pels and cen tred at dx,dy relative% to the centre of the image. The edge of the ci

14、rcle is a cos ine shaped% edge of width w (from 10 to 90% poi nts).x = on es(N,1) * (1:N - (N+1)/2 - dx)/r);y = (1:N' - (N+1)/2 - dy)/r) * on es(1,N);p = 0.5 + 0.5 * sin (mi n( max(exp(-0.5 * (x + y ) - exp(-0.5)*(r*3/w), -pi/2), pi/2);return(3) function Z = dtwaveifm2(YI,Yh,biort,qshift,gain_ma

15、sk);% Function to perform an n-level dual-tree complex wavelet (DTCWT)% 2-D recon structi on.% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Yl -> The real lowpass image from the fin al level% Yh -> A cell array containing the 6 complex highpass subimages for each level.% biort ->'an t

16、o nini'=> Antonini 9,7 tap filters.%'legall' => LeGall 5,3 tap filters.%'n ear_sym_a' => Near-Symmetric 5,7 tap filters.%'n ear_sym_b' => Near-Symmetric 13,19 tap filters.% qshift -> 'qshift_06' => Quarter Sample Shift Orthogo nal (Q-Shift) 10,10

17、 tap filters,%(only 6,6 non-zero taps).'qshift_a' => Q-shift 10,10 tap filters,(with 10,10 non-zero taps, un like qshift_06).'qshift_b' => Q-Shift 14,14 tap filters.'qshift_c' => Q-Shift 16,16 tap filters.'qshift_d' => Q-Shift 18,18 tap filters.%gain _mask

18、 -> Gain to be applied to each subba nd.gain _mask(d,l) is gain for subba nd with direct ion d at level l.If gain _mask(d,l) = 0, no computatio n is performed for band (d,l).Default gain _mask = on es(6,le ngth(Yh).% Z -> Rec on structed real image matrix% For example: Z = dtwaveifm2(YI,Yh,

19、9; near_sym_b','qshift_b');% performs a 3-level recon structio n from Yl,Yh using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.% Nick Kin gsbury and Cian Shaffrey% Cambridge Un iversity. May 2002a = len gth(Yh); % No of levels.if nargin < 5, gai

20、n _mask = on es(6,a); end % Default gain _mask.if isstr(biort) & isstr(qshift)%Check if the in puts are stri ngsbiort_exist = exist(biort '.mat');qshift_exist = exist(qshift '.mat');if biort_exist = 2 & qshift_exist = 2; %Check to see if the in puts exist as .mat files load (

21、biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTWAVEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTWAVEIFM2.');endcurre nt_level = a;Z = Yl;while curr

22、e nt_level >= 2; ; %this en sures that for level -1 we n ever do the followi nglh = c2q(Yhcurre nt_level(:，:,1 6),gain_mask(1 6,curre nt_level);hl = c2q(Yhcurre nt_level(:,:,3 4),gain_mask(3 4,curre nt_level);hh = c2q(Yhcurre nt_level(:,:,2 5),gain_mask(2 5,curre nt_level);% Do even Qshift filter

23、s on colu mns.yl = colifilt(Z,gOb,gOa) + colifilt(lh,g1b,g1a);y2 = colifilt(hl,gOb,gOa) + colifilt(hh,g1b,g1a);% Do even Qshift filters on rows.Z = (colifilt(y1.',gOb,gOa) + colifilt(y2.',g1b,g1a).'% Check size of Z and crop as requiredrow_size col_size = size(Z);S = 2*size(Yhcurre nt_le

24、vel-1);if row_size = S(1)%check to see if this result n eeds to be cropped for the rowsZ = Z(2:row_size-1,:);endif col_size = S(2)%check to see if this result n eeds to be cropped for the colsZ = Z(:,2:col_size-1);endif an y(size(Z) = S(1:2), error('Sizes of subbands are not valid for DTW AVEIFM

25、2');endcurre nt_level = curre nt_level - 1;endif curre nt_level = 1;lh = c2q(Yhcurre nt_level(:,:,1 6),gain_mask(1 6,curre nt_level);hl = c2q(Yhcurre nt_level(:,:,3 4),gain_mask(3 4,curre nt_level);hh = c2q(Yhcurre nt_level(:,:,2 5),gain_mask(2 5,curre nt_level);% Do odd top-level filters on col

26、u mns.y1 = colfilter(Z,g0o) + colfilter(lh,g1o);y2 = colfilter(hl,g0o) + colfilter(hh,g1o);% Do odd top-level filters on rows.Z = (colfilter(y1.',g0o) + colfilter(y2.',g1o).'endreturn*INTERNAL FUNCTION*function x = c2q(w,ga in)% fun ctio n z = c2q(w,ga in)% Scale by gain and convert from

27、 complex w(:,:,1:2) to real quad-nu mbers in 乙% Arrange pixels from the real and imag parts of the 2 subba nds% into 4 separate subimages .%ABReIm of w(:,:,1)%1 1%1 1%CDReIm of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2);if any (w(:) & any (ga in)sc = sqrt(0.5) * gain;P = w(:,:,1)*sc(1) + w(:,:,2)*

28、sc(2);Q = w(:,:,1)*sc(1) - w(:,:,2)*sc(2);t1 = 1:2:size(x,1);t2 = 1:2:size(x,2);% Recover each of the 4 corners of the quads.x(t1,t2)= real(P); % a = (A+C)*sc;x(t1,t2+1)= imag(P); % b = (B+D)*sc;x(t1+1,t2)= imag(Q); % c = (B-D)*sc;x(t1+1,t2+1) = -real(Q); % d = (C-A)*sc;endreturn(3) function Yl,Yh,Y

29、scale = dtwavexfm2(X,nlevels,biort,qshift);% Function to perform a n-level DTCWT-2D decomposti on on a 2D matrix X% Yl,Yh,Yscale = dtwavexfm2(X, nlevels,biort,qshift);% X -> 2D real matrix/Image%nl evels -> No. of levels of wavelet decompositi on% biort ->'an to nini'=> Antonini

30、9,7 tap filters.%'legall' => LeGall 5,3 tap filters.%'n ear_sym_a' => Near-Symmetric 5,7 tap filters.%'n ear_sym_b' => Near-Symmetric 13,19 tap filters.% qshift -> 'qshift_06' => Quarter Sample Shift Orthogo nal (Q-Shift) 10,10 tap filters,%(only 6,6 no

31、n-zero taps).%'qshift_a' => Q-shift 10,10 tap filters,%(with 10,10 non-zero taps, uni ike qshift_06).%'qshift_b' => Q-Shift14,14 tap filters.%'qshift_c' => Q-Shift16,16 tap filters.%'qshift_d' => Q-Shift18,18 tap filters.% Yl -> The real lowpass image f

32、rom the fin al level% Yh -> A cell array containing the 6 complex highpass subimages for each level.%Yscale -> This is an OPTIONAL output argume nt, that is a cell array containing%real lowpass coefficie nts for every scale.% Example: YI,Yh = dtwavexfm2(X,3,' near_sym_b','qshift_b&

33、#39;);% performs a 3-level tran sform on the real image X using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.% Nick Kin gsbury and Cian Shaffrey% Cambridge Uni versity, Sept 2001if isstr(biort) & isstr(qshift)%Check if the in puts are stri ngsbiort_exist =

34、 exist(biort '.mat');qshift_exist = exist(qshift '.mat');if biort_exist = 2 & qshift_exist = 2;%Check to see if the in puts existas .mat filesload (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTWAVEXFM2 for d

35、etails.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTWAVEXFM2.');end orgi nal_size = size(X);if ndims(X) >= 3;error(sprintf('The entered image is %dx%dx%d, please enter each image sliceseparately.',org in al_size(1),org in al

36、_size(2),org in al_size(3); end% The n ext few lines of code check to see if the image is odd in size, if so an extra .% row/colu mn will be added to the bottom/right of the imagein itial_row_exte nd = 0;%in itialisein itial_col_exte nd = 0;if any(rem(orgin al_size(1),2), %if sx(1) is not divisable

37、by 2 the n we n eed to exte nd X by addi ng a row at the bottomX = X; X(e nd,:); in itial_row_exte nd = 1;endif any(rem(orgin al_size(2),2), addi ng a col to the leftX = X X(:,e nd);in itial_col_exte nd = 1;endexte nded_size = size(X);%Any further exte nsion will be done in due course.%if sx(2) is n

38、ot divisable by 2 then we need to extend X by%Any further exte nsion will be done in due course.if nl evels = 0, retur n; end%in itialiseYh=cell( nlevels,1);if n argout = 3Yscale=cell( nl evels,1);comp onent.end%this is only required if the user specifies a third outputS =; sx = size(X); if nl evels

39、 >= 1,% Do odd top-level filters on cols.Lo = colfilter(X,h0o).'Hi = colfilter(X,h1o).'% Do odd top-level filters on rows.LoLo = colfilter(Lo,hOo).'Yh1 = zeros(size(LoLo)/2 6);Yh1(:,:,1 6) = q2c(colfilter(Hi,h0o).'); Yh1(:，:,3 4) = q2c(colfilter(Lo,h1o).'); Yh1(:，:,2 5) = q2c(

40、colfilter(Hi,h1o).'); S = size(LoLo) ;S;% LoLo% Horiz on tal pair% Vertical pair% Diag onal pairif n argout = 3Yscale1 = LoLo;endendif nl evels >= 2;for level = 2: nlevels;row_size col_size = size(LoLo);if any(rem (row_size,4),% Exte nd by 2 rows if no. of rows of LoLo are divisableby 4;LoLo

41、= LoLo(1,:); LoLo; LoLo(e nd,:);endif any( rem(col_size,4),% Exte nd by 2 cols if no. of cols of LoLo are divisable by4;LoLo = LoLo(:,1) LoLo LoLo(:,e nd);end% Do eve n Qshift filters on rows.Lo = coldfilt(LoLo,hOb,hOa).'Hi = coldfilt(LoLo,h1b,h1a).'% Do eve n Qshift filters on colu mns.LoLo

42、 = coldfilt(Lo,hOb,hOa).' %LoLoYhlevel = zeros(size(LoLo)/2 6);Yhlevel(:，:,1 6) = q2c(coldfilt(Hi,h0b,h0a).');%Horiz ontalYhlevel(:，:,3 4) = q2c(coldfilt(Lo,h1b,h1a).');%VerticalYhlevel(:，:,2 5) = q2c(coldfilt(Hi,h1b,h1a).');%Diago nalS = size(LoLo) ;S;if n argout = 3Yscalelevel = Lo

43、Lo;endendendYl = LoLo;if in itial_row_exte nd = 1 & in itial_col_exte nd = 1;warnin g(spri ntf(' rr The image en tered is now a %dx%d NOT a %dx%d r The bottom row and rightmost colu mn have bee n duplicated, prior to decompositi on. rr ',.exte nded_size(1),exte nded_size(2),org in al_siz

44、e(1),org in al_size(2);endif in itial_row_exte nd = 1 ;warnin g(spri ntf(' rr The image en tered is now a %dx%d NOT a %dx%d r Row nu mber %d has bee n duplicated, and added to the bottom of the image, prior to decompositi on. rr',.exte nded_size(1),exte nded_size(2),orgi nal_size(1),orgi nal

45、_size(2),orgi nal_size(1);end if in itial_col_exte nd = 1;warnin g(spri ntf(' rr The image en tered is now a %dx%d NOT a %dx%d r Col nu mber %d has bee n duplicated, and added to the right of the image, prior to decompositi on. rr',.exte nded_size(1),exte nded_size(2),org in al_size(1),org i

46、n al_size(2),org in al_size(2); end return%=*INTERNAL FUNCTION*%= fun ctio n z = q2c(y)% fun ctio n z = q2c(y)% Convert from quads in y to complex nu mbers in 乙sy = size(y);t1 = 1:2:sy(1); t2 = 1:2:sy(2);j2 = sqrt(0.5 -0.5);% Arrange pixels from the corners of the quads into % 2 subimages of alterna

47、te real and imag pixels.% ab% II% II%c-d% Comb ine (a,b) and (d,c) to form two complex subimages.p = y(t1,t2)*j2(1) + y(t1,t2+1)*j2 (2);% p = (a + jb) / sqrt(2)q = y(t1+1,t2+1)*j2(1) - y(t1+1, t2)*j2(2); % q = (d - jc) / sqrt(2)% Form the 2 subba nds in z.z = cat(3,p-q,p+q);return(4) function Z = dt

48、waveifm2(Yl,Yh,biort,qshift,gain_mask);% Function to perform an n-level dual-tree complex wavelet (DTCWT)% 2-D recon structi on.% Z = dtwaveifm2(YI,Yh,biort,qshift,gain_mask);%Yl -> The real lowpass image from the fin al levelYh -> A cell array containing the 6 complex highpass subimages for e

49、ach level.%biort ->'an to nini'=> Antonini 9,7 tap filters.'legall' => LeGall 5,3 tap filters.'n ear_sym_a' => Near-Symmetric 5,7 tap filters.'n ear_sym_b' => Near-Symmetric 13,19 tap filters.qshift -> 'qshift_06' => Quarter Sample Shift O

50、rthogo nal (Q-Shift) 10,10 tap filters,(only 6,6 non-zero taps).'qshift_a' => Q-shift 10,10 tap filters,(with 10,10 non-zero taps, uni ike qshift_06).'qshift_b' => Q-Shift 14,14 tap filters.'qshift_c' => Q-Shift 16,16 tap filters.'qshift_d' => Q-Shift 18,1

51、8 tap filters.gain _mask -> Gain to be applied to each subba nd.gain _mask(d,l) is gain for subba nd with direct ion d at level l.If gain _mask(d,l) = 0, no computati on is performed for band (d,l).Default gain _mask = on es(6,le ngth(Yh).Z -> Recon structed real image matrix % For example: Z

52、= dtwaveifm2(Yl,Yh,' near_sym_b','qshift_b');% performs a 3-level recon structio n from Yl,Yh using the 13,19-tap filters% for level 1 and the Q-shift 14-tap filters for levels >= 2.% Nick Kin gsbury and Cian Shaffrey% Cambridge Un iversity. May 2002 a = len gth(Yh); % No of level

53、s.if nargin < 5, gain _mask = on es(6,a); end% Default gain _mask.if isstr(biort) & isstr(qshift)%Check if the in puts are stri ngsbiort_exist = exist(biort '.mat');qshift_exist = exist(qshift '.mat');if biort_exist = 2 & qshift_exist = 2; %Check to see if the in puts exis

54、t as .mat files load (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTWAVEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTWAVEIFM2.');endcurre nt_level

55、 = a;Z = Yl;while curre nt_level >= 2; ;%this en sures that for level -1 we n ever do the followi nglh = c2q(Yhcurre nt_level(:，:,1 6),gain_mask(1 6,curre nt_level);hl = c2q(Yhcurre nt_level(:,:,3 4),gain_mask(3 4,curre nt_level);hh = c2q(Yhcurre nt_level(:,:,2 5),gain_mask(2 5,curre nt_level);%

56、Do even Qshift filters on colu mns.y1 = colifilt(Z,gOb,gOa) + colifilt(lh,g1b,g1a);y2 = colifilt(hl,gOb,gOa) + colifilt(hh,g1b,g1a);% Do even Qshift filters on rows.Z = (colifilt(y1.',g0b,g0a) + colifilt(y2.',g1b,g1a).'% Check size of Z and crop as requiredrow_size col_size = size(Z);S =

57、 2*size(Yhcurre nt_level-1);if row_size = S(1)%check to see if this result n eeds to be cropped for the rowsZ = Z(2:row_size-1,:);endif col_size = S(2)%check to see if this result n eeds to be cropped for the colsZ = Z(:,2:col_size-1);endif an y(size(Z) = S(1:2),error('Sizes of subbands are not

58、valid for DTW AVEIFM2');endcurre nt_level = curre nt_level - 1;endif curre nt_level = 1;lh = c2q(Yhcurre nt_level(:,:,1 6),gain_mask(1 6,curre nt_level);hl = c2q(Yhcurre nt_level(:,:,3 4),gain_mask(3 4,curre nt_level);hh = c2q(Yhcurre nt_level(:,:,2 5),gain_mask(2 5,curre nt_level);% Do odd top-level filters