# The range of the coefficients of 2 level Wavelet LeGall 5/3 2D transform

As in the title, I am confused about the range of the coefficients of Wavelet LeGall 5/3 (has to be exact this filter) 2D transform (only for a 8*8 block) if the value of the input matrix are within the range from 0-255. The link is here:wavelet legall 5/3

Here is what I did for now:

Minus 128 for all value (easier to calculate the low frequency values, see later); Do the transform in horizontal direction. This will generate all coefficients in all lines: the first 4 are low frequency and last 4 are high frequency. It is easy to find the range for high frequency is -255 to +255 (double the range of input). And the range for low frequency is actually -192 to +192 (1.5* the range of input). Do the transform in vertical direction. This will do the same in vertical directions. And there are four blocks generated: LL (lowlow), LH (low high), HL, HH. It is easy to calculate the range for the HH is largest: -511 to +511 and the range for LL is 1.5*1.5 = 2.25 times of the input range (-128 to +127). Then, here comes the question. What if I do this wavelet again for the LL block? Theoretically the range of the HH (2nd level coefficients) of the LL block should be 4 times of the LL range which becomes 10 times of the input range (-128 to +127) which is -1280 to +1270.

However, I tried many times random calculation, the max value never exceed -511 to +511 (see code at the end). I guess it is because the theoretical values cannot be reached because all calculations are based on the previous one. But this seemly easy question is difficult for me to prove theoretically. So can someone help me get out please?

Code I used (put two files in one directory and run the test file for any times you want but just the max value will not exceed 512...):

1. Function of waveletlegall53:

function X = waveletlegall53(X, Level)
%WAVELETLEGALL53  Le Gall 5/3 (Spline 2.2) wavelet transform.
%   Y = WAVELETLEGALL53(X, L) decomposes X with L stages of the
%   Le Gall 5/3 wavelet.  For the inverse transform,
%   WAVELETLEGALL53(X, -L) inverts L stages.  Filter boundary
%   handling is half-sample symmetric.
%
%   X may be of any size; it need not have size divisible by 2^L.
%   For example, if X has length 9, one stage of decomposition
%   produces a lowpass subband of length 5 and a highpass subband
%   of length 4.  Transforms of any length have perfect
%   reconstruction (exact inversion).
%
%   If X is a matrix, WAVELETLEGALL53 performs a (tensor) 2D
%   wavelet transform.  If X has three dimensions, the 2D
%   transform is applied along the first two dimensions.
%
%   Example:
%   Y = waveletlegall53(X, 5);    % Transform image X using 5 stages
%   R = waveletlegall53(Y, -5);   % Reconstruct from Y

X = double(X);
if nargin < 2, error('Not enough input arguments.'); end
if ndims(X) > 3, error('Input must be a 2D or 3D array.'); end
if any(size(Level) ~= 1), error('Invalid transform level.'); end

N1 = size(X,1);
N2 = size(X,2);

% Lifting scheme filter coefficients for Le Gall 5/3
LiftFilter = [-1/2,1/4];
ScaleFactor =1; sqrt(2);

LiftFilter = LiftFilter([1,1],:);

if Level >= 0   % Forward transform
for k = 1:Level
M1 = ceil(N1/2);
M2 = ceil(N2/2);

%%% Transform along columns %%%
if N1 > 1
RightShift = [2:M1,M1];
X0 = X(1:2:N1,1:N2,:);

% Apply lifting stages
if rem(N1,2)
X1 = [X(2:2:N1,1:N2,:);X0(M1,:,:)]...
+floor(filter(LiftFilter(:,1),1,X0(RightShift,:,:),...
X0(1,:,:)*LiftFilter(1,1),1));
else
X1 = X(2:2:N1,1:N2,:) ...
+floor(filter(LiftFilter(:,1),1,X0(RightShift,:,:),...
X0(1,:,:)*LiftFilter(1,1),1));
end

X0 = X0 + floor(filter(LiftFilter(:,2),1,...
X1,X1(1,:,:)*LiftFilter(1,2),1)+0.5);

if rem(N1,2)
X1(M1,:,:) = [];
end

X(1:N1,1:N2,:) = [X0*ScaleFactor;X1/ScaleFactor];

end

%%% Transform along rows %%%
if N2 > 1
RightShift = [2:M2,M2];
X0 = permute(X(1:N1,1:2:N2,:),[2,1,3]);

% Apply lifting stages
if rem(N2,2)
X1 = permute([X(1:N1,2:2:N2,:),X(1:N1,N2,:)],[2,1,3])...
+ floor(filter(LiftFilter(:,1),1,X0(RightShift,:,:),...
X0(1,:,:)*LiftFilter(1,1),1));
else
X1 = permute(X(1:N1,2:2:N2,:),[2,1,3]) ...
+ floor(filter(LiftFilter(:,1),1,X0(RightShift,:,:),...
X0(1,:,:)*LiftFilter(1,1),1));
end

X0 = X0 +floor( filter(LiftFilter(:,2),1,...
X1,X1(1,:,:)*LiftFilter(1,2),1)+0.5);

if rem(N2,2)
X1(M2,:,:) = [];
end

X(1:N1,1:N2,:) = permute([X0*ScaleFactor;X1/ScaleFactor],[2,1,3]);
end

N1 = M1;
N2 = M2;
end
else           % Inverse transform
for k = 1+Level:0
M1 = ceil(N1*pow2(k));
M2 = ceil(N2*pow2(k));

%%% Inverse transform along rows %%%
if M2 > 1
Q = ceil(M2/2);
RightShift = [2:Q,Q];
X1 = permute(X(1:M1,Q+1:M2,:)*ScaleFactor,[2,1,3]);

if rem(M2,2)
X1(Q,1,1) = 0;
end

% Undo lifting stages
X0 = permute(X(1:M1,1:Q,:)/ScaleFactor,[2,1,3]) ...
- floor(filter(LiftFilter(:,2),1,X1,X1(1,:,:)*LiftFilter(1,2),1)+0.5);
X1 = X1 - floor(filter(LiftFilter(:,1),1,X0(RightShift,:,:),...
X0(1,:,:)*LiftFilter(1,1),1));

if rem(M2,2)
X1(Q,:,:) = [];
end

X(1:M1,[1:2:M2,2:2:M2],:) = permute([X0;X1],[2,1,3]);
end

%%% Inverse transform along columns %%%
if M1 > 1
Q = ceil(M1/2);
RightShift = [2:Q,Q];
X1 = X(Q+1:M1,1:M2,:)*ScaleFactor;

if rem(M1,2)
X1(Q,1,1) = 0;
end

% Undo lifting stages
X0 = X(1:Q,1:M2,:)/ScaleFactor ...
- floor(filter(LiftFilter(:,2),1,X1,X1(1,:,:)*LiftFilter(1,2),1)+0.5);
X1 = X1 - floor(filter(LiftFilter(:,1),1,X0(RightShift,:,:),...
X0(1,:,:)*LiftFilter(1,1),1));

if rem(M1,2)
X1(Q,:,:) = [];
end

X([1:2:M1,2:2:M1],1:M2,:) = [X0;X1];
end
end
end

2. The test .m file:

clear all
close all
clc

n=100000;
maxt=zeros(1,n);
maxt2=zeros(1,n);
for it=1:n
X=floor(rand(8,8)*256);
X = X-128;
a = waveletlegall53(X,2);
maxt(it)=max(max(abs(a)));
if max(max(abs(a))) > 470
max(max(abs(a)))
end

end
[fr ind]=hist(maxt,length(unique(maxt)));
pr = length(find(maxt>512))/n
fr=fr/n;

figure()
plot(ind, fr)
grid on

Maxvalue = max(maxt)

• link I have answered my own question... – qiuhan1989 Jul 18 '16 at 12:58