function [x,level] = ldwt(x, operation, wavelet, extension, minsize)
if (nargin <= 1) operation = 'forward';end;
if (nargin <= 2) wavelet = 'legall53';end;
if (nargin <= 3) extension = 'zeropadded';end;
if (nargin <= 4) minsize = 8;end;
ls = lanalysis(wavelet,extension);
if strcmpi(operation,'forward');
[x,level] = wavedecnd(x,ls,minsize);
elseif strcmpi(operation,'inverse');
x = wavedecnd_inv(x,ls,minsize);
elseif strcmpi(operation,'adjoint');
x = wavedecnd_star(x,ls,minsize);
elseif strcmpi(operation,'inverseadjoint') || strcmpi(operation,'adjointinverse')
x = wavedecnd_star_inv(x,ls,minsize);
else
error('Unknown operation');
end
end
function [x, l] = wavedecnd(x,ls,m)
level = sublevel_sizes_all_sublevels(size(x),m);
l = level;
while size(level,1) > 1
N1vec = level(end,:);
I = selector(N1vec);
level = level(1:end-1,:);
x(I{:}) = lwtnd(x(I{:}),ls,m);
end
end
function x = wavedecnd_inv(x,ls,m)
level = sublevel_sizes_all_sublevels(size(x),m);
while size(level,1) > 1
level = level(2:end,:);
N1vec = level(1,:);
I = selector(N1vec);
x(I{:}) = lwtnd_inv(x(I{:}),ls,m);
end
end
function x = wavedecnd_star(x,ls,m)
level = sublevel_sizes_all_sublevels(size(x),m);
while size(level,1) > 1
level = level(2:end,:);
N1vec = level(1,:);
I = selector(N1vec);
x(I{:}) = lwtnd_star(x(I{:}),ls,m);
end
end
function x = wavedecnd_star_inv(x,ls,m)
level = sublevel_sizes_all_sublevels(size(x),m);
while size(level,1) > 1
N1vec = level(end,:);
I = selector(N1vec);
level = level(1:end-1,:);
x(I{:}) = lwtnd_star_inv(x(I{:}),ls,m);
end
end
function [x,N1vec] = lwtnd(x,ls,m)
N1vec = size(x);
for i =1:ndims(x)
s = size(x);
x = reshape(x,s(1),prod(s(2:end)));
if s(1) > m
[x,N1vec(i)] = lwt1d(x,ls);
end;
x = reshape(x,s);
x = shiftdim(x,1);
end;
end
function [x,N1vec] = lwtnd_inv(x,ls,m)
N1vec = size(x);
for i =1:ndims(x)
x = shiftdim(x,ndims(x)-1);
s = size(x);
x = reshape(x,s(1),prod(s(2:end)));
if s(1) > m
[x,N1vec(i)] = lwt1d_inv(x,ls);
end;
x = reshape(x,s);
end;
end
function [x,N1vec] = lwtnd_star(x,ls,m)
N1vec = size(x);
for i =1:ndims(x)
x = shiftdim(x,ndims(x)-1);
s = size(x);
x = reshape(x,s(1),prod(s(2:end)));
if s(1) > m
[x,N1vec(i)] = lwt1d_star(x,ls);
end;
x = reshape(x,s);
end;
end
function [x,N1vec] = lwtnd_star_inv(x,ls,m)
N1vec = size(x);
for i =1:ndims(x)
x = shiftdim(x,1);
s = size(x);
x = reshape(x,s(1),prod(s(2:end)));
if s(1) > m
[x,N1vec(i)] = lwt1d_star_inv(x,ls);
end;
x = reshape(x,s);
end;
end
function ls = lanalysis(scheme,extension)
if strcmpi(scheme,'haar');
K=sqrt(2.0);
ls = [0.0 -1.0 1;...
0.5 0.0 0;...
K K 0];
elseif strcmpi(scheme,'legall53')
a = -0.5;
b = 0.25;
K = sqrt(2.0);
ls = [a a 1;...
b b 0;...
K K 0];
elseif strcmpi(scheme,'cdf97')
a = -1.586134342E+0;
b = -5.298011854E-2;
c = 0.8829110762E+0;
d = 0.4435068522E+0;
K = 1.149604398E0;
ls = [a a 1;...
b b 0;...
c c 1;...
d d 0;...
K K 0];
else
error('Unknown lifting scheme.');
end
if strcmpi(extension,'symmetric')
ls(end,end) = 1;
elseif strcmpi(extension, 'zeropadded')
ls(end,end) = 0;
else
error('Unknown extension');
end
end
function [x,N1] = lwt1d(x,ls)
N = size(x,1);
N1 = sublevel_sizes(N);
extension = ls(end,end);
for i = 1:size(ls,1)-1
type = ls(i,3);
filt = ls(i,1:2);
if (extension == 0)
eb = 0; ee = 0;
elseif (extension == 1)
if (type == 0)
[eb,ee] = extension_high(filt(1),filt(2),N);
else
[eb, ee] = extension_low(filt(1),filt(2),N);
end;
else
error('Unknown extension');
end;
if (type == 0)
x = lift_row_high(x, eb, filt(1), filt(2), ee);
elseif (type == 1)
x = lift_row_low(x, eb, filt(1), filt(2), ee);
else
error('Unknown step');
end;
end
K = ls(end,1); L = ls(end,2); x = scale_row(x, K,L);
x = pack_row(x);
end
function [x,N1] = lwt1d_star(x,ls)
N = size(x,1);
N1 = sublevel_sizes(N);
x = pack_row_star(x);
K = ls(end,1); L = ls(end,2); x = scale_row_star(x, K,L);
extension = ls(end,end);
for i = size(ls,1)-1:-1:1
type = ls(i,3);
filt = ls(i,1:2);
if (extension == 0)
eb = 0; ee = 0;
elseif (extension == 1)
if (type == 0)
[eb,ee] = extension_high(filt(1),filt(2),N);
else
[eb, ee] = extension_low(filt(1),filt(2),N);
end;
else
error('Unknown extension');
end;
if (type == 0)
x = lift_row_high_star(x, eb, filt(1), filt(2), ee);
elseif (type == 1)
x = lift_row_low_star(x, eb, filt(1), filt(2), ee);
else
error('Unknown step');
end;
end
end
function [x,N1] = lwt1d_inv(x,ls)
N = size(x,1);
N1 = sublevel_sizes(N);
x = pack_row_inv(x);
K = ls(end,1); L = ls(end,2); x = scale_row_inv(x, K,L);
extension = ls(end,end);
for i = size(ls,1)-1:-1:1
type = ls(i,3);
filt = ls(i,1:2);
if (extension == 0)
eb = 0; ee = 0;
elseif (extension == 1)
if (type == 0)
[eb,ee] = extension_high(filt(1),filt(2),N);
else
[eb, ee] = extension_low(filt(1),filt(2),N);
end;
else
error('Unknown extension');
end;
if (type == 0)
x = lift_row_high_inv(x, eb, filt(1), filt(2), ee);
elseif (type == 1)
x = lift_row_low_inv(x, eb, filt(1), filt(2), ee);
else
error('Unknown step');
end;
end
end
function [x,N1] = lwt1d_star_inv(x,ls)
N = size(x,1);
N1 = sublevel_sizes(N);
extension = ls(end,end);
for i = 1:size(ls,1)-1
type = ls(i,3);
filt = ls(i,1:2);
if (extension == 0)
eb = 0; ee = 0;
elseif (extension == 1)
if (type == 0)
[eb,ee] = extension_high(filt(1),filt(2),N);
else
[eb, ee] = extension_low(filt(1),filt(2),N);
end;
else
error('Unknown extension');
end;
if (type == 0)
x = lift_row_high_star_inv(x, eb, filt(1), filt(2), ee);
elseif (type == 1)
x = lift_row_low_star_inv(x, eb, filt(1), filt(2), ee);
else
error('Unknown step');
end;
end
K = ls(end,1); L = ls(end,2); x = scale_row_star_inv(x, K,L);
x = pack_row_star_inv(x);
end
function [eb,ee] = extension_low(a,b,N)
eb = 0; if (mod(N,2) == 0) ee = a; else ee = 0; end;
end
function [eb,ee] = extension_high(c,d,N)
eb = d; if (mod(N,2) == 0) ee = 0; else ee = c; end;
end
function x = pack_row_star_inv( x)
x = pack_row(x);
end
function x = pack_row_star( x)
x = pack_row_inv(x);
end
function x = pack_row_inv( x)
[I1, I2, N1, N2, N] = split_index(x);
t = x(N1+1:end,:);
x(I1,:) = x(1:N1,:);
x(I2,:) = t;
end
function x = pack_row( x)
[I1, I2, N1, N2, N] = split_index(x);
x = [x(I1,:);x(I2,:)];
end
function x = scale_row_star_inv( x, K,L)
x = scale_row(x, 1/K,1/L);
end
function x = scale_row_star( x, K,L)
x = scale_row(x, K,L);
end
function x = scale_row_inv( x, K,L)
x = scale_row(x, 1/K,1/L);
end
function x = scale_row( x, K,L)
[I1, I2, N1, N2, N] = split_index(x);
x(I1,:) = x(I1,:)*K;
x(I2,:) = x(I2,:)/L;
end
function x = lift_row_low_star_inv( x, eb,a,b,ee)
x = lift_row_high(x, -eb,-b,-a,-ee);
end
function x = lift_row_low_star( x, eb,a,b,ee)
x = lift_row_high(x, eb,b,a,ee);
end
function x = lift_row_low_inv( x, eb,a,b,ee)
x = lift_row_low(x, -eb,-a,-b,-ee);
end
function x = lift_row_low( x, eb,a,b,ee)
[I1, I2, N1, N2, N] = split_index(x);
if N > 2
x(I2(1),:) = x(I2(1),:) + (b+eb)*x(I1(1),:) + a*x(I1(2),:);
x(I2(2:N2-1),:) = x(I2(2:N2-1),:) + b*x(I1(2:N2-1),:) + a*x(I1(3:N2),:);
if (mod(N,2) == 0)
x(I2(N2),:) = x(I2(N2),:) + (b+ee)*x(I1(N1),:);
else
x(I2(N2),:) = x(I2(N2),:) + b*x(I1(N1-1),:) + (a+ee)*x(I1(N1),:);
end;
else
x(I2(1),:) = x(I2(1),:) + (b+eb+ee)*x(I1(1),:);
end
end
function x = lift_row_high_star_inv( x, eb,c,d,ee)
x = lift_row_low(x, -eb,-d,-c,-ee);
end
function x = lift_row_high_star( x, eb,c,d,ee)
x = lift_row_low(x, eb,d,c,ee);
end
function x = lift_row_high_inv( x, eb,c,d,ee)
x = lift_row_high(x, -eb,-c,-d,-ee);
end
function x = lift_row_high( x, eb,c,d,ee)
[I1, I2, N1, N2, N] = split_index(x);
if N > 2
x(I1(1),:) = x(I1(1),:) + (c+eb)*x(I2(1),:);
x(I1(2:N1-1),:) = x(I1(2:N1-1),:) + d*x(I2(1:N1-2),:) + c*x(I2(2:N1-1),:);
if (mod(N,2) == 0)
x(I1(N1),:) = x(I1(N1),:) + d*x(I2(N1-1),:) +(c+ee)*x(I2(N1),:);
else
x(I1(N1),:) = x(I1(N1),:) + (d+ee)*x(I2(N1-1),:);
end;
else
x(I1(1),:) = x(I1(1),:) + (c+eb+ee)*x(I2(1),:);
end
end
function [I1, I2, N1, N2, N] = split_index(x)
N = size(x,1);
[N1,N2] = sublevel_sizes(N);
I1 = 1:2:2*N1-1;
I2 = 2:2:2*N2;
end
function [N1,N2] = sublevel_sizes(N)
N1 = ceil(N/2);
N2 = N - N1;
end
function level = sublevel_sizes_all_sublevels(s, m)
level = s;
while max(s) > m
for i = 1:length(s)
if s(i) > m
s(i) = sublevel_sizes(s(i));
end;
end;
level = [s;level];
end
end
function I = selector(N1vec)
for i = 1:length(N1vec)
I{i} = 1:N1vec(i);
end
end