%  Multirresolución mediante ENO-DB4 para funciones con ruido con las
%  modificaciones:
%            - Las extrapolaciones se realizan mediante AP de grado 1 y stencils de longitud 4.
%            - La detección de discontinuidades se realiza mediante el umbral VISU.
%        capenowavedb4TInoiseamcsc(x,f,fsin,pasos)
%        Entradas:
%            x     =  nodos iniciales
%            f     =  f(x) con ruido
%            fsin  =  función inicial sin ruido
%            pasos =  niveles que bajamos
%        Salidas:
%            plot de la reconstruccion (-) y de la funcion inicial sin
%            ruido (--) y con ruido (:)
%
% Autor: José Jaime Noguera.   +INFO: www.uv.es/anims/noguera

function  capenowavedb4TInoisevisu3(x,f,fsin,pasos)
[res,num]=size(f);

alfa(pasos+1,:)=f(:)';

c(1)=(1+sqrt(3))/(4*sqrt(2));
c(2)=(3+sqrt(3))/(4*sqrt(2));
c(3)=(3-sqrt(3))/(4*sqrt(2));
c(4)=(1-sqrt(3))/(4*sqrt(2));
h(1)=c(4); h(2)=-c(3); h(3)=c(2); h(4)=-c(1);
r=[c(2),c(1),c(4),c(3)]; 
d=[h(2),h(1),h(4),h(3)];
 
m=length(alfa(pasos +1,:));
s=zeros(pasos+1,num);

for l=1:m/2-1
    beta(pasos,l)=h(1)*alfa(pasos+1,2*l-1)+h(2)*alfa(pasos+1,2*l)+h(3)*alfa(pasos+1,2*l+1)+h(4)*alfa(pasos+1,2*l+2);
end
desviacion=median(abs(beta(pasos,1:m/2-1)))/0.6745; %MAD/0.6745
landa=desviacion*sqrt(2*log(num)); 

for j=pasos:-1:1
   m=m/2;   
   i=5;
   for l=1:m-1
      alfa(j,l)=c(1)*alfa(j+1,2*l-1)+c(2)*alfa(j+1,2*l)+c(3)*alfa(j+1,2*l+1)+c(4)*alfa(j+1,2*l+2);
      beta(j,l)=h(1)*alfa(j+1,2*l-1)+h(2)*alfa(j+1,2*l)+h(3)*alfa(j+1,2*l+1)+h(4)*alfa(j+1,2*l+2); 
   end    
           
   while(i<=m-3)    
        if (abs(beta(j,i))>=landa && abs(beta(j,i+1))<landa )
            %un solo stencil cruza la discontinuidad
            extpl=extrapola_lsa([0,1,2,3],[alfa(j,i-4),alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],1,4);
            beta(j,i)=(1/d(2))*(alfa(j+1,2*i-1)-r(2)*extpl-r(4)*alfa(j,i-1)-d(4)*beta(j,i-1));
            alfa(j,i)=(1/r(3))*(alfa(j+1,2*i+2)-r(1)*alfa(j,i+1)-d(1)*beta(j,i+1));
            s(j,i)=1;
        else
            if(abs(beta(j,i))>=landa && abs(beta(j,i+1))>=landa && abs(beta(j,i+2))<landa)
                %tenemos dos stencis que cruzan discontinuidad
                extpl=extrapola_lsa([0,1,2,3],[alfa(j,i-4),alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],1,4);
                %alf_left=(1/c(4))*(extpl-c(1)*alfa(j+1,2*i-1)- c(2)*alfa(j+1,2*i)-c(3)*alfa(j+1,2*i+1));
                beta(j,i)=(1/d(2))*(alfa(j+1,2*i-1)-r(2)*extpl-r(4)*alfa(j,i-1)-d(4)*beta(j,i-1));
                extpl2=extrapola_lsa([0,1,2,3],[alfa(j,i-4),alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],1,5);         
                beta(j,i+1)=(1/d(2))*(alfa(j+1,2*i+1)-r(2)*extpl2-r(4)*extpl-d(4)*beta(j,i));
                %lado derecho
                alfa(j,i+1)=(1/r(3))*(alfa(j+1,2*i+4)-r(1)*alfa(j,i+2)-d(1)*beta(j,i+2));
                alfa(j,i)=(1/r(3))*(alfa(j+1,2*i+2)-r(1)*alfa(j,i+1));
                s(j,i)=1;
                s(j,i+1)=1;
                i=i+1; 
            else
                %si no pues no hay discontinuidad y calcula lo normal (la beta ya esta de antes hecho)
                alfa(j,i)=c(1)*alfa(j+1,2*i-1)+c(2)*alfa(j+1,2*i)+c(3)*alfa(j+1,2*i+1)+c(4)*alfa(j+1,2*i+2);
            end    
        end
        i=i+1;       
   end
   %Fronteras
   alfa(j,m-1)=c(1)*alfa(j+1,2*(m-1)-1)+c(2)*alfa(j+1,2*(m-1))+c(3)*alfa(j+1,2*(m-1)+1)+c(4)*alfa(j+1,2*(m-1)+2);
   beta(j,m-1)=h(1)*alfa(j+1,2*(m-1)-1)+h(2)*alfa(j+1,2*(m-1))+h(3)*alfa(j+1,2*(m-1)+1)+h(4)*alfa(j+1,2*(m-1)+2);
    
   alfa(j,m-2)=c(1)*alfa(j+1,2*(m-2)-1)+c(2)*alfa(j+1,2*(m-2))+c(3)*alfa(j+1,2*(m-2)+1)+c(4)*alfa(j+1,2*(m-2)+2);
   beta(j,m-2)=h(1)*alfa(j+1,2*(m-2)-1)+h(2)*alfa(j+1,2*(m-2))+h(3)*alfa(j+1,2*(m-2)+1)+h(4)*alfa(j+1,2*(m-2)+2);
   
   beta(j,m)=h(1)*alfa(j+1,2*m-1)+h(2)*alfa(j+1,2*m)+h(3)*alfa(j+1,2*m)+h(4)*alfa(j+1,2*m); % esto lo extiendo para que salga mejor la frontera derecha
   alfa(j,m)=c(1)*alfa(j+1,2*m-1)+c(2)*alfa(j+1,2*m)+c(3)*alfa(j+1,2*m)+c(4)*alfa(j+1,2*m); % esto lo extiendo para que salga mejor la frontera derecha
   
   betaa(j,1)=h(1)*alfa(j+1,1)+h(2)*alfa(j+1,1)+h(3)*alfa(j+1,1)+h(4)*alfa(j+1,2); % esto lo estiendo para que la frontera izquerda salga mejor
   alfaa(j,1)=c(1)*alfa(j+1,1)+c(2)*alfa(j+1,1)+c(3)*alfa(j+1,1)+c(4)*alfa(j+1,2); % esto lo estiendo para que la frontera izquierda salga mejor
   %Fin_fronteras
end
 
 [s1,s2]=size(beta);
 beta=zeros(s1,s2);
 
 p=m/2;
 alfaf=zeros(pasos+1,num);
 alfaf(1,:)=alfa(1,:);

for j=1:pasos
     p=2*p; 
     %Fronteres izquierda
     alfaf(j+1,1)=r(2)*alfaf(j,1)  +d(2)*beta(j,1) + r(4)*alfaa(j,1) +d(4)*betaa(j,1);       %En la frontera.
     alfaf(j+1,2)=r(1)*alfaf(j,1)  +d(1)*beta(j,1) +  r(3)*alfaa(j,1)+d(3)*betaa(j,1);
     for i=2:p-1
        if(i>1 && s(j,i)==0 && s(j,i-1)==0 )
            alfaf(j+1,2*i-1)=r(2)*alfaf(j,i)+d(2)*beta(j,i)+ r(4)*alfaf(j,i-1)+d(4)*beta(j,i-1);
            alfaf(j+1,2*i)=r(1)*alfaf(j,i)+d(1)*beta(j,i)+ r(3)*alfaf(j,i-1)+d(3)*beta(j,i-1);
        end
        if(i>1 && s(j,i)==1 && s(j,i-1)==0)      
            extpl=extrapola_lsa([-1,0,1,2],[alfaf(j,i-4),alfaf(j,i-3),alfaf(j,i-2),alfaf(j,i-1)],1,3);
            alfaf(j+1,2*i-1)=r(2)*extpl+d(2)*beta(j,i)+ r(4)*alfaf(j,i-1)+d(4)*beta(j,i-1);
            alfaf(j+1,2*i)=r(1)*extpl+d(1)*beta(j,i)+ r(3)*alfaf(j,i-1)+d(3)*beta(j,i-1);
        end   
        if (i>1 && s(j,i)==0 && s(j,i-1)==1)  
            alfaf(j+1,2*i-1)=r(2)*alfaf(j,i)+d(2)*beta(j,i)+ r(4)*alfaf(j,i-1);
            alfaf(j+1,2*i)=r(1)*alfaf(j,i)+d(1)*beta(j,i)+ r(3)*alfaf(j,i-1);                              
        end
        if (i>1 && s(j,i)==1 && s(j,i-1)==1)
            ex1=extrapola_lsa([-1,0,1,2],[alfaf(j,i-5),alfaf(j,i-4),alfaf(j,i-3),alfaf(j,i-2)],1,3);
            ex2=extrapola_lsa([-1,0,1,2],[alfaf(j,i-5),alfaf(j,i-4),alfaf(j,i-3),alfaf(j,i-2)],1,4);
            alfaf(j+1,2*i-1)=r(2)*ex2+d(2)*beta(j,i)+ r(4)*ex1+d(4)*beta(j,i-1);
            alfaf(j+1,2*i)=r(1)*alfaf(j,i)+r(3)*alfaf(j,i-1);
        end
     end  
     %Frontera derecha
     alfaf(j+1,2*p-1)=r(2)*alfa(j,p)+d(2)*beta(j,p)+r(4)*alfa(j,p-1)+d(4)*beta(j,p-1);
     alfaf(j+1,2*p)=r(1)*alfa(j,p)+d(1)*beta(j,p)+r(3)*alfa(j,p-1)+d(3)*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 2 ]);
%set(gca,'Xtick',[-1 -0.5 0 0.5 1], 'FontSize',15)
%set(gca,'Ytick',[-1 -0.5 0 0.5 1 1.5 2], '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
%%%%%%%%%%%%%%%%
    
