function f=BAMS_alone(sign)

% This is stand-alone BAMS. 
% only dwtr.m and idwtr.m are required. 



%Take the filter H, in this case this is SYMMLET 4
filt = [ -0.07576571478934  -0.02963552764595  ...
          0.49761866763246   0.80373875180522  ...
          0.29785779560554  -0.09921954357694  ...
         -0.01260396726226   0.03222310060407]; 


sw = dwtr(sign, 5, filt);

n=length(sign);

coarsest=floor(log2(log(n)))+1;
J = log2(n);
swt = sw;
 %----------setting the parameters: mu----------------------
 indx=(2^(J-1)+1):2^J;
 niz=sw(indx);
 nizo=sort(niz); mm=length(nizo);
 low=floor(1/4 * mm); high=floor(3/4 * mm);
 pseudos = abs(nizo(high)-nizo(low))/1.5; %tukey's book
 mu = 1/pseudos^2;

 %----------setting pars: tau and epsilon (levelwise)-------
 finest_lev=J-1;
 coarsest = 5;
 shape=ones(1, finest_lev);
 tauu= sqrt((std(sign).^2 - 1/mu));
 for i = finest_lev:-1:coarsest
  eps(i) = 1 - 1/(i - coarsest +1)^1.5;
  tau(i) = tauu;
  
 %----------------------------------------------------------
 ind=(2^i+1):2^(i+1);
  a = sw(ind);
  b = 0 .* a;
  nn=length(a);
   for  j = 1:nn,
       [zz,  b(j)] = bayesrule(a(j), mu, tau(i), eps(i));
   end
   swt(ind)= b; 
 end
 

f=idwtr(swt, coarsest, filt);


%----------------------to be called------------------

function [me, bre]=bayesrule(d, mu, tau, eps)
%
%
%  Bayesrule: Calculates Marginal Distribution and Shrinkage Bayes Rule
%  (needed by BAMS)
%  Usage
%   [me, br]= bayesrule(d, mu, tau, eps)
%  Inputs
%    d      Wavelet Coefficient
%    mu     Parameter - concerning prior on the scale of noise
%    tau    Scale of the prior on location of d (signal)
%    eps    Weight of point mass at zero in the prior on signal
%  Outputs
%    me     value of the marginal distribution at d
%    br     value of the bayes rule at d
%
de = 1/2 * sqrt(2 * mu) * exp(- abs(d)*sqrt(2 * mu) );

m = (tau * exp(-abs(d)/tau) - 1/sqrt(2 * mu) * ...
   exp( - sqrt(2 * mu) * abs(d) ) )/(2 * tau^2 - 1/mu);
me = eps * de + (1-eps) * m;

if d >= 0,
    br= ( tau * (tau^2 - 1/(2*mu)) * d * exp(-d/tau) + ...
    tau^2/mu * ( exp(- sqrt(2*mu)*d) - exp(-d/tau)))/  ...
    ( (tau^2 - 1/(2*mu)) * (tau *  exp(-d/tau) -  ...
    1/sqrt(2*mu)*  exp(- sqrt(2*mu)*d) ) ); 
else
    br= ( tau * (tau^2 - 1/(2*mu)) * d * exp(d/tau) - ...
    tau^2/mu * ( exp(sqrt(2*mu)*d) - exp(d/tau)))/  ...
    ( (tau^2 - 1/(2*mu)) * (tau *  exp(d/tau) -  ...
    1/sqrt(2*mu)*  exp(sqrt(2*mu)*d) ) ); 
end

bre= (1-eps)*m*br / ( (1-eps)*m + eps * de );

%
% Copyright (c) 2000  B. Vidakovic and F. Ruggeri
%