% Transformada wavelet ENO-DB6 usando la TRANSFORMACION INVERSA. Igualamos todos los coeficentes beta a 0.
%       capenowavedb6TI(x,f,pasos)
%         Entradas:
%            x   =  nodos iniciales
%            f   =  f(x)
%            pasos= niveles que bajamos
%         Salidas:
%            plot de la reconstruccion (-) y de la funcion inicial  (--)
%
% Autor: José Jaime Noguera.   +INFO: www.uv.es/anims/noguera

function  capenowavedb6TI(x,f,pasos)
[res,num]=size(f);
   %Tomamos por defecto estos valores, aunque puede modificarse
   a=2;
   bmayor=0.0001;

alfa(pasos+1,:)=f(:)';
beta(1:pasos,1:length(x))=zeros(pasos,length(x));

[c,h]=wfilters('db3','d');

p=3;
for i=1:p
    r(2*i)=c(2*i-1);
    d(2*i)=h(2*i-1);
    r(2*i-1)=c(2*i);
    d(2*i-1)=h(2*i);
end    

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

for j=pasos:-1:1
   m=m/2; 
   
    % se utiliza despues, pero se sobreescribe lo que se va calculando
    for l=1:m-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)+h(5)*alfa(j+1,2*l+3)+h(6)*alfa(j+1,2*l+4);    
         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)+c(5)*alfa(j+1,2*l+3)+c(6)*alfa(j+1,2*l+4);   
    end
   
    i=1;        
    while(i<=m-6)    
       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+2))<=a*abs(beta(j,i+3)))
         
            %dos stencils stencils cruzan la discontinuidad
             %lado izquierdo
            alfa(j,i)=extrapola([1,2,3],[alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],4);
                 
            beta(j,i)=(1/d(2))*(alfa(j+1,2*i-1)-r(2)*alfa(j,i)-r(4)*alfa(j,i-1)-d(4)*beta(j,i-1)-r(6)*alfa(j,i-2)-d(6)*beta(j,i-2));
  
            alfa(j,i+1)=extrapola([1,2,3],[alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],5);
            beta(j,i+1)=(1/d(2))*(alfa(j+1,2*i+1)-r(2)*alfa(j,i+1)-r(4)*alfa(j,i)-d(4)*beta(j,i)-r(6)*alfa(j,i-1)-d(6)*beta(j,i-1));
             %lado derecho
            alfa(j,i+1)=(1/r(5))*(alfa(j+1,2*i+6)-r(1)*alfa(j,i+3)-d(1)*beta(j,i+3)-r(3)*alfa(j,i+2)-d(3)*beta(j,i+2)); 
            alfa(j,i)=(1/r(5))*(alfa(j+1,2*i+4)-r(1)*alfa(j,i+2)-d(1)*beta(j,i+2)-r(3)*alfa(j,i+1));
           
            s(j,i)=1;
            s(j,i+1)=1;
            i=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+2))>a*abs(beta(j,i+3)))
               %tenemos tres stencils que cruzan discontinuidad
               %lado izquierdo
               alfa(j,i)=extrapola([1,2,3],[alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],4);
               beta(j,i)=(1/d(2))*(alfa(j+1,2*i-1)-r(2)*alfa(j,i)-r(4)*alfa(j,i-1)-d(4)*beta(j,i-1)-r(6)*alfa(j,i-2)-d(6)*beta(j,i-2));
               alfa(j,i+1)=extrapola([1,2,3],[alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],5);
               beta(j,i+1)=(1/d(2))*(alfa(j+1,2*i+1)-r(2)*alfa(j,i+1)-r(4)*alfa(j,i)-d(4)*beta(j,i)-r(6)*alfa(j,i-1)-d(6)*beta(j,i-1));
               alfa(j,i+2)=extrapola([1,2,3],[alfa(j,i-3),alfa(j,i-2),alfa(j,i-1)],6);
               beta(j,i+2)=(1/d(2))*(alfa(j+1,2*i+3)-r(2)*alfa(j,i+2)-r(4)*alfa(j,i+1)-d(4)*beta(j,i+1)-r(6)*alfa(j,i)-d(6)*beta(j,i));
   
               %lado derecho
               alfa(j,i+2)=(1/r(5))*(alfa(j+1,2*i+8)-r(1)*alfa(j,i+4)-d(1)*beta(j,i+4)-r(3)*alfa(j,i+3)-d(3)*beta(j,i+3)); 
               alfa(j,i+1)=(1/r(5))*(alfa(j+1,2*i+6)-r(1)*alfa(j,i+3)-d(1)*beta(j,i+3)-r(3)*alfa(j,i+2));
               alfa(j,i)=(1/r(5))*(alfa(j+1,2*i+4)-r(1)*alfa(j,i+2)-r(3)*alfa(j,i+1));
             
               s(j,i)=1;
               s(j,i+1)=1;
               s(j,i+2)=1;
               i=i+2; 
             
            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)+c(5)*alfa(j+1,2*i+3)+c(6)*alfa(j+1,2*i+4);

            end    
       end
       i=i+1;       
    end
    %Fronteras
    alfa(j,m-5)=c(1)*alfa(j+1,2*(m-5)-1)+c(2)*alfa(j+1,2*(m-5))+c(3)*alfa(j+1,2*(m-5)+1)+c(4)*alfa(j+1,2*(m-5)+2)+c(5)*alfa(j+1,2*(m-5)+3)+c(6)*alfa(j+1,2*(m-5)+4);
    beta(j,m-5)=h(1)*alfa(j+1,2*(m-5)-1)+h(2)*alfa(j+1,2*(m-5))+h(3)*alfa(j+1,2*(m-5)+1)+h(4)*alfa(j+1,2*(m-5)+2)+h(5)*alfa(j+1,2*(m-5)+3)+h(6)*alfa(j+1,2*(m-5)+4);
   
    alfa(j,m-4)=c(1)*alfa(j+1,2*(m-4)-1)+c(2)*alfa(j+1,2*(m-4))+c(3)*alfa(j+1,2*(m-4)+1)+c(4)*alfa(j+1,2*(m-4)+2)+c(5)*alfa(j+1,2*(m-4)+3)+c(6)*alfa(j+1,2*(m-4)+4);
    beta(j,m-4)=h(1)*alfa(j+1,2*(m-4)-1)+h(2)*alfa(j+1,2*(m-4))+h(3)*alfa(j+1,2*(m-4)+1)+h(4)*alfa(j+1,2*(m-4)+2)+h(5)*alfa(j+1,2*(m-4)+3)+h(6)*alfa(j+1,2*(m-4)+4);
   
    alfa(j,m-3)=c(1)*alfa(j+1,2*(m-3)-1)+c(2)*alfa(j+1,2*(m-3))+c(3)*alfa(j+1,2*(m-3)+1)+c(4)*alfa(j+1,2*(m-3)+2)+c(5)*alfa(j+1,2*(m-3)+3)+c(6)*alfa(j+1,2*(m-3)+4);
    beta(j,m-3)=h(1)*alfa(j+1,2*(m-3)-1)+h(2)*alfa(j+1,2*(m-3))+h(3)*alfa(j+1,2*(m-3)+1)+h(4)*alfa(j+1,2*(m-3)+2)+h(5)*alfa(j+1,2*(m-3)+3)+h(6)*alfa(j+1,2*(m-3)+4);
   
    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)+c(5)*alfa(j+1,2*(m-2)+3)+c(6)*alfa(j+1,2*(m-2)+4);
    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)+h(5)*alfa(j+1,2*(m-2)+3)+h(6)*alfa(j+1,2*(m-2)+4);
  
    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)+c(5)*alfa(j+1,2*m)+c(6)*alfa(j+1,2*m);
    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)+h(5)*alfa(j+1,2*m)+h(6)*alfa(j+1,2*m);
    
    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)+h(5)*alfa(j+1,2*m)+h(6)*alfa(j+1,2*m); % esto lo extiendo 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)+c(5)*alfa(j+1,2*m)+c(6)*alfa(j+1,2*m); % esto lo extiendo 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)+h(5)*alfa(j+1,3)+h(6)*alfa(j+1,4); % esto lo extiendo 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)+c(5)*alfa(j+1,3)+c(6)*alfa(j+1,4); % esto lo extiendo para que la frontera izquierda salga mas bonita
    betaaa(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,1)+h(5)*alfa(j+1,1)+h(6)*alfa(j+1,2); % esto lo extiendo para que la frontera izquerda salga mas bonita
    alfaaa(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,1)+c(5)*alfa(j+1,1)+c(6)*alfa(j+1,2); % esto lo extiendo para que la frontera izquierda salga mas bonita
    %Fin_fronteras

end

 [s1,s2]=size(beta);
 beta=zeros(s1,s2);

 p=m/2; %esta p no tiene nada que ver con la p de DB-2p
 alfaf=zeros(pasos+1,num);
 alfaf(1,:)=alfa(1,:);

for j=1:pasos
   p=2*p;    
      %Fronteras 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)+ r(6)*alfaaa(j,1) +d(6)*betaaa(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)+ r(5)*alfaaa(j,1) +d(5)*betaaa(j,1);
      alfaf(j+1,3)=r(2)*alfaf(j,2)+d(2)*beta(j,2)+ r(4)*alfaf(j,1)+d(4)*beta(j,1)+ r(6)*alfaa(j,1)+d(6)*betaa(j,1);
      alfaf(j+1,4)=r(1)*alfaf(j,2)+d(1)*beta(j,2)+ r(3)*alfaf(j,1)+d(3)*beta(j,1)+ r(5)*alfaa(j,1)+d(5)*betaa(j,1);
      %Fin_fronteras izquierda
      for i=3:p-1
        if(s(j,i)==0 && s(j,i-1)==0 && s(j,i-2)==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)+ r(6)*alfaf(j,i-2)+d(6)*beta(j,i-2);
            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)+ r(5)*alfaf(j,i-2)+d(5)*beta(j,i-2);
        end
        if (s(j,i)==0 && s(j,i-1)==0 && s(j,i-2)==1)
            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)+ r(6)*alfaf(j,i-2);
            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)+ r(5)*alfaf(j,i-2);
        end   
        if (s(j,i)==0 && s(j,i-1)==1 && s(j,i-2)==1)
            alfaf(j+1,2*i-1)=r(2)*alfaf(j,i)+d(2)*beta(j,i)+ r(4)*alfaf(j,i-1)+ r(6)*alfaf(j,i-2);
            alfaf(j+1,2*i)=r(1)*alfaf(j,i)+d(1)*beta(j,i)+ r(3)*alfaf(j,i-1)+ r(5)*alfaf(j,i-2);
        end  
        if (s(j,i)==1 && s(j,i-1)==1 && s(j,i-2)==1)
            %primero recuperamos lo de la derecha. los beta deben ser
            %anulados
            alfaf(j+1,2*i)=r(1)*alfaf(j,i)+ r(3)*alfaf(j,i-1)+ r(5)*alfaf(j,i-2);
            %ahora lo de la izquierda. los alfa deben ser extrapolados 
            alfafeji_2=extrapola([1,2,3],[alfaf(j,i-5),alfaf(j,i-4),alfaf(j,i-3)],4);
            alfafeji_1=extrapola([1,2,3],[alfaf(j,i-5),alfaf(j,i-4),alfaf(j,i-3)],5);
            alfafeji=extrapola([1,2,3],[alfaf(j,i-5),alfaf(j,i-4),alfaf(j,i-3)],6);
            alfaf(j+1,2*i-1)=r(2)*alfafeji+d(2)*beta(j,i)+ r(4)*alfafeji_1+d(4)*beta(j,i-1)+ r(6)*alfafeji_2+d(6)*beta(j,i-2);
        end
        if (s(j,i)==1 && s(j,i-1)==1 && s(j,i-2)==0 ) 
            alfafeji_1=extrapola([1,2,3],[alfaf(j,i-4),alfaf(j,i-3),alfaf(j,i-2)],4);
            alfafeji=extrapola([1,2,3],[alfaf(j,i-4),alfaf(j,i-3),alfaf(j,i-2)],5);
            alfaf(j+1,2*i-1)=r(2)*alfafeji+d(2)*beta(j,i)+ r(4)*alfafeji_1+d(4)*beta(j,i-1)+ r(6)*alfaf(j,i-2)+d(6)*beta(j,i-2);
            alfaf(j+1,2*i)=r(1)*alfafeji+d(1)*beta(j,i)+ r(3)*alfafeji_1+d(3)*beta(j,i-1)+ r(5)*alfaf(j,i-2)+d(5)*beta(j,i-2);                           
        end
        if(s(j,i)==1 && s(j,i-1)==0 && s(j,i-2)==0)      
            alfafeji=extrapola([1,2,3],[alfaf(j,i-3),alfaf(j,i-2),alfaf(j,i-1)],4);
            alfaf(j+1,2*i-1)=r(2)*alfafeji+d(2)*beta(j,i)+ r(4)*alfaf(j,i-1)+d(4)*beta(j,i-1)+ r(6)*alfaf(j,i-2)+d(6)*beta(j,i-2);
            alfaf(j+1,2*i)=r(1)*alfafeji+d(1)*beta(j,i)+ r(3)*alfaf(j,i-1)+d(3)*beta(j,i-1)+ r(5)*alfaf(j,i-2)+d(5)*beta(j,i-2);
        end         
      end 
      
      %Frontera derecha
      alfaf(j+1,2*p-1)=r(2)*alfaf(j,p)+d(2)*beta(j,p)+r(4)*alfaf(j,p-1)+d(4)*beta(j,p-1)+r(6)*alfaf(j,p-2)+d(6)*beta(j,p-2);
      alfaf(j+1,2*p)=r(1)*alfaf(j,p)+d(1)*beta(j,p)+r(3)*alfaf(j,p-1)+d(3)*beta(j,p-1)+r(5)*alfaf(j,p-2)+d(5)*beta(j,p-2);
      %Fin_frontera derecha
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 discretizada con 257 pts 
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
    
