% Multirresolución mediante la transformada ENO-wavelet DB-2p utilizando la TRANSFORMACION INVERSA. Igualamos totos los coeficentes beta a 0.
%         capenowavedb2ptifinal(x,f,pasos,p,a,bmayor)
%         Entradas:
%            x      =  nodos iniciales
%            f      =  f(x)
%            pasos  =  niveles que bajamos
%            p      =  el p de DB-2p
%            a      =  condicion para localizar las discontinuidades (|B_{j,i}|>=a|B_{j,i-1}|) 
%            bmayor =  condicion para localizar las discontinuidades (B_{j,i}>bmayor) 
%         Salidas:
%            plot de la reconstruccion (-), de la funcion inicial (--) y de
%            la diferencia entre ambas (:), desplazado 0.1 en el eje Y
%
% Autor: José Jaime Noguera.   +INFO: www.uv.es/anims/noguera

function  [normainfinito,norma2]=capenowavedb2ptifinal(x,f,pasos,p,a,bmayor)
%Este algoritmo no está optimizado en cuanto a coste computacional y de
%almacenamiento

%OJO con estos valores. En ocasiones aparecen malas aproximaciones debidas
%a que esta elección no es idónea para detectar las discontinuidades de la
%función propuesta.
if nargin ==4
   a=2;
   bmayor=0.0001;
end   
num=length(f);
alfa(pasos+1,:)=f(:)';
beta(1:pasos,1:length(x))=zeros(pasos,length(x));

[c,h]=wfilters(['db' num2str(p)],'d');

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;   
   % Datos que se utilizarán después, con independencia que algunos se
   % sobreescribirán

   for i=1:m-3*(p-1)
           suma1=0;
           suma2=0;
           for k=2*i-1:2*i+2*p-2
              suma1=suma1+h(k-2*i+2)*alfa(j+1,k);
              suma2=suma2+c(k-2*i+2)*alfa(j+1,k);
           end 
           beta(j,i)=suma1;
           betadecide(j,i)=suma1; %este se utilizara para localizar la discontinuidad. Así aseguramos que no se sobreescriben datos.
           alfa(j,i)=suma2;
           suma=0;
   end    
   
   %Reiniciemos el índice i
   i=1;  
   
   while(i<=m-3*(p-1))         
         if(i>p  && abs(betadecide(j,i))>=bmayor && abs(betadecide(j,i-1))<=a*abs(betadecide(j,i-2))&& abs(betadecide(j,i))>a*abs(betadecide(j,i-1)) && abs(betadecide(j,i+p-1))<=a*abs(betadecide(j,i+p)))
             %hay p-1 stencils afectados
             %lado izquierdo
             for k=i:(i+p-2) 
                alfa(j,k)=extrapola((1:p),alfa(j,(i-p):(i-1)),p+k-i+1); 
             end 
             for l=0:p-2 
                sum=0; 
                for k=2:p 
                    sum=sum+r(2*k)*alfa(j,i-k+l+1)+d(2*k)*beta(j,i-k+l+1); 
                end 
                beta(j,i+l)=(1/d(2))*(alfa(j+1,2*i-1+2*l)-r(2)*alfa(j,i+l)-sum); 
             end 
             %lado derecho
             betaaux=beta; %para evitar que se sobrescriban datos necesarios
             for l=(p-2):-1:0 
                  betaaux(j,i+l)=0; 
             end
             for l=(p-2):-1:0 
                 sum=0; 
                 for k=0:p-2 
                      sum=sum+ r(2*k+1)*alfa(j,i-k+p+l-1)+d(2*k+1)*betaaux(j,i-k+p+l-1); 
                 end 
                 alfa(j,i+l)=(1/r(2*p-1))*(alfa(j+1,2*(i+l)-1+2*p-1)-sum); %este se puede sobreescribir
             end 
             for k=i:i+p-2 
                 s(j,k)=1; 
             end 
             i=i+p-2;                   
         else
             if( i>2 && abs(betadecide(j,i))>=bmayor && abs(betadecide(j,i-1))<=a*abs(betadecide(j,i-2)) && abs(betadecide(j,i))>a*abs(betadecide(j,i-1)) && abs(betadecide(j,i+p-1))>a*abs(betadecide(j,i+p)))
                 %p stencils afectados
                 %lado izquierdo  
                 for k=i:(i+p-1) 
                     alfa(j,k)=extrapola([1:p],[alfa(j,(i-p):(i-1))],p+k-i+1); 
                 end
                 for l=0:p-1
                    sum=0; 
                    for k=2:p 
                        sum=sum+r(2*k)*alfa(j,i-k+l+1)+d(2*k)*beta(j,i-k+l+1); 
                    end 
                    beta(j,i+l)=(1/d(2))*(alfa(j+1,2*i-1+2*l)-r(2)*alfa(j,i+l)-sum);
                 end
                 %lado derecho
                 betaaux=beta; %para evitar que se sobreescriba las beta anteriores
                 for l=(p-1):-1:0  
                     betaaux(j,i+l)=0; 
                 end 
                 for l=(p-1):-1:0 
                     sum=0; 
                     for k=0:p-2 
                        sum=sum+ r(2*k+1)*(alfa(j,i-k+p+l-1))+d(2*k+1)*betaaux(j,i-k+p+l-1); 
                     end 
                     alfa(j,i+l)=(1/r(2*p-1))*(alfa(j+1,2*(i+l)-1+2*p-1)-sum); %este si se puede sobreescribir
                 end  
                 for k=i:i+p-1 
                      s(j,k)=1; 
                 end 
                 i=i+p-1;
             
             else
                 %si no pues no hay discontinuidad y calcula lo normal (la beta ya esta de antes hecho)                
                 suma=0;
                 for k=2*i-1:2*i+2*p-2
                       suma=suma+c(k-2*i+2)*alfa(j+1,k);
                 end 
                 alfa(j,i)=suma; 
             end    
         end
       i=i+1;       
   end
    
       %TODO ESTO ES PARA ARREGLAR FRONTERAS
       for l=m-3*(p-1)+1:m-p+1
               suma1=0;
               suma2=0;
               for k=1:2*p
                  suma1=suma1+c(k)*alfa(j+1,2*l-1+k-1);
                  suma2=suma2+h(k)*alfa(j+1,2*l-1+k-1);
               end 
               alfa(j,l)=suma1;
               beta(j,l)=suma2;
        end 
          
        %frontera derecha
        for l=m-p+2:m
               suma1=0;
               suma2=0;          
               for k=1:2*p
                  if(2*l-1+k-1>2*m)
                      suma1=suma1+c(k)*alfa(j+1,2*m);
                      suma2=suma2+h(k)*alfa(j+1,2*m);
                  else    
                      suma1=suma1+c(k)*alfa(j+1,2*l-1+k-1);
                      suma2=suma2+h(k)*alfa(j+1,2*l-1+k-1);
                  end
               end 
               alfa(j,l)=suma1;
               beta(j,l)=suma2;
         end

         %frontera izquierda
         index=1;
         for l=-(p-2):0
               suma1=0;
               suma2=0;          
               for k=1:2*p
                  if(2*l-1+k-1<=0)
                      suma1=suma1+c(k)*alfa(j+1,1);
                      suma2=suma2+h(k)*alfa(j+1,1);
                  else    
                      suma1=suma1+c(k)*alfa(j+1,2*l-1+k-1);
                      suma2=suma2+h(k)*alfa(j+1,2*l-1+k-1);
                  end
               end 
               alfaa(j,index)=suma1;
               betaa(j,index)=suma2;
               index=index+1;
         end 
          
end
 
[s1,s2]=size(beta);
beta=zeros(s1,s2); %APROXIMACION LINEAL

longpas=m/2; 
alfaf=zeros(pasos+1,num);
alfaf(1,:)=alfa(1,:);
 
betaaux=beta; 
S=zeros(2*p,p); S(2,1)=1; 
S=full(spdiags(ones(p,p),[-p:-1]',S)); 

for j=1:pasos
    longpas=2*longpas; 
    alfaaux(j,:)=alfaf(j,:);
   
    %Frontera izquierda
    for i=1:p-1
          sum1=0; 
          sum2=0;
          index=p-1;
          for k=0:(p-1) 
                  if(i-k<=0)
                    sum1=sum1+r(2*k+2)*alfaa(j,index)+d(2*k+2)*betaa(j,index); 
                    sum2=sum2+r(2*k+1)*alfaa(j,index)+d(2*k+1)*betaa(j,index);
                    index=index-1;
                  else    
                    sum1=sum1+r(2*k+2)*alfaf(j,i-k)+d(2*k+2)*beta(j,i-k); 
                    sum2=sum2+r(2*k+1)*alfaf(j,i-k)+d(2*k+1)*beta(j,i-k); 
                  end  
          end 
          alfaf(j+1,2*i-1)=sum1; 
          alfaf(j+1,2*i)=sum2; 
    end     
    %Frontera derecha
    for i=longpas-p+2:longpas
          sum1=0; 
          sum2=0;
          for k=0:(p-1) 
   
                    sum1=sum1+r(2*k+2)*alfaf(j,i-k)+d(2*k+2)*beta(j,i-k); 
                    sum2=sum2+r(2*k+1)*alfaf(j,i-k)+d(2*k+1)*beta(j,i-k); 
              
          end 
          alfaf(j+1,2*i-1)=sum1; 
          alfaf(j+1,2*i)=sum2; 
    end  
    
    %La parte importante:
    for i=p:longpas-p+1
      for l=1:2*p 
              if( s(j,i-p+1:i)==S(l,:)) 
                    if(l==1) 
                           sum1=0; 
                           sum2=0; 
                           for k=0:(p-1) 
                                 sum1=sum1+r(2*k+2)*alfaf(j,i-k)+d(2*k+2)*beta(j,i-k); 
                                 sum2=sum2+r(2*k+1)*alfaf(j,i-k)+d(2*k+1)*beta(j,i-k); 
                           end 
                           alfaf(j+1,2*i-1)=sum1; 
                           alfaf(j+1,2*i)=sum2; 
                    end 
               
                    if( l>1 && l<=p) 
                           for ii=(i-p+1):(i-p+l-1) 
                                  beta(j,ii)=0; 
                           end 
                           sum1=0; 
                           sum2=0; 
                           for k=0:(p-1) 
                                  sum1=sum1+r(2*k+2)*alfaf(j,i-k)+d(2*k+2)*beta(j,i-k); 
                                  sum2=sum2+r(2*k+1)*alfaf(j,i-k)+d(2*k+1)*beta(j,i-k);  
                           end 
                           alfaf(j+1,2*i-1)=sum1; 
                           alfaf(j+1,2*i)=sum2; 
                           alfaf(j,:)=alfaaux(j,:); 
                           beta(j,:)=betaaux(j,:); 
                    end 
    
                    if(l==p+1) 
                           for k=(i-p+1):i 
                                    alfaf(j,k)=extrapola([1:p],[alfaf(j,(i-2*p+1):(i-p))],p+k-i+p); 
                           end 
                           sum1=0;                   
                           for k=0:(p-1) 
                                      sum1=sum1+r(2*k+2)*alfaf(j,i-k)+d(2*k+2)*beta(j,i-k); 
                           end 
                           alfaf(j+1,2*i-1)=sum1; 
                           alfaf(j,:)=alfaaux(j,:); 
                           beta(j,:)=betaaux(j,:); 
                           for k=(i-p+1):i 
                                  beta(j,k)=0; 
                           end 
                           sum2=0; 
                           for k=0:(p-1) 
                                  sum2=sum2+r(2*k+1)*alfaf(j,i-k)+d(2*k+1)*beta(j,i-k); 
                           end 
                           alfaf(j+1,2*i)=sum2; 
                           alfaf(j,:)=alfaaux(j,:); 
                           beta(j,:)=betaaux(j,:); 
                    end 
   
                    if(l>p+1) 
                           for k=(i-2*p+l):i 
                                   alfaf(j,k)=extrapola([1:p],[alfaf(j,(k-p):(k-1))],p+1); 
                           end 
                           sum1=0; 
                           sum2=0; 
                           for k=0:(p-1)  
                                   sum1=sum1+r(2*k+2)*alfaf(j,i-k)+d(2*k+2)*beta(j,i-k); 
                                   sum2=sum2+r(2*k+1)*alfaf(j,i-k)+d(2*k+1)*beta(j,i-k); 
                           end 
                           alfaf(j+1,2*i-1)=sum1; 
                           alfaf(j+1,2*i)=sum2; 
                           alfaf(j,:)=alfaaux(j,:); 
                           beta(j,:)=betaaux(j,:); 
                    end 
              end   
      end
         
   end 

end

%disp('Norma infinito');
%normainfinito=norm(alfaf(pasos+1,:)-f(:)',inf);
lleva=p*2^(pasos+1); %esto son los valores de la frontera que no consideramos para el cálculo de órdenes, ya que en la frontera no se mantiene el orden.
normainfinito=norm(alfaf(pasos+1,lleva:num-lleva)-f(lleva:num-lleva),inf);
%disp('Norma 2');
%norma2=norm(alfaf(pasos+1,:)-f(:)',2);
norma2=norm(alfaf(pasos+1,lleva:num-lleva)-f(lleva:num-lleva),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
legend('Original','Reconstrución')
%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