%  ELIMINACIÓN DE RUIDO mediante AMC entre discontinuidades. Las discontinuidades se hallan mediante AP de stencil creciente.
%  Los grados de la AMC entre las discontinuidades detectadas se hallan
%  mediante el test de Durbin-Watson. 
%        capelimina1escamcentredwtest(x,f,fsin)
%        Entradas:
%            x    = nodos funcion de entrada
%            f    = f(x)
%            fsin = f pero sin ruido
%        Salidas:
%            plot de funcion original sin ruido (--), con ruido (.) y la reconstruccion (-)
%            se muestra por pantalla diferentes combinaciones para ECM y SNR
%            se muestran los grados de cada trozo
%
% Autor: José Jaime Noguera.  +INFO: www.uv.es/anims/noguera

function capelimina1escamcentredwtest(x,f,fsin)
[mx,nx]=size(x);
grau=3; %es el grado para hacer la deteccion de discontinuidades
hold on
inicio=1;
i=10;

%%%% DETECCION DE DISCONTINUIDADES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
salto(1)=1;
k=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(k)=i;
                k=k+1;
                inicio=i;
                i=i+9;
         end
         i=i+1;
end
salto(k)=nx+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[ns,ms]=size(salto);
DU=[1.684 1.693 1.704 1.715 1.725 1.735 1.746 1.757 1.768 1.779]; %esto es para 200 samples pero no he encontrado para más
                                                                  %en principiofunciona bien
for i=1:ms-1
        disp('salto')    
        coefi=polyfit(x(salto(i):(salto(i+1)-1)),f(salto(i):(salto(i+1)-1)),1);
        solu1(salto(i):(salto(i+1)-1))=polyval(coefi,x(salto(i):(salto(i+1)-1)));        
        em_1=solu1(salto(i):(salto(i+1)-1))-f(salto(i):(salto(i+1)-1));
        
        %CALCULO DE LOS GRADOS CON DWTEST
        m=2;
        entra=0;
        while (m<=9)   
            coefi=polyfit(x(salto(i):(salto(i+1)-1)),f(salto(i):(salto(i+1)-1)),m);
            solu1(salto(i):(salto(i+1)-1))=polyval(coefi,x(salto(i):(salto(i+1)-1)));
            em=solu1(salto(i):(salto(i+1)-1))-f(salto(i):(salto(i+1)-1));
            if estadistico_d(em) > DU(m)
                if abs(estadistico_d(em_1)-2)<abs(estadistico_d(em)-2)
                    m0=m-1;
                    m=11;
                    entra=1;
                else
                    m=m+1;
                    em_1=em;
                end    
            else
                m=m+1;
                em_1=em;
            end  
        end 
        if(entra==0)
           m0=10;
        end 
       disp('Grados hallados')
       m0
        %la aproximacion definitiva con el grado m0 es
        coefi=polyfit(x(salto(i):(salto(i+1)-1)),f(salto(i):(salto(i+1)-1)),m0);
        solu1(salto(i):(salto(i+1)-1))=polyval(coefi,x(salto(i):(salto(i+1)-1)));       
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,solu1);
    error
 
    disp('Entre fr y fn')
    error=mse(f,solu1);
    error
    
    disp('Entre fr y fn')
    error=calculosnr(solu1,f);
    error
    
    disp('Entre fin y fr')
    error=calculosnr(fsin,solu1);
    error
   
hold on
    
%axis([-1,1,-0.1,1.5])
plot(x, f,'k.','LineWidth',5,'MarkerSize',5)   
plot(x, fsin,'--k','LineWidth',2) 
plot(x(1:nx),solu1,'-k','LineWidth',3);

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

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


%SUBRUTINA PARA EL ESTADÍSTICO d
function solu=estadistico_d(em)
    n=length(em);
    numerador=0;
    denominador=0;
    for j=2:n
        numerador=numerador+(em(j)-em(j-1))^2;
        denominador=denominador+(em(j))^2;
    end
solu=numerador/denominador;    