function y = Retrasos(x,N,TIPO)


%------------------Funcion para la separacion ciega de fuentes con retardos---------------------
%
%
%       Esta funcion implementa las dos versiones del algoritmo de separacion de fuentes con
%   retardos. Recibe como entrada el vector de mezclas de fuentes con retardos, el numero de
%   fuentes que existen en las mezclas y la constante TIPO que indica la version del algoritmo
%   de separacion que se desea usar para la separacion de las fuentes.
%
%------------------------------------------------------------------------------------------------
%
%   ARGUMENTOS DE SALIDA
%
%   y       Vector de fuentes estimadas    y=W*x
%
%   ARGUMENTOS DE ENTRADA
%
%   x       Vector de mezclas
%   N       Numero de fuentes presentes en las mezclas
%   TIPO    Parametro que indica la version del algoritmo (TIPO=1 V. Basica, TIPO=2 V. Avanzada)
%
%-----------------------------------------------------------------------------------------------

%Constantes

Nit=1400;       %Numero de iteraciones que vamos a usar para el calculo de los autovectores
NF=2;           %Variable que indica el numero de filtros que vamos a utilizar
Umbral=.5e-3;   %Constante para determinar si la matriz de mezclas se genero bien

%Matriz que contiene los filtros que vamos a utilizar para obtener cada una de las 
%matrices de autocovarianzas filtradas
%Si se desea añadir un filtro, solo hay que añadir una fila con el filtro deseado e 
%indicar el nuevo valor de la variable NF

Filtros=[0,ones(1,15);0 1 1 ones(1,10) 0  0 0;0 1 0 1 0 ones(1,11)];

%Matriz en la que almacenaremos las submatrices Ri obtenidas con cada filtro
                       
Rtotal=zeros(2*N,2*N*NF);

%Determinamos el numero de receptores y la longitud de las muestras

[MR,T]=size(x);

%Si TIPO=1, aplicamos la version basica del algoritmo de separacion

if TIPO==1
    
    Rx=(1/T)*x*x';   %Determinamos Rx(0)
    [U,L,V]=svd(Rx); %Descomponemos en valores singulares la matriz Rx


    %Seguidamente obtenemos la matriz de autocovarianzas filtrada. Para ello filtramos
    %cada una de las señales recibidas en los sensores por separado

    for i=1:MR
    
        a=filter(Filtros(1,:),1,x(i,1:T-1));
        if i==1
            k=a;
        else
            k=[k;a];
        end    
    end


    Rxh=(1/T)*k*k';
    U=U(1:MR,1:2*N); %Generamos la matriz U deseada
    L=L-trace(L(2*N+1:MR,2*N+1:MR))*(1/(MR-2*N)); %Estimamos la varianza del ruido y la restamos
    A=L(1:2*N,1:2*N);   %Computamos la matriz A 

    %A partir de ella, determinamos la matriz auxiliar R

    R=inv(A)*U'*Rxh*U;

    %Obtenemos los autovalores de R

    [E,D]=eig(R);      

                   
    %Matriz de separacion. Tenemos que multiplicar las señales por esta matriz para obtener las 
    %originales

    W=E'*U';         


y=W*x;

%Por defecto, usaremos la version avanzada del algoritmo

else
    
    y=zeros(2*N,T);
    
    for t=1:2*N

   
    Rx=(1/T)*x*x';      %Determinamos Rx(0)
    [U,L,V]=svd(Rx);    %Descomponemos en valores singulares la matriz Rx

     
    U=U(1:MR,1:2*N);    %Generamos la matriz U deseada
    
    %Estimamos la varianza del ruido y la restamos
    
    L=L-trace(L(2*N+1:MR,2*N+1:MR))*(1/(MR-(2*N))); 
    A=L(1:2*N,1:2*N);   %Computamos la matriz A 

    %A partir de ella, determinamos la matriz R

    
    %Calculamos cada una de las matrices de autocovarianzas filtrada. Esto lo hacemos
    %filtrando todas las observaciones con cada filtro.
    %Las matrices de autocovarianzas obtenidas se almacenan en la matriz Rtotal

    for j=0:NF-1
        
        for i=1:MR
    
            a=filter(Filtros(j+1,:),1,x(i,1:T));
            if i==1
                k=a;
            else
                k=[k;a];
            end    
        end
        
        Rxh=(1/T)*k*k';
        D=eig(A);
        
        if sum(D)<Umbral
            
            A=eye(2*N);
            
        end
        
        R=inv(A)*U'*Rxh*U;
        
        %Aqui vamos almacenando cada una de las matrices Ri auxiliares
        
        Rtotal(:,j*2*N+1:(j+1)*2*N)=R; 
    end
    
    [E1,D1]=eig(R); %Obtenemos los autovalores de R1
        
   %Asi damos un valor inicial a v1 proximo al real que queremos estimar

    v1=E1(:,t);           
    r=zeros(2*N,1);
    
    %A partir de las diversas matrices de autocorrelacion, obtenemos una estimacion mas
    %precisa de v1
    
            for j=1:Nit
                
                for i=0:NF-1
                    
                    r=r+conj(v1'*Rtotal(:,2*N*i+1:2*N*(i+1))*v1)*Rtotal(:,2*N*i+1:2*N*(i+1))*v1;
                    
                end
    
                
                v1=r/norm(r);
                r=zeros(2*N,1);
                
            end
     
      
        W=v1'*U';   %Matriz de separacion. Tenemos que multiplicar por esta matriz para
                    %obtener las fuentes originales
                    
    
        y(t,:)=W*x; %Obtenemos las fuentes originales y las almacenamos en la matriz C.
        
        %Seguidamente realizamos la operacion de deflacion para elimnar de la mezcla la fuente
        %que previamente hemos extraido
        
        y1=W*x;
        d=x*y1'/T;
        Pa_=eye(MR)-d*pinv(d); %matriz de proyeccion ortogonal
        x=Pa_*x;                         

    end
end