📄 perform_wavelet_transform.m
字号:
% Description
% If wc is the result of a forward 2d wavelet transform, with
% wc = FWT2_PO(x,L,qmf), then x = IWT2_PO(wc,L,qmf) reconstructs x
% exactly if qmf is a nice qmf, e.g. one made by MakeONFilter.
%
% See Also
% FWT2_PO, MakeONFilter
%
[n,J] = quadlength(wc);
x = wc;
nc = 2^(L+1);
for jscal=L:J-1,
top = (nc/2+1):nc; bot = 1:(nc/2); all = 1:nc;
for iy=1:nc,
x(all,iy) = UpDyadLo(x(bot,iy)',qmf)' ...
+ UpDyadHi(x(top,iy)',qmf)';
end
for ix=1:nc,
x(ix,all) = UpDyadLo(x(ix,bot),qmf) ...
+ UpDyadHi(x(ix,top),qmf);
end
nc = 2*nc;
end
%
% Copyright (c) 1993. David L. Donoho
%
%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
%
function [n,J] = quadlength(x)
% quadlength -- Find length and dyadic length of square matrix
% Usage
% [n,J] = quadlength(x)
% Inputs
% x 2-d image; size(n,n), n = 2^J (hopefully)
% Outputs
% n length(x)
% J least power of two greater than n
%
% Side Effects
% A warning message is issue if n is not a power of 2,
% or if x is not a square matrix.
%
s = size(x);
n = s(1);
if s(2) ~= s(1),
disp('Warning in quadlength: nr != nc')
end
k = 1 ; J = 0; while k < n , k=2*k; J = 1+J ; end ;
if k ~= n ,
disp('Warning in quadlength: n != 2^J')
end
%
% Copyright (c) 1993. David L. Donoho
%
%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
%
function wp = FWT_TI(x,L,qmf)
% FWT_TI -- translation invariant forward wavelet transform
% Usage
% TIWT = FWT_TI(x,L,qmf)
% Inputs
% x array of dyadic length n=2^J
% L degree of coarsest scale
% qmf orthonormal quadrature mirror filter
% Outputs
% TIWT stationary wavelet transform table
% formally same data structure as packet table
%
% See Also
% IWT_TI
%
[n,J] = dyadlength(x);
D = J-L;
wp = zeros(n,D+1);
x = ShapeAsRow(x);
%
wp(:,1) = x';
for d=0:(D-1),
for b=0:(2^d-1),
s = wp(packet(d,b,n),1)';
hsr = DownDyadHi(s,qmf);
hsl = DownDyadHi(rshift(s),qmf);
lsr = DownDyadLo(s,qmf);
lsl = DownDyadLo(rshift(s),qmf);
wp(packet(d+1,2*b ,n),d+2) = hsr';
wp(packet(d+1,2*b+1,n),d+2) = hsl';
wp(packet(d+1,2*b ,n),1 ) = lsr';
wp(packet(d+1,2*b+1,n),1 ) = lsl';
end
end
%
% Copyright (c) 1994. David L. Donoho
%
%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
%
function tiwt = FWT2_TI(x,L,qmf)
% FWT_TI -- 2-D translation invariant forward wavelet transform
% Usage
% TIWT = FWT2_TI(x,L,qmf)
% Inputs
% x 2-d image (n by n real array, n dyadic)
% L degree of coarsest scale
% qmf orthonormal quadrature mirror filter
% Outputs
% TIWT translation-invariant wavelet transform table, (3*(J-L)+1)*n by n
%
% See Also
% IWT2_TI, IWT2_TIMedian
%
[n,J] = quadlength(x);
D = J-L;
tiwt = zeros((3*D+1)*n, n);
lastx = (3*D*n+1):(3*D*n+n); lasty = 1:n;
tiwt(lastx,lasty) = x;
%
for d=0:(D-1),
l = J-d-1; ns = 2^(J-d);
for b1=0:(2^d-1), for b2=0:(2^d-1),
s = tiwt(3*D*n+packet(d,b1,n), packet(d,b2,n));
wc00 = FWT2_PO(s,l,qmf);
wc01 = FWT2_PO(CircularShift(s,0,1),l,qmf);
wc10 = FWT2_PO(CircularShift(s,1,0),l,qmf);
wc11 = FWT2_PO(CircularShift(s,1,1),l,qmf);
index10 = packet(d+1,2*b1,n); index20 = packet(d+1,2*b2,n);
index11 = packet(d+1,2*b1+1,n); index21 = packet(d+1,2*b2+1,n);
% horizontal stuff
tiwt(3*d*n + index10 , index20) = wc00(1:(ns/2),(ns/2+1):ns);
tiwt(3*d*n + index11, index20) = wc01(1:(ns/2),(ns/2+1):ns);
tiwt(3*d*n + index10 , index21) = wc10(1:(ns/2),(ns/2+1):ns);
tiwt(3*d*n + index11 , index21) = wc11(1:(ns/2),(ns/2+1):ns);
% vertical stuff
tiwt((3*d+1)*n + index10 , index20) = wc00((ns/2+1):ns,1:(ns/2));
tiwt((3*d+1)*n + index11, index20) = wc01((ns/2+1):ns,1:(ns/2));
tiwt((3*d+1)*n + index10 , index21) = wc10((ns/2+1):ns,1:(ns/2));
tiwt((3*d+1)*n + index11 , index21) = wc11((ns/2+1):ns,1:(ns/2));
% diagonal stuff
tiwt((3*d+2)*n + index10 , index20) = wc00((ns/2+1):ns,(ns/2+1):ns);
tiwt((3*d+2)*n + index11, index20) = wc01((ns/2+1):ns,(ns/2+1):ns);
tiwt((3*d+2)*n + index10 , index21) = wc10((ns/2+1):ns,(ns/2+1):ns);
tiwt((3*d+2)*n + index11 , index21) = wc11((ns/2+1):ns,(ns/2+1):ns);
% low freq stuff
tiwt(3*D*n + index10 , index20) = wc00(1:(ns/2),1:(ns/2));
tiwt(3*D*n + index11, index20) = wc01(1:(ns/2),1:(ns/2));
tiwt(3*D*n + index10 , index21) = wc10(1:(ns/2),1:(ns/2));
tiwt(3*D*n + index11 , index21) = wc11(1:(ns/2),1:(ns/2));
end, end
end
%
% Copyright (c) 1995. David L. Donoho and Thomas P.Y. Yu
%
%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
%
function p=packet(d,b,n)
% packet -- Packet table indexing
% Usage
% p = packet(d,b,n)
% Inputs
% d depth of splitting in packet decomposition
% b block index among 2^d possibilities at depth d
% n length of signal
% Outputs
% p linear indices of all coeff's in that block
%
npack = 2^d;
p = ( (b * (n/npack) + 1) : ((b+1)*n/npack ) ) ;
%
% Copyright (c) 1993. David L. Donoho
%
%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
%
function x = IWT_TI(pkt,qmf)
% IWT_TI -- Invert translation invariant wavelet transform
% Usage
% x = IWT_TI(TIWT,qmf)
% Inputs
% TIWT translation-invariant wavelet transform table
% qmf quadrature mirror filter
% Outputs
% x 1-d signal reconstructed from translation-invariant
% transform TIWT
%
% See Also
% FWT_TI
%
[n,D1] = size(pkt);
D = D1-1;
J = log2(n);
L = J-D;
%
wp = pkt;
%
sig = wp(:,1)';
for d= D-1:-1:0,
for b=0:(2^d-1)
hsr = wp(packet(d+1,2*b ,n),d+2)';
hsl = wp(packet(d+1,2*b+1,n),d+2)';
lsr = sig(packet(d+1,2*b ,n) );
lsl = sig(packet(d+1,2*b+1,n) );
loterm = (UpDyadLo(lsr,qmf) + lshift(UpDyadLo(lsl,qmf)))/2;
hiterm = (UpDyadHi(hsr,qmf) + lshift(UpDyadHi(hsl,qmf)))/2;
sig(packet(d,b,n)) = loterm+hiterm;
end
end
x = sig;
%
% Copyright (c) 1994. David L. Donoho
%
%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
%
function x = IWT2_TI(tiwt,L,qmf)
% IWT2_TI -- Invert 2-d translation invariant wavelet transform
% Usage
% x = IWT2_TI(TIWT,qmf)
% Inputs
% TIWT translation-invariant wavelet transform table, (3*(J-L)+1)*n by n
% L degree of coarsest scale
% qmf quadrature mirror filter
% Outputs
% x 2-d image reconstructed from translation-invariant transform TIWT
%
% See Also
% FWT2_TI, IWT2_TIMedian
%
[D1,n] = size(tiwt);
J = log2(n);
D = J-L;
%
lastx = (3*D*n+1):(3*D*n+n); lasty = 1:n;
x = tiwt(lastx,lasty);
for d=(D-1):-1:0,
l = J-d-1; ns = 2^(J-d);
for b1=0:(2^d-1), for b2=0:(2^d-1),
index10 = packet(d+1,2*b1,n); index20 = packet(d+1,2*b2,n);
index11 = packet(d+1,2*b1+1,n); index21 = packet(d+1,2*b2+1,n);
wc00 = [x(index10,index20) tiwt(3*d*n+index10,index20) ; ...
tiwt((3*d+1)*n+index10,index20) tiwt((3*d+2)*n+index10,index20)];
wc01 = [x(index11,index20) tiwt(3*d*n+index11,index20) ; ...
tiwt((3*d+1)*n+index11,index20) tiwt((3*d+2)*n+index11,index20)];
wc10 = [x(index10,index21) tiwt(3*d*n+index10,index21) ; ...
tiwt((3*d+1)*n+index10,index21) tiwt((3*d+2)*n+index10,index21)];
wc11 = [x(index11,index21) tiwt(3*d*n+index11,index21) ; ...
tiwt((3*d+1)*n+index11,index21) tiwt((3*d+2)*n+index11,index21)];
x(packet(d,b1,n), packet(d,b2,n)) = ( IWT2_PO(wc00,l,qmf) + ....
CircularShift(IWT2_PO(wc01,l,qmf),0,-1) + ...
CircularShift(IWT2_PO(wc10,l,qmf),-1,0) + ...
CircularShift(IWT2_PO(wc11,l,qmf),-1,-1) ) / 4;
end, end
end
%
% Copyright (c) 1995. David L. Donoho and Thomas P.Y. Yu
%
%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
%
function result = CircularShift(matrix, colshift, rowshift)
% CIRCULARSHIFT: Circular shifting of a matrix/image, i.e., pixels that get
% shifted off one side of the image are put back on the other side.
%
% result = circularShift(matrix, colshift, rowshift)
%
% EPS, DJH '96
lastrow = size(matrix, 1);
lastcol = size(matrix, 2);
result = matrix;
% Shift the cols
if (colshift>0)
result = [result(:,[lastcol-colshift+1:lastcol]) ...
result(:,[1:lastcol-colshift])];
else
colshift = -colshift;
result = [result(:,[colshift+1:lastcol]) ...
result(:,[1:colshift])];
end
% Shift the rows
if (rowshift>0)
result = [result([lastrow-rowshift+1:lastrow],:) ; ...
result([1:lastrow-rowshift],:)];
else
rowshift = -rowshift;
result = [result([rowshift+1:lastrow],:) ; ...
result([1:rowshift],:)];
end
function x = reverse(x)
x = x(end:-1:1);
function StatWT = TI2Stat(TIWT)
% TI2Stat -- Convert Translation-Invariant Transform to Stationary Wavelet Transform
% Usage
% StatWT = TI2Stat(TIWT)
% Inputs
% TIWT translation invariant table from FWT_TI
% Outputs
% StatWT stationary wavelet transform table table as FWT_Stat
%
% See Also
% Stat2TI, FWT_TI, FWT_Stat
%
StatWT = TIWT;
[n,D1] = size(StatWT);
D = D1-1;
J = log2(n);
L = J-D;
%
index = 1;
for d=1:D,
nb = 2^d;
nk = n/nb;
index = [ (index+nb/2); index];
index = index(:)';
for b= 0:(nb-1),
StatWT(d*n + (index(b+1):nb:n)) = TIWT(d*n + packet(d,b,n));
end
end
for b= 0:(nb-1),
StatWT((index(b+1):nb:n)) = TIWT(packet(d,b,n));
end
%
% Copyright (c) 1994. Shaobing Chen
%
function TIWT = Stat2TI(StatWT)
% Stat2TI -- Convert Stationary Wavelet Transform to Translation-Invariant Transform
% Usage
% TIWT = Stat2TI(StatWT)
% Inputs
% StatWT stationary wavelet transform table as FWT_Stat
% Outputs
% TIWT translation-invariant transform table as FWT_TI
%
% See Also
% Stat2TI, FWT_TI, FWT_Stat
%
TIWT = StatWT;
[n,D1] = size(StatWT);
D = D1-1;
index = 1;
for d=1:D,
nb = 2^d;
nk = n/nb;
index = [ (index+nb/2); index];
index = index(:)';
for b= 0:(nb-1),
TIWT(d*n + packet(d,b,n)) = StatWT(d*n + (index(b+1):nb:n));
end
end
for b= 0:(nb-1),
TIWT(packet(d,b,n)) = StatWT((index(b+1):nb:n));
end
%
% Copyright (c) 1994. Shaobing Chen
%
%
% Part of Wavelab Version 850
% Built Tue Jan 3 13:20:40 EST 2006
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail wavelab@stat.stanford.edu
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -