function [wtc,vthr]=wdenmdg(wt,L,sorh,approx,u);
%
% wavelet denoising, MinD algorithm, global (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

wtc=wt;

MM=length(L)-1;
 
if approx==1,
   indici=[1:L(MM+1)]; % denoise the approximation
else
   indici=[L(1)+1:L(MM+1)]; % keep the approximation
end
wtc(indici)=0;
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);


lm=length(wtl);
swt=sort(abs(wtl));
Ma0=max(swt);
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

wtl=[zeros(1,L(MM+1)-length(wtl)) wtl];
if sorh=='h',       
    wtc(abs(wtl)>=thr)=wtl(abs(wtl)>=thr)+meanwt;
else
    wtc(abs(wtl)>=thr)=wtl(abs(wtl)>=thr)+meanwt-sign(wtl(abs(wtl)>=thr))*thr;
end

vthr=thr*ones(1,MM);
vthr(1)=abs(approx-2)*thr;
% if approx==2,
%     wtc(1:L(1))=wt(1:L(1))+meanwt;
% end