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

一、题目：双树复小波变换二、目的：双树复小波和实小波变换的比较三、算法及其实现：提取阶梯型边界点1．算法。

幅值:相位：2．代码实现。

我采用Matlab函数编程实现。

具体程序见shift_test_2D.m，drawcirc.m，setfig.m，dtwavexfm2.m，dtwaveifm2.m,waveifm2.m,wavexfm2.mSkelMap.m。

设和分别是双正交对偶尺度函数与对偶小波，，，和实部：虚部：双树复小波变换可以通过离散小波变换DWT实现：一个DWT产生实部，另一个产生虚部。

四、实现工具：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 invariance properties.%% Nick Kingsbury, Cambridge University, May 2002.clear allclose all% Draw a circular disc.x = round((drawcirc(64,1,0,0,256) - 0.5) * 200);setfig(1);colormap(gray(256))image(min(max(x+128,1),256));set(gca,'position',[0.1 0.25 .25 .5]);()ctψ=,,j,c h gf fψψψ=+|(,)|cd j n=(,)arctan(,)gchd j nd j n∠= ⎪⎝⎭,h hφφ%,h hψψ%()h n()h n%()1h n()1h n%11()()(2)()()(2)()()(2)()()(2)h hnh hnh hnh hnt h n t nt h n t nt h n t nt h n t nφφψφφφψφ=-=-=-=-%%%%%%11()()(2)()()(2)()()(2)()()(2)g gng gng gng gnt g n t nt g n t nt g n t nt g n t nφφψφφφψφ=-=-=-=-%%%%%%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 turn.% 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);drawnowyy(yt1,yt2) = z;yt2 = yt2 + size(x,2)/2;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;draw(z);drawnowyy(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);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 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(:,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 ','FontSize',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 (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 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.% '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 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 performed 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.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, May 2002a = length(Yh); % No of levels.if nargin < 5, gain_mask = ones(6,a); end % Default gain_mask.if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_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 (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEIFM2.');endcurrent_level = a;Z = Yl;while current_level >= 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_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 crop as required[row_size col_size] = size(Z);S = 2*size(Yh{current_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,:);endif col_size ~= S(2) %check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1);endif any(size(Z) ~= S(1:2)),error('Sizes of subbands are not valid for DTW A VEIFM2');endcurrent_level = current_level - 1;endif current_level == 1;lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_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)% function 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(gain)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(t1+1,t2+1) = -real(Q); % d = (C-A)*sc;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,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 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 -> This 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 levels >= 2.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, Sept 2001if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_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 filesload (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEXFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEXFM2.');endorginal_size = size(X);if ndims(X) >= 3;error(sprintf('The entered image is %dx%dx%d, please enter each image sliceseparately.',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; %initialiseinitial_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 bottomX = [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 2 then we need to extend X by adding a col to the leftX = [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 == 3Yscale=cell(nlevels,1); %this is only required if the user 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).'; % LoLoYh{1} = zeros([size(LoLo)/2 6]);Yh{1}(:,:,[1 6]) = q2c(colfilter(Hi,h0o).'); % Horizontal pairYh{1}(:,:,[3 4]) = q2c(colfilter(Lo,h1o).'); % Vertical pairYh{1}(:,:,[2 5]) = q2c(colfilter(Hi,h1o).'); % Diagonal pairS = [ size(LoLo) ;S];if nargout == 3Yscale{1} = 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. of rows of LoLo are divisable by 4;LoLo = [LoLo(1,:); LoLo; LoLo(end,:)];endif 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).';% Do even Qshift filters on columns.LoLo = coldfilt(Lo,h0b,h0a).'; %LoLoYh{level} = zeros([size(LoLo)/2 6]);Yh{level}(:,:,[1 6]) = q2c(coldfilt(Hi,h0b,h0a).'); % HorizontalYh{level}(:,:,[3 4]) = q2c(coldfilt(Lo,h1b,h1a).'); % VerticalYh{level}(:,:,[2 5]) = q2c(coldfilt(Hi,h1b,h1a).'); % DiagonalS = [ size(LoLo) ;S];if nargout == 3Yscale{level} = LoLo;endendendYl = LoLo;if initial_row_extend == 1 & initial_col_extend == 1;warning(sprintf(' \r\r 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. \r\r ',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2)));endif initial_row_extend == 1 ;warning(sprintf(' \r\r 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. \r\r',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(1)));endif initial_col_extend == 1;warning(sprintf(' \r\r 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. \r\r',...extended_size(1),extended_size(2),orginal_size(1),orginal_size(2),orginal_size(2)));endreturn%==================================================================== ======================% ********** INTERNAL 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% Combine (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 to 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.% '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 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 performed 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.%% Nick Kingsbury and Cian Shaffrey% Cambridge University, May 2002a = length(Yh); % No of levels.if nargin < 5, gain_mask = ones(6,a); end % Default gain_mask.if isstr(biort) & isstr(qshift) %Check if the inputs are stringsbiort_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 (biort);load (qshift);elseerror('Please enter the correct names of the Biorthogonal or Q-Shift Filters, see help DTW A VEIFM2 for details.');endelseerror('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTW A VEIFM2.');endcurrent_level = a;Z = Yl;while current_level >= 2; ; %this ensures that for level -1 we never do the following lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_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 crop as required[row_size col_size] = size(Z);S = 2*size(Yh{current_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,:);endif col_size ~= S(2) %check to see if this result needs to be cropped for the cols Z = Z(:,2:col_size-1);endif any(size(Z) ~= S(1:2)),error('Sizes of subbands are not valid for DTW A VEIFM2');endcurrent_level = current_level - 1;endif current_level == 1;lh = c2q(Yh{current_level}(:,:,[1 6]),gain_mask([1 6],current_level));hl = c2q(Yh{current_level}(:,:,[3 4]),gain_mask([3 4],current_level));hh = c2q(Yh{current_level}(:,:,[2 5]),gain_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)% function 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(gain)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(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 of 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 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 2002if isstr(biort) % Check if the biort input is a stringbiort_exist = exist([biort '.mat']);if biort_exist == 2, % Check to see if the filter exists as a .mat fileload (biort);elseerror('Please enter the correct name of the Biorthogonal Filter, see help W A VEXFM2 for details.');endelseerror('Please enter the name of the Biorthogonal Filter as shown in help W A VEXFM2.');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 == 3Yscale=cell(nlevels,1); % This is 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,:)];endif rem(size(LoLo,2),4), % Check to see if height of LoLo is divisable by 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).'; %LoLoYh{level} = zeros([size(LoLo) 3]);Yh{level}(:,:,1) = coldwtfilt(Hi,h0o,0).'; % HorizontalYh{level}(:,:,3) = coldwtfilt(Lo,h1o,1).'; % VerticalYh{level}(:,:,2) = coldwtfilt(Hi,h1o,1).'; % Diagonalif nargout == 3Yscale{level} = LoLo;endendYl = LoLo;return(6)function Z = waveifm2(Yl,Yh,biort,gain_mask);% Function to perform an n-level DWT 2-D reconstruction.%% Z = waveifm2(Yl,Yh,biort,gain_mask);%% Yl -> The real lowpass subband from the final level% Yh -> A cell array containing the complex highpass subband for each level.%% 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.%% gain_mask -> Gain to be applied to each subband.% gain_mask(l) is gain for wavelet subband at level l.% If gain_mask(l) == 0, no computation is performed for band (l).% Default gain_mask = ones(1,length(Yh)).%% Z -> Reconstructed real signal vector (or matrix).%%% For example: Z = waveifm2(Yl,Yh,'near_sym_b');% performs a reconstruction from Yl,Yh using the 13,19-tap filters.%% Nick Kingsbury, Cambridge University, May 2002nlevels = length(Yh); % No of levels.if nargin < 4, gain_mask = ones(3,nlevels); end % Default gain_mask.if isstr(biort) % Check if the biort input is a stringbiort_exist = exist([biort '.mat']);if biort_exist == 2, % Check to see if the filter exists as a .mat fileload (biort);elseerror('Please enter the correct name of the Biorthogonal Filter, see help W A VEIFM2 for details.');endelseerror('Please enter the name of the Biorthogonal Filter as shown in help W A VEIFM2.');endLoLo = Yl;for level = nlevels:-1:1, % Reconstruct levels in reverse order.if size(LoLo,1) ~= size(Yh{level},1) % If LoLo is not the same height as the next Yh => t1 was extended.LoLo = LoLo(2:size(LoLo,1)-1,:); % Therefore we have to clip LoLo so it is the same height as the next Yh.endif size(LoLo,2) ~= size(Yh{level},2) % If LoLo is not the same width as the next Yh => t1was extended.LoLo = LoLo(:,2:size(LoLo,2)-1); % Therefore we have to clip LoLo so it is the same width as the next Yh.endif any([size(LoLo) 3] ~= size(Yh{level})),error('Yh sizes are not valid for WA VEIFM2');endlh = Yh{level}(:,:,1) * gain_mask(1,level);hl = Yh{level}(:,:,3) * gain_mask(3,level);hh = Yh{level}(:,:,2) * gain_mask(2,level);% Do even Qshift filters on columns.y1 = coliwtfilt(LoLo,2*g0o,0) + coliwtfilt(lh,2*g1o,1);y2 = coliwtfilt(hl,2*g0o,0) + coliwtfilt(hh,2*g1o,1);% Do even Qshift filters on rows.LoLo = (coliwtfilt(y1.',2*g0o,0) + coliwtfilt(y2.',2*g1o,1)).';endZ = LoLo;return(7)function setfig(f)% function setfig(f)% Set figure dimensions for correct display of text to agree with% print -deps.figure(f)set(gcf,'position',[220 40 790 526],'DefaultTextFontSize',14);return六、运行结果：七、结果分析：。