Как и в заголовке, меня смущает диапазон коэффициентов вейвлета ЛеГалла 5/3 (должен быть именно этот фильтр) 2D-преобразование (только для блока 8 * 8), если значение входной матрицы находится в пределах диапазона от 0-255.
Ссылка на формулы здесь: Wavelet LeGall 5/ 3
Вот что я сделал на данный момент:
Минус 128 для всех значений (проще вычислить низкочастотные значения, см. далее);
Сделайте преобразование в горизонтальном направлении. Это сгенерирует все коэффициенты во всех строках: первые 4 — низкочастотные, а последние 4 — высокочастотные. Легко найти диапазон для высоких частот от -255 до +255 (удвоить диапазон ввода). И диапазон для низких частот на самом деле составляет от -192 до +192 (1,5 * диапазон ввода).
Сделайте преобразование в вертикальном направлении. Это будет делать то же самое в вертикальном направлении. Всего сгенерировано четыре блока: LL (lowlow), LH (low high), HL, HH. Легко рассчитать диапазон для HH, который является самым большим: от -511 до +511, а диапазон для LL составляет 1,5 * 1,5 = 2,25 раза от входного диапазона (от -128 до +127).
Тогда возникает вопрос. Что, если я повторю этот вейвлет для блока LL? Теоретически диапазон HH (коэффициенты 2-го уровня) блока LL должен в 4 раза превышать диапазон LL, который становится 10-кратным диапазоном ввода (от -128 до +127), то есть от -1280 до +1270.
Однако я много раз пробовал случайные вычисления, максимальное значение никогда не превышало от -511 до +511 (см. код в конце). Я думаю, это потому, что теоретические значения не могут быть достигнуты, потому что все расчеты основаны на предыдущем. Но этот, казалось бы, простой вопрос мне трудно доказать теоретически. Так может кто-нибудь помочь мне выбраться, пожалуйста?
Код, который я использовал (поместите два файла в один каталог и запускайте тестовый файл в любое время, но только максимальное значение не будет превышать 512...):
Функция вейвлетлегалла53:
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Тестовый файл .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)
