% Eliminación de ruido mediante AP punto a punto con stencil centrado excepto en las discontinuidades (se detectan mediante AP de stencil creciente).
%   capelimina1escamcpapcentexcdisc(x,f,fsin,k,n)
%   Entradas
%      x    = nodos iniciales
%      f    = f(x)
%      fsin = f sin ruido
%      k    = grado para el polinomio de la AP
%      n    = tamaño del stencil (debe ser impar)
%   Salidas
%      Plot de la funcion con ruido, sin ruido y la reconstruccion
%      ECM y SNR segun varias combinaciones
%
% Autor: José Jaime Noguera.   +INFO: www.uv.es/anims/noguera

function capelimina1escamcpapcentexcdisc(x,f,fsin,k,n)
[mx,nx]=size(x);

if(mod(n,2)==0)
    disp('El tamaño del stencil debe ser impar')
end    

%%%%%%%%%%%%%%%%SUBPROGRAMA AP STENCIL CRECIENTE%%%%%%%%%%%%%%%%%%%
hold on
inicio=1;
i=10;
grau=3;
 salto(1)=1;
 l=2;
while i<nx-1
         coefi=qk(x(inicio:i),f(inicio:i),grau);
         [nc,mc]=size(coefi);
         coefi1=coefi(nc:-1:1);
         solu1(inicio:i)=polyval(coefi1,x(inicio:i));
         suma=0;
         for j=inicio:i
            suma=suma+(solu1(j)-f(j))^2;
         end
         error(i)=suma;
         if (i>inicio)
            errmed(i)=abs(error(i)-error(i-1));
         end  
         if (i>10 && i~=inicio+11 && i~=inicio+10 && errmed(i)>3*max(errmed(inicio+1:i-1)))
              salto(l)=i;
              %plot(x(i),0,'k*')
                   l=l+1;
                   inicio=i;
                   i=i+9;
         end   
         i=i+1;
end
salto(l)=nx+1; % El último salto es el último de los datos iniciales
  
%%FIN DETECCION DISCONTINUIDADS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%SUBPR0GRAMA PARA AP CENTRADO EXCEPTO EN DISCONTINUIDADES%%%%%%%%%
[ns,ms]=size(salto);

for j=2:ms
   for i=salto(j-1)+(n-1)/2:salto(j)-(n-1)/2-1
        coef=qk(x(i-(n-1)/2:i+(n-1)/2),f(i-(n-1)/2:i+(n-1)/2),k);
        [c1,c2]=size(coef);
        coef1=coef(c1:-1:1);
        fa(i)=polyval(coef1,x(i));  
   end
   kk=0;
   for i=salto(j-1):salto(j-1)+(n-1)/2-1
        coef=qk(x(i-kk:i+n-1-kk),f(i-kk:i+n-1-kk),k);
        [c1,c2]=size(coef);
        coef1=coef(c1:-1:1);
        fa(i)=polyval(coef1,x(i));
        kk=kk+1;
   end   
   kk=0;
   for i=salto(j)-(n-1)/2:salto(j)-1
        kk=kk+1;
        coef=qk(x(i-(n-1)/2-kk:i+(n-1)/2-kk),f(i-(n-1)/2-kk:i+(n-1)/2-kk),k);
        [c1,c2]=size(coef);
        coef1=coef(c1:-1:1);
        fa(i)=polyval(coef1,x(i));
   end   
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %fr la reconstruccion
    %fsin la sin ruido , el ideal a reconstruir
    %fn la funcion original con ruido

    disp('Entre fr y fsin')
    error=mse(fsin,fa);
    error
 
    disp('Entre fn y fr')
    error=mse(f,fa);
    error
    
    disp('Entre fr y fn')
    error=calculosnr(fa,f);
    error
    
    disp('Entre fin y fr')
    error=calculosnr(fsin,fa);
    error
   
hold on   
plot(x, f,'k.','LineWidth',5,'MarkerSize',5)   
plot(x, fsin,'--k','LineWidth',2) 
plot(x,fa,'-k','LineWidth',3);

%para trabajar con funcio
%axis([-1 1 -1.2 1.2 ]);
%set(gca,'Xtick',[ -1 -0.5 0 0.5 1], 'FontSize',15)
%set(gca,'Ytick',[ -1 -0.5 0 0.5 1], 'FontSize',15)

%para trabajar con funtreballtesi
%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