Noise analysis
Part of the paper addresses the effect of the OGS algorithm on white Gaussian noise. This MATLAB program shows the calculations corresponding to that part of the paper. Specically, this program measure the standard deviation of unit variance Gaussian noise after OGS for group size K = (1,1) and group size K = (3,3). The calculations are also carried out for soft thresholding.
Contents
clear close all printme = @(txt) print('-dpdf', sprintf('figures/noise_analysis_%s', txt)); Fs = 8; % font size for titles. N1 = 500; % N1, N2 : size of 2D signal N2 = 600;
Output std versus threshold (soft-thresholding)
Find output std by computation:
L = 20; Tmax = 4; T = linspace(0,Tmax,L); s = zeros(L,1); for i = 1:L y = randn(N1,N2); a = soft(y, T(i)); s(i) = std(a(:)); % display signal before and after figure(1), clf imagesc([y; a], [-4 4]) axis image colormap(gray) colorbar title(sprintf('Soft thresholding, Thresh = %.2f', T(i))) drawnow % display histogram of values figure(2), clf del = 0.02; bins = -Tmax:del:Tmax; h = hist(a(:), bins); h = h/sum(h)/del; stairs(bins, h) title(sprintf('Theshold = %.2f', T(i))) ylim([0 1]) % printme(sprintf('hist_%d_%d_%d_%2d', K1*K2, K1, K2,round(100*lambda(i)))) end


Compare empirical std with formula (verify that formula is correct).
f = -2/sqrt(2*pi)* T .*exp(-T.^2/2) + 2 * (T.^2 + 1) .* qfunc(T); figure(2) clf plot(T, s,'o-', T, sqrt(f), '.-') legend('data','formula') xlabel('Threshold') title('standard deviation after soft thresholding') ylim([0 1.1]) slope = @(k) -sqrt(2) * gamma(k/2+1/2)./gamma(k/2); K = 1; % group size line([0 -1/slope(K)], [1 0], 'color', 'red') drawnow

Plot output std versus threshold using formula
T = 0:0.05:4; f = -2/sqrt(2*pi)* T .*exp(-T.^2/2) + 2 * (T.^2 + 1) .* qfunc(T); figure(3) clf plot(T, sqrt(f)); xlabel('Threshold') title('A) Soft thresholding') ylim([0 1.1]) ylabel('\sigma_x') set(gca,'ytick', 0:0.25:1); K = 1; % group size line([0 -1/slope(K)], [1 0], 'linestyle', '--') printme('soft_line')

Output histogram, group size K = 1
T = 1.3757; f1 = -2/sqrt(2*pi)* T .*exp(-T.^2/2); f2 = 2 * (T.^2) .* qfunc(T); f3 = 2 * qfunc(T); f = f1 + f2 + f3; S = sqrt(f) y = randn(N1,N2); a = soft(y, T); s = std(a(:)) % display signal before and after figure(1), clf imagesc([y; a], [-4 4]) axis image colormap(gray) colorbar title(sprintf('Soft thresholding, Thresh = %.2f', T)) drawnow % display histogram of values figure(2), clf del = 0.02; bins = -Tmax:del:Tmax; h = hist(a(:), bins); h = h/sum(h)/del; stairs(bins, h) title(sprintf('Theshold = %.2f', T)) ylim([0 1]) % printme(sprintf('hist_%d_%d_%d_%2d', K1*K2, K1, K2,round(100*lambda(i))))
S = 0.2500 s = 0.2518


Compute output histogram by the formula
a = -4:del:4; ac = abs(a) + T; p = 1/sqrt(2*pi) * exp(-ac.^2/2); P0 = 1 - 2*qfunc(T); % point mass at origin % check: sum(p)*del + P0 % should be 1 (area under pdf) figure(3) clf stem(0, P0, 'marker','^','MarkerFaceColor','black'); line(a, p) xlabel('a') txt = sprintf('T = %.2f', T); text(0.6, 0.9, txt, 'units', 'normalized') txt = sprintf('\\sigma_x = %.2f', S); text(0.6, 0.6, txt, 'units', 'normalized') text(0.1, 0.9, 'A) Soft thresholding', 'units', 'normalized') xlim([-3 3]) printme('soft_pdf')
ans = 1.0000

3-sigma
The threshold corresponding to 3-sigma
f = @(T) 2 * ((T.^2) + 1).*qfunc(T) - 2/sqrt(2*pi) * T .*exp(-T.^2/2); S = sqrt(f(3)) sqrt(f(3.36))
S = 0.0202 ans = 0.0100
Output std for group size K = (3,3)
This is precomputed and saved in a file as it takes a while to compute.
lambda = 0:0.02:0.6; L = length(lambda); K1 = 3; K2 = 3; Nit = 30; if exist('noise_analysis_K33.txt') s = load('noise_analysis_K33.txt'); else s = zeros(L,1); for i = 1:L y = randn(N1,N2); [a, cost] = ogshrink2(y, K1, K2, lambda(i), Nit); s(i) = std(a(:)); % display signal before and after figure(1), clf imagesc([y; a], [-4 4]) axis image colormap(gray), colorbar title(sprintf('OGS, lambda = %.2f', lambda(i))) drawnow % display histogram of values figure(2), clf bins = -3:del:3; h = hist(a(:), bins); h = h/sum(h)/del; stairs(bins, h) title(sprintf('lambda = %.2f', lambda(i))) end save noise_analysis_K33.txt s -ascii end
figure(2) clf plot(lambda, s); ylim([0 1.1]) xlim([0 lambda(end)]) xlabel('\lambda') ylabel('\sigma_x') title('B) OGS algorithm') txt = sprintf('Group size %d x %d', K1, K2); text(0.5, 0.8, txt, 'units', 'normalized', 'horizontalalignment', 'center'); % , 'fontsize', 14 set(gca,'ytick', 0:0.25:1); printme('ogs') line([0 -1/slope(K1*K2)], [1 0], 'linestyle', '--') printme('ogs_line')

Output histogram, group size 3x3. lambda = 0.265
lambda = 0.265; K1 = 3; K2 = 3; Nit = 30; y = randn(N1,N2); [a, cost] = ogshrink2(y, K1, K2, lambda, Nit); s = std(a(:)) % display signal before and after figure(1), clf imagesc([y; a], [-4 4]) axis image colormap(gray), colorbar title(sprintf('OGS, lambda = %.2f', lambda)) drawnow % display histogram of values del = 0.02; bins = -3:del:3; h = hist(a(:), bins); h = h/sum(h)/del; figure(2), clf plot(bins, h) txt = sprintf('\\lambda = %.2f', lambda); text(0.6, 0.9, txt, 'units', 'normalized') ylim([0 4]) xlabel('a') txt = sprintf('\\sigma_x = %.2f', s); text(0.6, 0.6, txt, 'units', 'normalized') text(0.1, 0.9, 'B) OGS algorithm', 'units', 'normalized') printme('ogs_pdf')
s = 0.2488


Output histogram, group size 3x3. lambda = 0.40
lambda = 0.4; K1 = 3; K2 = 3; Nit = 30; y = randn(N1,N2); [a, cost] = ogshrink2(y, K1, K2, lambda, Nit); s = std(a(:)) % display signal before and after figure(1), clf imagesc([y; a], [-4 4]) axis image colormap(gray), colorbar title(sprintf('OGS, lambda = %.2f', lambda)) drawnow % display histogram of values del = 0.03; bins = -3:del:3; h = hist(a(:), bins); h = h/sum(h)/del; P0 = max(h)*del; [tmp, k1] = max(h); h(k1) = h(k1+1); figure(2), clf stem(0, P0, 'marker','^','MarkerFaceColor','black'); line(bins, h) txt = sprintf('\\lambda = %.2f', lambda); text(0.6, 0.9, txt, 'units', 'normalized') xlabel('a') txt = sprintf('\\sigma_x = %.2f', s); text(0.6, 0.6, txt, 'units', 'normalized') text(0.1, 0.9, 'C) OGS algorithm', 'units', 'normalized') printme('ogs_pdf_2')
s = 0.0196

