function [c,cnz]=threshold(d,dn,kind,div)
% [c,cnz]=threshold(d,dn,kind,div)
%
% Thresholds wavelet or scaling coefficients
% by selection WITHIN EACH SCALE (see Mallat p. 437 ff)
%
% INPUT:
%
% d       Coefficient cell array
% dn      Length of the cell array
% kind    'hard' or 'soft' [default]
%         (May need to look into halving T for soft thresholding.)
% div     Divides threshold level, if at all needed [default: 1]
%         Higher values means lower threshold [div=999: no thresholding] 
%
% OUTPUT:
%
% c       Thresholded coefficients
% cnz     Number of nonzero coefficients in every cell
%
% SEE ALSO: THRESHOLD2
%
% Last modified by fjsimons-at-alum.mit.edu, 09/19/2007

defval('kind','soft')
defval('div',1)

if div==999
  c=d;
  cnz=[];
else
  % Mallat Equation 10.53, Donoho and Johnstone 1994

  for index=1:length(dn)    
    % I think 25.04.2006 that this is good. It is SCALE-DEPENDENT
    % THRESHOLDING where the VARIABILITY at every scale is determined by
    % the MEDIAN absolute deviation from the MEDIAN of the coefficients
    % at every scale. See Johnstone and Silverman.
    % Calculate location estimate
    % Calculate different scale estimate for coefficients at ALL
    % resolution scales 
    % Use median absolute deviation from the median, DJ p 446
    % But they do this only at the smallest scale;
    % We do it at all scales
    % Also, not all coefficients should be shrunk - only the number of
    % vanishing moments... those have an expectation zero
    
    % On second thought, it's sort of keeping with the philosophy, to use
    % one threshold across ALL scales. Should provide this here, too.
    dloc(index)=median(d{index});
    dscal(index)=median(abs(d{index}-dloc(index)))/0.6745;
    T(index)=dscal(index)*sqrt(2*log(dn(index)))/div;
    switch kind
     case 'hard'
      dist=abs(d{index})>T(index);
      c{index}=d{index}.*dist;
      if index==1 ; disp('Hard thresholding') ; end
     case 'soft'
      dist=abs(d{index})-T(index);
      dist=(dist+abs(dist))/2;
      c{index}=sign(d{index}).*dist;
      if index==1 ; disp('Soft thresholding') ; end
    end
    % Return effective number of (nonzero) coefficients
    cnz(index)=sum(~~dist);
  end
end