0

タイトルにあるように、入力行列の値が範囲内にある場合、Wavelet LeGall 5/3 (正確にこのフィルターである必要があります) 2D 変換 (8*8 ブロックの場合のみ) の係数の範囲について混乱しています。 0 から 255 まで。

数式については、リンクは次のとおりです: Wavelet LeGall 5/3

これが私が今やったことです:

  1. すべての値に対してマイナス 128 (低頻度の値を計算する方が簡単です。後で参照してください)。

  2. 水平方向に変換します。これにより、すべての行のすべての係数が生成されます。最初の 4 つは低周波で、最後の 4 つは高周波です。高周波の範囲は、-255 ~ +255 (入力範囲の 2 倍) であることが簡単にわかります。また、低周波の範囲は実際には -192 ~ +192 (入力範囲の 1.5 倍) です。

  3. 垂直方向に変換します。これは、垂直方向にも同じことを行います。生成されるブロックは、LL (lowlow)、LH (low high)、HL、HH の 4 つです。HH の範囲を計算するのは簡単です: -511 から +511 まで、LL の範囲は入力範囲の 1.5*1.5 = 2.25 倍 (-128 から +127) です。

それでは、ここで質問です。LL ブロックに対してこのウェーブレットをもう一度実行するとどうなりますか? 理論的には、LL ブロックの HH (第 2 レベル係数) の範囲は、入力範囲 (-128 ~ +127) である -1280 ~ +1270 の 10 倍になる LL 範囲の 4 倍になるはずです。

ただし、ランダムな計算を何度も試みましたが、最大値は -511 から +511 を超えることはありません (最後のコードを参照)。前回の計算に基づいて計算しているため、理論値に到達できないためだと思います。しかし、この一見簡単そうな問題を理論的に証明するのは私には困難です。誰か私が脱出するのを手伝ってくれませんか?

私が使用したコード (1 つのディレクトリに 2 つのファイルを配置し、必要なときにいつでもテスト ファイルを実行しますが、最大値だけが 512 を超えることはありません...):

  1. 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. テスト .m ファイル:

    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)
    
4

0 に答える 0