% Multirresolución mediante transformada wavelet ENO-DB4 igualando totos
% los coeficentes beta a 0. (Se hace al estilo Chan-Zhou )
%       capenowavedb4(x,f,pasos)
%       Entradas:
%            x     = nodos iniciales
%            f     = f(x)
%            pasos = niveles de multirresolución
%       Salidas:
%            plot de la reconstruccion (-), de la funcion inicial  (--) y
%            de la inicial menos la recosntrucción desplazada -0.1 en eje Y
%
% Autor: José Jaime Noguera.    +INFO: www.uv.es/anims/noguera

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

 %Tomamos estos valores, pero pueden modificarse
 a=2;
 bmayor=0.001;

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);

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

s=zeros(pasos+1,num);

for j=pasos:-1:1
   m=m/2;   
   
    % esto es para luego determinar las discontinuidades
    for l=1:m-3
        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 
    
    i=1;
    while(i<=m-3)    
         if(i>2 && i<m-1 && abs(beta(j,i))>bmayor && abs(beta(j,i-1))<=a*abs(beta(j,i-2))&& abs(beta(j,i))>a*abs(beta(j,i-1)) && abs(beta(j,i+1))<=a*abs(beta(j,i+2)))
            %un solo stencil cruza la discontinuidad
            extpl=extrapola([1,2],[alfa(j,i-2),alfa(j,i-1)],3);
            beta(j,i)=h(1)*alfa(j+1,2*i-1)+h(2)*alfa(j+1,2*i)+(h(3)/c(3))*(extpl-c(1)*alfa(j+1,2*i-1)-c(2)*alfa(j+1,2*i));
            alfa(j,i)=-(c(1)/h(1))*(h(3)*alfa(j+1,2*i+1)+h(4)*alfa(j+1,2*i+2))+c(3)*alfa(j+1,2*i+1)+c(4)*alfa(j+1,2*i+2);
           
            s(j,i)=1;
         else
             if( i>2 && i<m-1 && abs(beta(j,i))>bmayor && abs(beta(j,i-1))<=a*abs(beta(j,i-2)) && abs(beta(j,i))>a*abs(beta(j,i-1)) && abs(beta(j,i+1))>a*abs(beta(j,i+2)))
                %tenemos dos stencis que cruzan discontinuidad
                extpl=extrapola([1,2],[alfa(j,i-2),alfa(j,i-1)],3);
                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)=h(1)*alfa(j+1,2*i-1)+h(2)*alfa(j+1,2*i)+h(3)*alfa(j+1,2*i+1)+h(4)*alf_left;
                % alfa(j,i)= despues.
                extpl=extrapola([1,2],[alfa(j,i-2),alfa(j,i-1)],4);         
                beta(j,i+1)=h(1)*alfa(j+1,2*i+1)+h(2)*alf_left+(h(4)/c(4))*(extpl-c(1)*alfa(j+1,2*i+1)-c(2)*alf_left);
                alf_rigth=-(1/h(1))*(h(2)*alfa(j+1,2*i+2)+h(3)*alfa(j+1,2*i+3)+h(4)*alfa(j+1,2*i+4));
                alfa(j,i+1)=c(1)*alf_rigth+c(2)*alfa(j+1,2*i+2)+c(3)*alfa(j+1,2*i+3)+c(4)*alfa(j+1,2*i+4);
                alfa(j,i)=-(c(1)/h(1))*(h(3)* alf_rigth+h(4)*alfa(j+1,2*i+2))+ c(3)*alf_rigth+c(4)*alfa(j+1,2*i+2);
             
                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 estiendo para que salga mas bonita 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 estiendo para que salga mas bonita 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 mas bonita
    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 mas bonita
    % Fin-fronteras
 end
 
 [s1,s2]=size(beta);
 beta=zeros(s1,s2);

 p=m/2;
 alfaf=zeros(pasos+1,num);
 alfaf(1,:)=alfa(1,:);
 r=[c(2),c(1),c(4),c(3)];
 d=[h(2),h(1),h(4),h(3)];
 
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([1,2],[alfaf(j,i-2),alfaf(j,i-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([1,2],[alfaf(j,i-3),alfaf(j,i-2)],3);
            ex2=extrapola([1,2],[alfaf(j,i-3),alfaf(j,i-2)],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

disp('Norma infinito')
norm(alfaf(pasos+1,:)-f(:)',inf)
disp('Norma 2')
norm(alfaf(pasos+1,:)-f(:)',2)
    
hold on
plot(x, f,'--k','LineWidth',2)                                       %original 
plot(x,(f-alfaf(pasos+1,:))-(0.1)*ones(1,m_ini),':k','LineWidth',3)  %original-reconstruccion
plot(x,alfaf(pasos+1,:),'-k','LineWidth',3)                          %reconstruccion
axis([0 1 -0.2 1.1 ]);
set(gca,'Xtick',[0 .5 1], 'FontSize',15)
set(gca,'Ytick',[0 .5 1], 'FontSize',15)
box
hold off


    
