testWavelet

Test discrete wavelet transform (incl. weight).

Contents

Syntax

Start this routine in the folder code by

init;
testWavelet

More About

See Also

Code

Initialization

init;

seti.dim = 2;

seti.nCD = 128;
seti.rCD = 0.5;

seti.nInv = 64;
seti.gscale = 0;

% -- setWavelet
seti.useWavelet = 1; % 1: yes; 0: no (wavelets...)
seti.wavWeight = 1;
seti.genWhiteNoiseNew = 0; % generate new white noise (even when it is stored)
seti.adaptWeightsMan = 1; % adapt weights manually (only available in 2D)
seti.wavelet = 'cdf97'; % wavelet type (documentation in ldwt.m): haar, legall53, cdf97
seti.extension = 'symmetric'; % symmetric, zeropadded
seti.smin = 8; % mimimal size of a wavelet level (in 1 dimension)
seti.samples = 10000; % number of samples to compute the white noise

seti = setGrid(seti);
seti = setReshapeVecMat(seti);
seti = setGridScale(seti);

seti = setWavelet(seti); % calls whiteNoiseCalc either

Set up Wavelets with or without Weight

% Notation:
%
% W = seti.W;
% iW = seti.iW;
% Wstar = seti.Wstar;
% iWstar = seti.iWstar;

% w = seti.wavWght; % weights
% If you want to test without weigths use the following and set w = 1
% seti.wW = @(x) w.*seti.W(x);
% seti.wiW = @(x) seti.iW(1./w.*x);
% seti.wWstar = @(x) seti.Wstar(w.*x);
% seti.wiWstar = @(x) 1./w.*seti.iWstar(x);

Tests

wW = seti.wW;
wiW = seti.wiW;
wWstar = seti.wWstar;
wiWstar = seti.wiWstar;

nnROI = seti.nROI^seti.dim;

n = 100; % iteration number of random vectors

% scalar product
S = @(x,y) x'*y;

% identities
% Id1 = @(x) iW(W(x));
% Id2 = @(x) W(iW(x));
% Id3 = @(x) iWstar(Wstar(x));
% Id4 = @(x) Wstar(iWstar(x));

Id1 = @(x) wiW(wW(x));
Id2 = @(x) wW(wiW(x));
Id3 = @(x) wiWstar(wWstar(x));
Id4 = @(x) wWstar(wiWstar(x));

disp('-- test wavelet transformation --')
disp('identities:')
disp('1) Id1( . ) = wiW(wW( . ))         | 2) Id2( . ) = wW(wiW( . ))')
disp('3) Id3( . ) = wiWstar(wWstar( . )) | 4) Id4( . ) = wWstar(wiWstar( . ))')
disp('adjoints:')
disp('5) <wW(x),y> = <x,wWstar(y)>       | 6) <wiW(x),y> = <x,wiWstar(y)>')

D = zeros(n,6); % D: difference...
for i = 1:n

    x = rand(nnROI,1)+1i*rand(nnROI,1);
    y = rand(nnROI,1)+1i*rand(nnROI,1);

    % identities
    D(i,1) = max(abs(transpose(x-Id1(x))));
    D(i,2) = max(abs(transpose(x-Id2(x))));
    D(i,3) = max(abs(transpose(x-Id3(x))));
    D(i,4) = max(abs(transpose(x-Id4(x))));

    % test adjoint via scalar product
    D(i,5) = abs( S(wW(x),y)-S(x,wWstar(y)) ); % <Wx,y> = <x,W* y> (<w.*W(x),y> = <x,W*(w.*y)>
    D(i,6) = abs( S(wiW(x),y)-S(x,wiWstar(y)) ); % <iW(x),y> = <x,iW*(y)> (<iW(1./w.*x),y> = <x,1./w.*iW*(y)>)

end
disp('Maximal absolute differences:')
max(D)

disp('--')
disp('Additional: test identities with unity vectors')
d = zeros(nnROI,1);
for i = 1:nnROI
    if i == 1
        fprintf('i = %g | ',i);
    elseif floor(i/100) == i/100
        fprintf('%g | ',i);
    elseif floor(i/1001) == i/1001
        fprintf('\n');
    end
    x = zeros(nnROI,1);
    x(i) = 1;
    d(i,1) = max(abs(transpose(x-Id1(x))));
    d(i,2) = max(abs(transpose(x-Id2(x))));
    d(i,3) = max(abs(transpose(x-Id3(x))));
    d(i,4) = max(abs(transpose(x-Id4(x))));
end
disp('Maximal absolute differences:')
max(d)
disp('--')