% Daubechies 4. Todos los coeficientes beta se anulan. Las fronteras se tratan como periódocas.
%     preliwavedb4u(x,f,pasos)
%     Entradas:
%         x   = nodos iniciales
%         f   = f(x)
%         pasos = numero de niveles que aplicaremos la transformada wavelet
%     Salidas:
%         Plot de la funcion original discretizada con 256 puntos (--), los (x,f(x)) iniciales (o), la reconstrucción (-) y la diferencia entre la original i la reconstrucción (:)

function  preliwavedb4u(x,f,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,:));

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*i+1)+h2*alfa(j+1,2*i+2)+h3*alfa(j+1,1)+h4*alfa(j+1,2);
    alfa(j,m)=c1*alfa(j+1,2*i+1)+c2*alfa(j+1,2*i+2)+c3*alfa(j+1,1)+c4*alfa(j+1,2);
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 periodicas
      alfaf(j+1,1)=c1*alfaf(j,1)  +h1*beta(j,1) + c3*alfaf(j,p) +h3*beta(j,p);       
      alfaf(j+1,2)=c2*alfaf(j,1)  +h2*beta(j,1) +  c4*alfaf(j,p)+h4*beta(j,p);
         
      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
      
      %Fronteres periodicas
      alfaf(j+1,2*p-1)=c1*alfaf(j,p)+h1*beta(j,p)+c3*alfaf(j,p-1)+h3*beta(j,p-1);
      alfaf(j+1,2*p)  =c2*alfaf(j,p)+h2*beta(j,p)+c4*alfaf(j,p-1)+h4*beta(j,p-1);
end
 
[xor,yor]=funtreballtesi(256,0); 

hold on
plot(xor, yor,'--k','LineWidth',2)                                   %original discretizada con 257 pts 
plot(x,f-alfaf(pasos+1,:),':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