1
$\begingroup$

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)
    
$\endgroup$
  • $\begingroup$ link I have answered my own question... $\endgroup$ – qiuhan1989 Jul 18 '16 at 12:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.