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

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

2、t % 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) - * 200);setfig(1); colormap(gray(256)image(min(max(x+128,1),256);set(gca,position, .25 .5);axis(off);a

3、xis(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); y

4、t2 = 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 with DWT.% Do 4 levels of Real DWT using antonini (9,7)-tap fil

5、ters.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(x,2)/2;endfigure;setfig(gcf); colormap(gray(256)image(min(ma

6、x(yy+128,1),256);set(gca,position, .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*size(yy,2),*size(yy,1),DT CWT,horiz,r);text*size(yy,2),*size(yy,1),wavelets:,horiz,r,vert,t);text*size(yy,2)

7、,*size(yy,1),DWT,horiz,r);for k=1:4, text(k*128-63,size(yy,1)*,sprintf(level %d,k),FontSize,14,horiz,c,vert,t); endtext(5*128+1,size(yy,1)*,level 4 scaling fn.,FontSize,14,horiz,c,vert,t);drawnow% print -deps disp(Press a key to see perfect reconstruction property .)pause% Accumulate the images from

8、 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, .8 .8);axis(off);axis(image);title(Accumulated reconstructions fr

9、om each level of DT CWT ,FontSize,14);text(size(yy,2)*,size(yy,1)*,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 p

10、els 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 = + * sin(min(max(exp * (x + y ) - exp)*(r*3/w), -pi/2), pi/2);return（3）

11、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 subima

12、ges 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 ta

13、ps).% 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 d

14、irection 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

15、 % 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 for level -1 we never do the following lh = c2q(Yhcurrent_level(:,:,1 6),gain_mask(1 6,current_level); hl

16、= 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) + colifilt(lh,g1b,g1a); y2 = colifilt(hl,g0b,g0a) + colifilt(hh,g1b,g1a); % Do even Qshift filters on rows. Z = (

17、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 needs to be cropped for the rows Z = Z(2:row_size-1,:); end if col_size = S(2) %check to see if this resu

18、lt 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_level = 1; lh = c2q(Yhcurrent_level(:,:,1 6),gain_mask(1 6,current_level); hl = c2q(Yhcurrent_level(:,:,3

19、 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(hl,g0o) + colfilter(hh,g1o); % Do odd top-level filters on rows. Z = (colfilter(y1.,g0o) + colfilter(y2.,

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

21、D Re Im of w(:,:,2)sw = size(w);x = zeros(2*sw(1:2);if any(w(:) & any(gain) sc = sqrt * 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)

22、 = 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);%

23、 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 Orthogo

24、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, 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

25、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 tran

26、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 Kingsbury and Cian Shaffrey% Cambridge University, Sept 2001if isstr(biort) & isstr(qshift) %Check if the inputs are strings biort_exist = exist(biort .mat); qshift_exist = exist(

27、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 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(The entered image is %dx%dx%d, pl