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