function [wtc,vthr]=wdenmd(wt,L,sorh,approx,u);
%
% wavelet denoising, MinD algorithm, scale by scale (see Ranta et al, IEEE SPL 12(8) 2005, p. 561-564) 
%
% wt = wavelet coefficients vector (see wavdecn wavelet toolbox Matlab)
% L = bookeeping vector (see wavdecn wavelet toolbox Matlab))
% sorh =  s or h (soft or hard)
% approx = 1 or 2 (1 threshold the approximation; 2 keep the approximation)
% u = shape coefficient of the generalized gaussian (if u=0, it is then empirically estimated)
%
% wtc = denoised wavelet coefficient vector
% vthr = vector of thresholds, by scale

%initialisation
wtc=wt;
MM=length(L)-1;
vthr=zeros(1,MM);

if approx==1,
   echappr=1; % denoise the approximation
   indici=[1:L(MM+1)];
else
   echappr=2; % keep the approximation
   indici=[L(1)+1:L(MM+1)];
end

wtl=wt(indici);
meanwt=mean(wtl);
wtl=wtl-meanwt;

if u==0,
    %gaussiennes generalisées
    [u,DEV,MEAN]=ggfitr(wtl);
end
ffa=3*gamma(1/u)/u*(u*exp(1))^(1/u);
Fa=sqrt(ffa);
   
for j=echappr:MM, %pour chaque echelle > 1
   % calcul seuil
   indici=(sum(L(1:j-1))+1:sum(L(1:j)));
   wtl=wt(indici);
   meanwt=mean(wtl);
   wtl=wtl-meanwt;
   swt=sort(abs(wtl));
   Ma0=max(swt);
   lm=L(j);

   thr=Fa*sqrt(sum(swt(1:lm).^2)/lm);
   while thr<Ma0,
        Ma0=thr;
        [Ma1,lm1]=max(swt(find(swt<=Ma0)));
        thr=Fa*sqrt(sum(swt(1:lm1).^2)/lm);
   end
   vthr(j)=thr;
   if sorh=='h',       
        wtc(indici)=(wtl+meanwt).*(abs(wtl)>=thr);
   else
        wtc(indici)=(wtl+meanwt-sign(wtl)*thr).*(abs(wtl)>=thr);
   end
end