function [xf, nround] = softthreshold(x, epsilon)
% Soft threshold filthrering
%           [xf, nround] = softthreshold(x, epsilon)
%
%       Input:
% x         Signal of length k * 2^J.
% epsilon   Treshold value. If epsilon is a scalar, the same value is used
%           for all frequency bands. If epsilon is a vector of length J+1,
%           the signal is assumed to be divided to J subbands.
%       Output:
% xf        The thresholded signal.
% nround    Number of signal values equal to 0.

nround = 0;
xf = zeros(size(x));
[nrows, ncols] = size(x);

% Specify number of subbands
if length( epsilon ) == 1
    % epsilon is a scalar
    J = 0;
else
    % J subbands
    J = length( epsilon ) - 1;
end

% Create a vector lenvec containing the lengths of subbands
lenvec = zeros(J+1,1);
lenvec(1) = ncols / 2^J;
% When k=2, exponent of 2 is J = J+2-2 = J+2-k
% When k=J+1, exponent of 2 is 1 = 1+J-J+1-1 = J+2-k
for k=2:(J+1)
    lenvec(k) = ncols / 2^(J+2-k);
end

previx = 1;
for k=1:(J+1)
    % Epsilon for this subband
    epsilonk = epsilon(k);
    % column indices for this subband
    ixs = previx:( previx - 1 + lenvec(k) );
    % update previx
    previx = previx + lenvec(k);
    % go through all rows and threshold subband given by ixs
    for n=1:nrows
        for m=ixs
            if abs(x(n,m)) <= epsilonk
                xf(n,m) = 0;
                nround = nround + 1;
            else
                xf(n,m) = sign( x(n,m) ) * ( abs(x(n,m)) - epsilonk );
            end
        end
    end
end