%%FUNCION QUE REALIZA EL ALGORITMO DE HARD WAVELET THRESHOLD
function x_est=hard_shrink(y,L,wfilt);
% L is number of Wavelet steps
% wfilt is for waveletfilter filter
[C S]=wavedec2(y,L,wfilt); %L depth 2d wavelet transform
n=1; %not hyperbolic
k=0; %soft

level_size=prod(S(end-1,:));

HH=C(3*level_size+1:4*level_size); %get detail coefficients of wavelet at current Level
%error estimate using current HH
sigma=median(abs(HH(:)))/0.6745; %sqrt(sum((HH-mean(HH)).^2)/(length(HH)-1));
% %error estimate using current HH
lambda=sqrt(2*log(prod(S(end,:))));
    
%determine threshold
threshold=0.5*lambda*sigma;
%determine threshold

offset=prod(S(1,:));

%wavelet thresholding
for(level=1:L)
    
    level_size=prod(S(level+1,:));
    
    %----------
    HH=C(offset+2*level_size+1:offset+3*level_size); %get detail coefficients of wavelet at current Level
    HH=sign(HH).*(thresholder(HH.^n,threshold^n,k)).^(1/n);
    
    %-----------
    HL=C(offset+1*level_size+1:offset+2*level_size);
    HL=sign(HL).*(thresholder(HL.^n,threshold^n,k)).^(1/n);
    
    %-------------
    LH=C(offset+1:offset+1*level_size);
    LH=sign(LH).*(thresholder(LH.^n,threshold^n,k)).^(1/n);
    
    %--------------
    
    C(offset+1:offset+1*level_size)=LH;
    C(offset+1*level_size+1:offset+2*level_size)=HL;    
    C(offset+2*level_size+1:offset+3*level_size)=HH;
    
    offset=offset+level_size*3;
    
    clear HH HL LH
    
end
%wavelet thresholding

 x_est=waverec2(C,S,wfilt);
