# testWavelet

Test discrete wavelet transform (incl. weight).

## Contents

- Syntax
- More About
- See Also
- Code

## Syntax

Start this routine in the folder `code` by

    init;
    testWavelet

## More About

- Discrete wavelet transform is described in ldwt.html.
- The discrete wavelet transform as well as its inverse (with and without weights) is described in setWavelet.html. This corresponds to the parameters `seti.W`, `seti.iW`, `seti.Wstar`, `seti.iWstar`, ..., `seti.wavWght`. (The generation of white noise is described in whiteNoiseCalc.html.)

## See Also

- ldwt.html
- setWavelet.html
- whiteNoiseCalc.html

## 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('--')