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

1、一、题目：双树复小波变换二、目旳：双树复小波和实小波变换旳比较三、算法及其实现：提取阶梯型边界点算法。幅值:相位：代码实现。我采用Matlab函数编程实现。具体程序见shift_test_2D.m，drawcirc.m，setfig.m，dtwavexfm2.m，dtwaveifm2.m,waveifm2.m,wavexfm2.mSkelMap.m。设 和 分别是双正交对偶尺度函数与对偶小波， ， ， 和 是相应旳低通滤波器和高通滤波器，即它们满足实部： 虚部：双树复小波变换可以通过离散小波变换DWT实现：一种DWT产生实部，另一种产生虚部。实现工具：Matlab程序代码： （1）shift_

2、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 invariance properties.% Nick Kingsbury, Cambridge University, May .clear allclose all% Draw a circular disc.x = round(drawcirc(64,1,0,0,256) - 0.5) *

3、 200);setfig(1); colormap(gray(256)image(min(max(x+128,1),256);set(gca,position,0.1 0.25 .25 .5);axis(off);axis(image);% draw(xx); title(Input (256 x 256),FontSize,14); drawnow% Do 4 levels of CWT.Yl,Yh = dtwavexfm2(x,4,near_sym_b,qshift_b);% Loop to reconstruct output from coefs at each level in tu

4、rn.% Starts with the finest level.titl = 1st;2nd;3rd;4th;Low;yy = zeros(size(x) .* 2 3);yt1 = 1:size(x,1); yt2 = 1:size(x,2);for mlev = 1:5, mask = zeros(6,5); mask(:,mlev) = 1; z = dtwaveifm2(Yl*mask(1,5),Yh,near_sym_b,qshift_b,mask); figure;draw(z);drawnow yy(yt1,yt2) = z; yt2 = yt2 + size(x,2)/2;

5、end% disp(Press a key .)% pause% Now do same with DWT.% Do 4 levels of Real DWT using antonini (9,7)-tap filters.Yl,Yh = wavexfm2(x,4,antonini);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,antonini,mask); figure;dr

6、aw(z);drawnow yy(yt1,yt2) = z; yt2 = yt2 + size(x,2)/2;endfigure;setfig(gcf); colormap(gray(256)image(min(max(yy+128,1),256);set(gca,position,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(Components of reconstructed disc images,FontSize,14

7、);text(-0.01*size(yy,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);for k=1:4, text(k*128-63,size(yy,1)*1.02,sprintf(level %d,k),FontSize,14,horiz,c,vert,t); endtext(5*128+1,size(yy,1)*1.02,level 4

8、scaling fn.,FontSize,14,horiz,c,vert,t);drawnow% print -deps circrecq.epsdisp(Press a key to see perfect reconstruction property .)pause% Accumulate the images from lowband upwards to show perfect reconstruction.sy = size(x,2)/2;for mlev = 4:-1:1, yt2 = 1:sy + (mlev-1)*sy; yy(:,yt2) = yy(:,yt2) + yy

9、(:,yt2+sy);endfigure;setfig(gcf); colormap(gray(256)image(min(max(yy+128,1),256);set(gca,position,0.1 0.1 .8 .8);axis(off);axis(image);title(Accumulated reconstructions from each level of DT CWT ,FontSize,14);text(size(yy,2)*0.5,size(yy,1)*1.02,Accumulated reconstructions from each level of DWT ,Fon

10、tSize,14,hor,c,vert,t);drawnowreturn（2）function p = drawcirc(r,w,dx,dy,N)% function p = drawcirc(r,w,dx,dy,N)% Generate an image of size N*N pels, containing a circle% radius r pels and centred at dx,dy relative% to the centre of the image. The edge of the circle is a cosine shaped% edge of width w

11、(from 10 to 90% points).x = ones(N,1) * (1:N - (N+1)/2 - dx)/r);y = (1:N - (N+1)/2 - dy)/r) * ones(1,N);p = 0.5 + 0.5 * sin(min(max(exp(-0.5 * (x + y ) - exp(-0.5)*(r*3/w), -pi/2), pi/2);return（3）function Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Function to perform an n-level dual-tree complex

12、 wavelet (DTCWT)% 2-D reconstruction.% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% % Yl - The real lowpass image from the final level% Yh - A cell array containing the 6 complex highpass subimages for each level.% biort - antonini = Antonini 9,7 tap filters.% legall = LeGall 5,3 tap filters.% nea

13、r_sym_a = Near-Symmetric 5,7 tap filters.% near_sym_b = Near-Symmetric 13,19 tap filters.% qshift - qshift_06 = Quarter Sample Shift Orthogonal (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, unlike qshift_06).% qshift_b = Q

14、-Shift 14,14 tap filters.% qshift_c = Q-Shift 16,16 tap filters.% qshift_d = Q-Shift 18,18 tap filters.% gain_mask - Gain to be applied to each subband. % gain_mask(d,l) is gain for subband with direction d at level l.% If gain_mask(d,l) = 0, no computation is performed for band (d,l).% Default gain

15、_mask = ones(6,length(Yh).% Z - Reconstructed real image matrix% % For example: Z = dtwaveifm2(Yl,Yh,near_sym_b,qshift_b);% performs a 3-level reconstruction from Yl,Yh using the 13,19-tap filters % for level 1 and the Q-shift 14-tap filters for levels = 2.% Nick Kingsbury and Cian Shaffrey% Cambrid

16、ge University, May a = length(Yh); % No of levels.if nargin = 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yhcurrent_level(:,:,1 6),gain_mask(1 6,current_level); hl = c2q(Yhcurrent_level(:,:,3 4),gain_mask(3 4,current_level); hh = c2q(Yhcurrent_level(:,:,2 5),gain_mask(2 5

17、,current_level); % Do even Qshift filters on columns. y1 = colifilt(Z,g0b,g0a) + colifilt(lh,g1b,g1a); y2 = colifilt(hl,g0b,g0a) + 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 required row_size col_size = s

18、ize(Z); S = 2*size(Yhcurrent_level-1); if row_size = S(1)%check to see if this result needs to be cropped for the rows Z = Z(2:row_size-1,:); end if col_size = S(2)%check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1); end if any(size(Z) = S(1:2), error(Sizes of subband

19、s are not valid for DTWAVEIFM2); end current_level = current_level - 1;endif current_level = 1; lh = c2q(Yhcurrent_level(:,:,1 6),gain_mask(1 6,current_level); hl = c2q(Yhcurrent_level(:,:,3 4),gain_mask(3 4,current_level); hh = c2q(Yhcurrent_level(:,:,2 5),gain_mask(2 5,current_level); % Do odd top

20、-level filters on columns. 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,gain)% function z = c2q(w,gain)% Scale by gain an

21、d convert from complex w(:,:,1:2) to real quad-numbers in z.% Arrange pixels from the real and imag parts of the 2 subbands% into 4 separate subimages .% A-B Re Im of w(:,:,1)% | |% | |% C-D Re Im of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2);if any(w(:) & any(gain) sc = sqrt(0.5) * gain; P = w(:,:,1)

22、*sc(1) + w(:,:,2)*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;

23、endreturn（3）function Yl,Yh,Yscale = dtwavexfm2(X,nlevels,biort,qshift);% Function to perform a n-level DTCWT-2D decompostion on a 2D matrix X% Yl,Yh,Yscale = dtwavexfm2(X,nlevels,biort,qshift);% X - 2D real matrix/Image% nlevels - No. of levels of wavelet decomposition% biort - antonini = Antonini 9

24、,7 tap filters.% legall = LeGall 5,3 tap filters.% near_sym_a = Near-Symmetric 5,7 tap filters.% near_sym_b = Near-Symmetric 13,19 tap filters.% qshift - qshift_06 = Quarter Sample Shift Orthogonal (Q-Shift) 10,10 tap filters, % (only 6,6 non-zero taps).% qshift_a = Q-shift 10,10 tap filters,% (with

25、 10,10 non-zero taps, unlike 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.% % Yl - The real lowpass image from the final level% Yh - A cell array containing the 6 complex highpass subimages for each level.% Yscale - T

26、his is an OPTIONAL output argument, that is a cell array containing % real lowpass coefficients for every scale.% % Example: Yl,Yh = dtwavexfm2(X,3,near_sym_b,qshift_b);% performs a 3-level transform on the real image X using the 13,19-tap filters % for level 1 and the Q-shift 14-tap filters for lev

27、els = 2.% Nick Kingsbury and Cian Shaffrey% Cambridge University, Sept if isstr(biort) & isstr(qshift)%Check if the inputs are strings biort_exist = exist(biort .mat); qshift_exist = exist(qshift .mat); if biort_exist = 2 & qshift_exist = 2; %Check to see if the inputs exist as .mat files load (bior

28、t); load (qshift); else error(Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTWAVEXFM2 for details.); endelse error(Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTWAVEXFM2.);end orginal_size = size(X);if ndims(X) = 3; error(sprintf(

29、The entered image is %dx%dx%d, please enter each image slice separately.,orginal_size(1),orginal_size(2),orginal_size(3);end% The next few lines of code check to see if the image is odd in size, if so an extra .% row/column will be added to the bottom/right of the imageinitial_row_extend = 0; %initi

30、aliseinitial_col_extend = 0;if any(rem(orginal_size(1),2), %if sx(1) is not divisable by 2 then we need to extend X by adding a row at the bottom X = X; X(end,:); %Any further extension will be done in due course. initial_row_extend = 1;endif any(rem(orginal_size(2),2), %if sx(2) is not divisable by

31、 2 then we need to extend X by adding a col to the left X = X X(:,end); %Any further extension will be done in due course. initial_col_extend = 1;endextended_size = size(X);if nlevels = 0, return; end%initialiseYh=cell(nlevels,1);if nargout = 3 Yscale=cell(nlevels,1); %this is only required if the u

32、ser specifies a third output component.endS = ;sx = size(X);if nlevels = 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,h0o).;% LoLo Yh1 = zeros(size(LoLo)/2 6); Yh1(:,:,1 6) = q2c(colfilter(Hi,h0o).);% H

33、orizontal pair Yh1(:,:,3 4) = q2c(colfilter(Lo,h1o).);% Vertical pair Yh1(:,:,2 5) = q2c(colfilter(Hi,h1o).); % Diagonal pair S = size(LoLo) ;S; if nargout = 3 Yscale1 = LoLo; endendif nlevels = 2; for level = 2:nlevels; row_size col_size = size(LoLo); if any(rem(row_size,4),% Extend by 2 rows if no

34、. of rows of LoLo are divisable by 4; LoLo = LoLo(1,:); LoLo; LoLo(end,:); end if any(rem(col_size,4),% Extend by 2 cols if no. of cols of LoLo are divisable by 4; LoLo = LoLo(:,1) LoLo LoLo(:,end); end % Do even Qshift filters on rows. Lo = coldfilt(LoLo,h0b,h0a).; Hi = coldfilt(LoLo,h1b,h1a).; % D

35、o even Qshift filters on columns. LoLo = coldfilt(Lo,h0b,h0a).;%LoLo Yhlevel = zeros(size(LoLo)/2 6); Yhlevel(:,:,1 6) = q2c(coldfilt(Hi,h0b,h0a).);% Horizontal Yhlevel(:,:,3 4) = q2c(coldfilt(Lo,h1b,h1a).);% Vertical Yhlevel(:,:,2 5) = q2c(coldfilt(Hi,h1b,h1a).);% Diagonal S = size(LoLo) ;S; if nar

36、gout = 3 Yscalelevel = LoLo; end endendYl = LoLo;if initial_row_extend = 1 & initial_col_extend = 1; warning(sprintf( rr The image entered is now a %dx%d NOT a %dx%d r The bottom row and rightmost column have been duplicated, prior to decomposition. rr ,. extended_size(1),extended_size(2),orginal_si

37、ze(1),orginal_size(2);endif initial_row_extend = 1 ; warning(sprintf( rr The image entered is now a %dx%d NOT a %dx%d r Row number %d has been duplicated, and added to the bottom of the image, prior to decomposition. rr,. extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size

38、(1);endif initial_col_extend = 1; warning(sprintf( rr The image entered is now a %dx%d NOT a %dx%d r Col number %d has been duplicated, and added to the right of the image, prior to decomposition. rr,. extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(2);endreturn%=%* IN

39、TERNAL FUNCTION *%=function z = q2c(y)% function z = q2c(y)% Convert from quads in y to complex numbers in z.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 alternate real and imag pixels.% a-b% | |% | |% c-d% Combin

40、e (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 subbands in z.z = cat(3,p-q,p+q);return（4）function Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Function t

41、o perform an n-level dual-tree complex wavelet (DTCWT)% 2-D reconstruction.% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% % Yl - The real lowpass image from the final level% Yh - A cell array containing the 6 complex highpass subimages for each level.% biort - antonini = Antonini 9,7 tap filters.%

42、 legall = LeGall 5,3 tap filters.% near_sym_a = Near-Symmetric 5,7 tap filters.% near_sym_b = Near-Symmetric 13,19 tap filters.% qshift - qshift_06 = Quarter Sample Shift Orthogonal (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

43、taps, unlike 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 - Gain to be applied to each subband. % gain_mask(d,l) is gain for subband with direction d at level l.% If gain_mask(d,l) = 0, no computation is p

44、erformed for band (d,l).% Default gain_mask = ones(6,length(Yh).% Z - Reconstructed real image matrix% % For example: Z = dtwaveifm2(Yl,Yh,near_sym_b,qshift_b);% performs a 3-level reconstruction from Yl,Yh using the 13,19-tap filters % for level 1 and the Q-shift 14-tap filters for levels = 2.% Nic

45、k Kingsbury and Cian Shaffrey% Cambridge University, May a = length(Yh); % No of levels.if nargin = 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yhcurrent_level(:,:,1 6),gain_mask(1 6,current_level); hl = c2q(Yhcurrent_level(:,:,3 4),gain_mask(3 4,current_level); hh = c2q(

46、Yhcurrent_level(:,:,2 5),gain_mask(2 5,current_level); % Do even Qshift filters on columns. y1 = colifilt(Z,g0b,g0a) + colifilt(lh,g1b,g1a); y2 = colifilt(hl,g0b,g0a) + colifilt(hh,g1b,g1a); % Do even Qshift filters on rows. Z = (colifilt(y1.,g0b,g0a) + colifilt(y2.,g1b,g1a).; % Check size of Z and

47、crop as required row_size col_size = size(Z); S = 2*size(Yhcurrent_level-1); if row_size = S(1)%check to see if this result needs to be cropped for the rows Z = Z(2:row_size-1,:); end if col_size = S(2)%check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1); end if any(si

48、ze(Z) = S(1:2), error(Sizes of subbands are not valid for DTWAVEIFM2); end current_level = current_level - 1;endif current_level = 1; lh = c2q(Yhcurrent_level(:,:,1 6),gain_mask(1 6,current_level); hl = c2q(Yhcurrent_level(:,:,3 4),gain_mask(3 4,current_level); hh = c2q(Yhcurrent_level(:,:,2 5),gain

49、_mask(2 5,current_level); % Do odd top-level filters on columns. 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,gain)% func

50、tion z = c2q(w,gain)% Scale by gain and convert from complex w(:,:,1:2) to real quad-numbers in z.% Arrange pixels from the real and imag parts of the 2 subbands% into 4 separate subimages .% A-B Re Im of w(:,:,1)% | |% | |% C-D Re Im of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2);if any(w(:) & any(gai

51、n) sc = sqrt(0.5) * gain; P = w(:,:,1)*sc(1) + w(:,:,2)*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(

52、t1+1,t2+1) = -real(Q); % d = (C-A)*sc;endreturn(5)function Yl,Yh,Yscale = wavexfm2(X,nlevels,biort);% Function to perform a n-level DWT-2D decompostion on a 2-D matrix X.% Yl,Yh,Yscale = dtwavexfm2(X,nlevels,biort);% X - real 1-D signal column vector (or matrix of vectors)% nlevels - No. of levels o

53、f wavelet decomposition% biort - antonini = Antonini 9,7 tap filters.% legall = LeGall 5,3 tap filters.% near_sym_a = Near-Symmetric 5,7 tap filters.% near_sym_b = Near-Symmetric 13,19 tap filters.% Yl - The lowpass subband from the final level.% Yh - A cell array containing the highpass subband for

54、 each level.% Yscale - This is an OPTIONAL output argument, that is a cell array containing % the lowpass coefficients at every scale.% % Example: Yl,Yh = wavexfm2(X,4,near_sym_b);% performs a 4-level 2-D DWT on the real image X using the 13,19-tap filters.% Nick Kingsbury, Cambridge University, May

55、 if isstr(biort)% Check if the biort input is a string biort_exist = exist(biort .mat); if biort_exist = 2, % Check to see if the filter exists as a .mat file load (biort); else error(Please enter the correct name of the Biorthogonal Filter, see help WAVEXFM2 for details.); endelse error(Please ente

56、r the name of the Biorthogonal Filter as shown in help WAVEXFM2.);endL = size(X);if any(rem(L,2), % ensure that X is an even length, thus enabling it to be extended if needs be. error(Size of X must be a multiple of 2);end%initialiseYh=cell(nlevels,1);if nargout = 3 Yscale=cell(nlevels,1); % This is

57、 only required if the user specifies a third output component.endLoLo = X;for level = 1:nlevels; if rem(size(LoLo,1),4),% Check to see if height of LoLo is divisable by 4, if not extend. LoLo = LoLo(1,:); LoLo; LoLo(end,:); end if rem(size(LoLo,2),4),% Check to see if height of LoLo is divisable by

58、4, if not extend. LoLo = LoLo(:,1) LoLo LoLo(:,end); end % Do filters on rows. Lo = coldwtfilt(LoLo,h0o,0).; Hi = coldwtfilt(LoLo,h1o,1).; % Do filters on columns. LoLo = coldwtfilt(Lo,h0o,0).;%LoLo Yhlevel = zeros(size(LoLo) 3); Yhlevel(:,:,1) = coldwtfilt(Hi,h0o,0).;% Horizontal Yhlevel(:,:,3) = coldwtfilt(Lo,h1o,1).;% Vertical Yhlevel(:,:,2) = coldwtfilt(Hi,h1o,1).;% Diagonal if nargout = 3 Yscalelevel = LoLo; e