% Multirresolución mediante transformada wavelet DB4 con corte SURE y soft threshold.
%         capeliminawavesure(x,f,fsin,pasos)
%         Entradas:
%            x   =  nodos iniciales
%            f   =  f(x)
%            fsin=  funcion inicial sin ruido
%            pasos= niveles que bajamos
%         Salidas:
%            plot de la reconstruccion (-), de la funcion inicial con ruido (.) y de la funcion inicial sin ruido (--)rapolacio per a que les fronteres funcionen
%            ECM y SNR con varias combinaciones
%
% Autor: José Jaime Noguea.    +INFO: www.uv.es/anims/noguera

function  capeliminawavesure(x,f,fsin,pasos)
[res,num]=size(f);
alfa(pasos+1,:)=f(:)';

c1=(1+sqrt(3))/(4*sqrt(2));
c2=(3+sqrt(3))/(4*sqrt(2));
c3=(3-sqrt(3))/(4*sqrt(2));
c4=(1-sqrt(3))/(4*sqrt(2));
h1=c4; h2=-c3; h3=c2; h4=-c1;

m=length(alfa(pasos +1,:));
m_ini=m;

for j=pasos:-1:1
    m=m/2;   
    for i=1:m-1             
         beta(j,i)=h1*alfa(j+1,2*i-1)+h2*alfa(j+1,2*i)+h3*alfa(j+1,2*i+1)+h4*alfa(j+1,2*i+2);
         alfa(j,i)=c1*alfa(j+1,2*i-1)+c2*alfa(j+1,2*i)+c3*alfa(j+1,2*i+1)+c4*alfa(j+1,2*i+2);     
    end     
    beta(j,m)=h1*alfa(j+1,2*m-1)+h2*alfa(j+1,2*m)+h3*alfa(j+1,2*m)+h4*alfa(j+1,2*m); % esto lo extiendo para que salga mejor la frontera derecha
    alfa(j,m)=c1*alfa(j+1,2*m-1)+c2*alfa(j+1,2*m)+c3*alfa(j+1,2*m)+c4*alfa(j+1,2*m); % esto lo extiendo para que salga mejor la frontera derecha
    betaa(j,1)=h1*alfa(j+1,1)+h2*alfa(j+1,1)+h3*alfa(j+1,1)+h4*alfa(j+1,2); % esto lo extiendo para que la frontera izquerda salga mejor
    alfaa(j,1)=c1*alfa(j+1,1)+c2*alfa(j+1,1)+c3*alfa(j+1,1)+c4*alfa(j+1,2); % esto lo extiendo para que la frontera izquierda salga mejor
end
mfin=m;
[s1,s2]=size(beta);

 m_aqui=m_ini;
 
 %%%%%%%SURE-SHRINK %%%%%%%%%%%%%%
 desviacion=median(abs(beta(pasos,1:m-1)))/0.6745; %MAD/0.6745
 landau=desviacion*sqrt(2*log(num));
 m=m_ini;
 for j=pasos:-1:1
     m=m/2;
     for k=1:m-1
        contador=0;
        suma=0;
        for i=1:m-1
            if( abs(beta(j,i)) < abs(beta(j,k)))
               contador=contador+1;
            end
            suma=suma+(min(abs(beta(j,i)),abs(beta(j,k))))^2;
        end
        if( abs(beta(j,k))>=landau )
            sure(k)=100000;
        else     
            sure(k)=m-1-2*contador+suma;
        end   
     end
     [minim,pos]=min(sure(1:m-1));
     landa(j)=abs(beta(j,pos));
 end
 for j=s1:-1:1
      m_aqui=m_aqui/2;      
      for i=1:m_aqui
          % HARD
          %if(abs(beta(j,i))<=landa)
          %   beta(j,i)=0;                                                        
          %end
          % SOFT
          beta(j,i)=sign(beta(j,i))*max(0,abs(beta(j,i))-landa(j));
      end
          betaa(j,1)=sign(betaa(j,1))*max(0,abs(betaa(j,1))-landa(j)); % frontera
    end      
 %FIN VISU-SHRINK
 
 p=mfin/2;
 alfaf=zeros(pasos+1,num);
 alfaf(1,:)=alfa(1,:);
 
 for j=1:pasos
     p=2*p; 
     %Fronteres izquierda
     alfaf(j+1,1)=c1*alfaf(j,1)  +h1*beta(j,1) + c3*alfaa(j,1) +h3*betaa(j,1);       
     alfaf(j+1,2)=c2*alfaf(j,1)  +h2*beta(j,1) +  c4*alfaa(j,1)+h4*betaa(j,1);
     for i=2:p-1       
         alfaf(j+1,2*i-1)=c1*alfaf(j,i)+h1*beta(j,i)+c3*alfaf(j,i-1)+h3*beta(j,i-1);  
         alfaf(j+1,2*i)= c2*alfaf(j,i)+h2*beta(j,i)+c4*alfaf(j,i-1)+h4*beta(j,i-1);
     end  
     %Frontera derecha
      alfaf(j+1,2*p-1)=c1*alfa(j,p)+h1*beta(j,p)+c3*alfa(j,p-1)+h3*beta(j,p-1);
      alfaf(j+1,2*p)  =c2*alfa(j,p)+h2*beta(j,p)+c4*alfa(j,p-1)+h4*beta(j,p-1);
 end

 fr=alfaf(pasos+1,:);
   
    disp('Entre fr y fsin')
    error=mse(fsin,fr);
    error
 
    disp('Entre fr y fn')
    error=mse(f,fr);
    error
    
    disp('Entre fr y fn')
    error=calculosnr(fr,f);
    error
    
    disp('Entre fin y fr')
    error=calculosnr(fsin,fr);
    error
   
hold on
   
%axis([-1,1,-0.1,1.5])
plot(x, f,'k.','LineWidth',5,'MarkerSize',7)   
plot(x, fsin,'--k','LineWidth',2) 
plot(x,fr,'-k','LineWidth',3);

%axis([-1 1 -1.2 1.2 ]);
%set(gca,'Xtick',[-1 -0.5 0 0.5 1], 'FontSize',15)
%set(gca,'Ytick',[-1 -0.5 0 0.5 1], 'FontSize',15)

%axis([0 1 -0.1 1.1 ]);
%set(gca,'Xtick',[0 .5 1], 'FontSize',15)
%set(gca,'Ytick',[0 .5 1], 'FontSize',15)
box
hold off



    
