testWavelet
Test discrete wavelet transform (incl. weight).
Contents
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
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('--')