0

如标题所示，如果输入矩阵的值在范围内，我对 Wavelet LeGall 5/3 的系数范围（必须准确地使用这个滤波器）2D 变换（仅适用于 8*8 块）感到困惑从 0 到 255。

对于公式，链接在这里：Wavelet LeGall 5/3

这是我现在所做的：

  1. 所有值减 128（更容易计算低频值，见下文）；

  2. 在水平方向进行变换。这将在所有行中生成所有系数：前 4 个是低频，后 4 个是高频。很容易找到高频的范围是-255到+255（输入范围的两倍）。低频范围实际上是-192到+192（1.5*输入范围）。

  3. 在垂直方向进行变换。这将在垂直方向上做同样的事情。并且产生了四个块：LL（lowlow），LH（low high），HL，HH。很容易计算出 HH 的范围最大：-511 到 +511，LL 的范围是 1.5*1.5 = 2.25 倍输入范围（-128 到 +127）。

那么，问题来了。如果我再次为 LL 块做这个小波怎么办？理论上，LL块的HH（第二级系数）的范围应该是LL范围的4倍，即输入范围（-128到+127）的10倍，即-1280到+1270。

但是，我尝试了很多次随机计算，最大值从不超过-511到+511（见最后的代码）。我猜这是因为无法达到理论值，因为所有计算都是基于前一个计算的。但这个看似简单的问题，我很难从理论上证明。那么有人可以帮我出去吗？

我使用的代码（将两个文件放在一个目录中并在任何时候运行测试文件，但最大值不会超过 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