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

一、题目：双树复小波变换

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

?c(t)??h(t)?j?g(t)f,?c?f,?h?jf,?g

22幅值: |dc(j,n)|?[dh(j,n)]?[dg(j,n)]相位： ?d(j,n)?arctan?dg(j,n)???cd(j,n)h??

2． 代码实现。

我采用Matlab函数编程实现。具体程序见shift_test_2D.m，drawcirc.m，setfig.m，dtwavexfm2.m，dtwaveifm2.m,waveifm2.m,wavexfm2.m

SkelMap.m。

n ，设 ? h , ? h和 ? h , ? h 分别是双正交对偶尺度函数与对偶小波， h 0 h 0 n ， h 1 n

和 h 1 n 是相应的低通滤波器和高通滤波器，即它们满足 g0(n)?g(2t?n)实部： ? ( t 2 h ( ? (2 t ? n ) 虚部：? g(t)?2h)?0n)hnn

?g(t)?2g1(n)?g(2t?n)?h(t)?2h1(n)?h(2t?n)nn

?h(t)?2h0(n)?h(2t?n)?g(t)?2g0(n)?g(2t?n)nn

?h(t)?2h1(n)?h(2t?n)?g(t)?2g1(n)?g(2t?n)nn

双树复小波变换可以通过离散小波变换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 all close 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]);

????????????????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);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 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(x,2)/2; end

figure;

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 on

plot(128*[[1;1]*[1:4] [0;6]]+1,128*[[0;4]*[1 1 1 1] [2;2]]+1,'-k'); hold off

title('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'); end

text(5*128+1,size(yy,1)*1.02,'level 4 scaling fn.','FontSize',14,'horiz','c','vert','t'); drawnow

% print -deps circrecq.eps

disp('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); end

figure;

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');

drawnow

return

（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 2002

a = 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 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 Filters, see help DTWAVEIFM2 for details.'); end else

error('Please enter the names of the Biorthogonal or Q-Shift Filters as shown in help DTWAVEIFM2.'); end

current_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);