function w = whiteNoiseCalc(seti,W,G,iG,nCD,nROI) % white noise calculation % two cases: temp. definitions in this file (via calling) % 1) nCD = seti.nCD, nROI = seti.nROI % 2) nCD = seti.nCDinv, nROI = seti.nInv % figure 4: weights from white noise % figure 5: weights from white noise manually adapted % figure 6: convergence of white noise plotWhiteNoise = 0; % plot figures 2 and 3 plotConv = 0; % to plot figure 4 samples = seti.samples; %numROI = seti.nROI.^seti.dim; % number of entries in ROI %numInv = seti.nInv.^seti.dim; % number of entries in nInv %n = zeros(numROI,1); %n = zeros(numInv,1); %numROI = nCD.^seti.dim; numROI = nROI.^seti.dim; n = zeros(numROI,1); %wn = n; % weight from white noise %w = n; % weight from white noise maybe manually adapted state = rng(); rng(0) %fileWhiteNoise = sprintf('mat/whiteNoise/wav_%s_ext_%s_smin_%d_nCD_%d_nROI_%d_dim_%d_samples_%d.mat',seti.wavelet,seti.extension,seti.smin,seti.nCD,seti.nROI,seti.dim,seti.samples); %fileWhiteNoise = sprintf('mat/whiteNoise/wav_%s_ext_%s_smin_%d_nCDinv_%d_nInv_%d_dim_%d_samples_%d.mat',seti.wavelet,seti.extension,seti.smin,seti.nCDinv,seti.nInv,seti.dim,seti.samples); fileWhiteNoise = sprintf('inmat/whiteNoise/wav_%s_ext_%s_smin_%d_nCD_%d_nROI_%d_dim_%d_samples_%d.mat',seti.wavelet,seti.extension,seti.smin,nCD,nROI,seti.dim,seti.samples); if (exist(fileWhiteNoise, 'file') == 2 && seti.genWhiteNoiseNew ~= 1) %file with name exists disp('Loading white noise matrix.'); load(fileWhiteNoise); % contains white noise matrix n else disp('Generating white noise...'); conv1vec = zeros(samples,1); conv2vec = zeros(samples,1); for i = 1:samples %noise = randn(size(n))+1i*randn(size(n)); noise = randn(size(n)); nBefore = n; %noise = noise/normws(noise,g)*normws(n,g)*delta; % optional... %n = n + abs(ldwt(noise,'forward',wavelet,extension,smin)).^2; %W = @(x) iG(ldwt(G(x),'forward',wavelet,extension,smin)); n = n + abs(W(noise)).^2; % W = @(x) iG(ldwt(G(x),'forward',wavelet,extension,smin)); %n = n + abs(W(n)).^2; %n = n + abs(ldwt(randn(sizeq))).^2; %original % -- convergence... experimental nDiff = n-nBefore; conv1 = norm(nDiff)/norm(n); conv2 = max(nDiff(:))/norm(n); % maximal entry in matrix / norm conv1vec(i) = conv1; conv2vec(i) = conv2; breakWithConv = 0; % experimental if breakWithConv fprintf('%d | %g | %g\n',i,conv1,conv2); if (conv1 < 1E-4 && conv2 < 1E-5) disp('convergence criterium is true. break.') break; end end % -- end if plotConv figure(6); xi = 100:samples; plot(xi,conv1vec(xi),xi,conv2vec(xi)); end n = 1/samples.*n; %white noise matrix % now: vector! fprintf('finished.\n'); rng(state); disp('Save white noise matrix.'); %now: vector! try mkdir inmat/whiteNoise/ catch end save(fileWhiteNoise,'n','-mat') end disp('calculating weight matrix...') % now: vector! wn = 1./n; %the resulting weights are % wn: weights from noise % now: use matrix... wn = G(wn); % matrix! w = zeros(size(wn)); %previously: nA = nCD... but here restricted to ROI % -- adapt weights manually (only in 2D) %nA = seti.nROI; %nA = seti.nInv; nA = nROI; %nA: n of Area in 1 dimension (previously: in first step nA = nCD; now: nA = nROI) if (nA < 2*seti.smin && seti.adaptWeightsMan) disp('ROI is too small to use manually adapted weights. So they are not used.') end if (nA >= 2*seti.smin && seti.adaptWeightsMan && ((strcmp(seti.wavelet,'legall53')) || (strcmp(seti.wavelet,'cdf97')))) % adapt weights manually AND use (wavelet legall53 OR wavelet cdf97) % ----- adapt weights manually mm = mean(mean(wn)); iwdepth = 0; %while nA >= smin while nA >= 2*seti.smin % to avoid problems with boundary %for iwdepth = 1:wdepth iwdepth = iwdepth+1; %i1 = 1:floor(nA/2); %i2 = (floor(nA/2)+1:nA); i1 = 1:floor(nA/2); i2 = (floor(nA/2)+1:nA); wn1 = wn(i1,i1); %top left wn2 = wn(i1,i2); %top right wn3 = wn(i2,i1); %bottom left wn4 = wn(i2,i2); %bottom right %automatic weight m = zeros(4,1); m(1) = mean(mean(wn1)); m(2) = mean(mean(wn2)); m(3) = mean(mean(wn3)); % m(3) = m(2) m(4) = mean(mean(wn4)); %wm1 = 2^iwdepth*m1/m2 %wm2 = 2^iwdepth*m2/m2 %wm3 = 2^iwdepth*m3/m2 wm = m/mm; % compared with complete w-Matrix (mm = mean(mean(w))) %wm = wm.^20; % adapt weight manually % wm = wm; %wm: weight manual %wm1 = 2^1; %wm2 = 2^2; %wm3 = 2^3; %multiplication: manual weight and noise weight w1 = wm(1).*wn1; w2 = wm(2).*wn2; w3 = wm(3).*wn3; w4 = wm(4).*wn4; w(1:nA,1:nA) = [w1, w2; w3, w4]; % resulting weight matrix % next iteration will overwrite w1 (this is intended) % (except there is no next iteration step) %plot in every iteration (progress can be seen) %figure(2+iwdepth); %imagesc(w); colorbar; %axis xy; nA = floor(nA/2); % new nA for next iteration end % ----- else %case: wavelet = 'haar'; OR no manually adapted weights %do not use manually adapted weights (use just the weights from whiteNoise) disp('else...'); w = wn; end % ----- if plotWhiteNoise figure(4); imagesc(wn); colorbar; axis xy; surf(wn);colorbar;shading interp; title('weights from white noise') figure(5); imagesc(w); colorbar; axis xy; surf(w);colorbar;shading interp; title('weights: manually adapted'); error('stop') end % make a vector of the matrix w = iG(w); %size(w) end