小波实验报告材料双树复小波变换.doc

《小波实验报告材料双树复小波变换.doc》由会员分享，可在线阅读，更多相关《小波实验报告材料双树复小波变换.doc（20页珍藏版）》请在三一办公上搜索。

1、word一、题目：双树复小波变换二、目的：双树复小波和实小波变换的比拟三、算法与其实现：提取阶梯型边界点1 算法。幅值:相位：2 代码实现。我采用Matlab函数编程实现。具体程序见shift_test_2D.m，drawcirc.m，setfig.m，dtwavexfm2.m。设 和 分别是双正交对偶尺度函数与对偶小波， ， ， 和 是相应的低通滤波器和高通滤波器，即它们满足实部： 虚部：双树复小波变换可以通过离散小波变换DWT实现：一个DWT产生实部，另一个产生虚部。4、 实现工具：Matlab5、 程序代码： 1shift_test_2D.m：% M-file to perform a

2、4-level wavelet transform on a circle using Q-shift % dual wavelet tree and DWT, and to pare 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

3、(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 turn.% Starts with the finest level.titl = 1st;2nd

4、;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;end% disp(Press a key .)% pause% Now do same wit

5、h 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);drawnow yy(yt1,yt2) = z; yt2 = yt2 + size(

6、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(ponents of reconstructed disc images,FontSize,14);text(-0.01*size(yy,2),0.25*size(yy,1),DT CWT,hori

7、z,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);drawnowdisp

8、(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),2

9、56);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);drawnowreturn2function p = drawcirc(r,w,dx,dy,N)% func

10、tion 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

11、)/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);return3function Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);% Function to perform an n-level dual-tree plex wavelet (DTCWT)% 2-D reconstruction.% Z = dtwaveifm2(Yl,Yh,biort,qshift,gain_mask);

12、% % Yl - The real lowpass image from the final level% Yh - A cell array containing the 6 plex 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 filter

13、s.% 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

14、,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 putation is performed for band (d,l).% Default gain_mask = ones(6,length(Yh).% Z - Reconstructed real image matrix% % For example: Z = dtwave

15、ifm2(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 = 2; ; %this ensures that

16、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,current_level); % Do even Qshift filters on columns. y1 = colifilt(Z,g0b,g0a) + colif

17、ilt(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(Yhcurrent_level-1); if row_size = S(1)%check to see if this result

18、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 subbands are not valid for DTWAVEIFM2); end current_level = current_level - 1;endif current_l

19、evel = 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-level filters on columns. y1 = colfilter(Z,g0o) + colfilter(lh,g1o); y2 = colfilter(h

20、l,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 plex w(:,:,1:2) to real quad-numbers in z.% Arrange pixels from the real

21、 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

22、: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;endreturn3function Yl,Yh,Yscale = dtwavexfm2(X,nlevels,biort,qshift);% Function to perform

23、 a n-level DTCWT-2D depostion on a 2D matrix X% Yl,Yh,Yscale = dtwavexfm2(X,nlevels,biort,qshift);% X - 2D real matrix/Image% nlevels - No. of levels of wavelet deposition% biort - antonini = Antonini 9,7 tap filters.% legall = LeGall 5,3 tap filters.% near_sym_a = Near-Symmetric 5,7 tap filters.% n

24、ear_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-Shif

25、t 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 plex highpass subimages for each level.% Yscale - This is an OPTIONAL output argument, that is a cell array containing % real lowpass coefficients for e

26、very 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(qs

27、hift)%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 (biort); load (qshift); else error(Please enter the correct names of the Biorthogonal or Q-Shift Filte

28、rs, 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(The entered image is %dx%dx%d, please enter each image slice separately.,orginal_size(1),orginal_

29、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 w

30、e 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 2 then we need to extend X by adding a col to the left X = X X(:,end); %Any further extension wi

31、ll 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 user specifies a third output ponent.endS = ;sx = size(X);if nlevels = 1, % Do odd top-level filte

32、rs 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).);% Horizontal pair Yh1(:,:,3 4) = q2c(colfilter(Lo,h1o).);% Vertical pair Yh1(:,:,2 5) = q2c(colfilter(H

33、i,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. of rows of LoLo are divisable by 4; LoLo = LoLo(1,:); LoLo; LoLo(end,:); end if any(rem(col_size,4

34、),% 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).;%LoLo Yhlevel = zeros(size(LoLo)/2 6)

35、; 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 nargout = 3 Yscalelevel = LoLo; end endendYl = LoLo;if initial_row_extend = 1 & initial_col_extend = 1;

36、 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 deposition. rr ,. extended_size(1),extended_size(2),orginal_size(1),orginal_size(2);endif initial_row_extend = 1 ; warning(sprintf( rr The image entered is now a %dx

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

38、umber %d has been duplicated, and added to the right of the image, prior to deposition. rr,. 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 plex 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% bine