
%----------------Algoritmo para la separacion de fuentes con retardos-----------------
%
%   Implementacion de la version avanzada del algoritmo de separacion de fuentes con  
% pequeños retardos.
%   En esta version, se representan las fuentes originales y sus estimaciones. El 
% usuario puede comprobar el comportamiento del algoritmo y realizar simulaciones 
% variando los parametros que intervienen
%
%
%--------------------------------------------------------------------------------------
%
%
%   Variables utlizadas en la simulacion
%
%   Fuentes    Matriz cuyas filas contiene las fuentes originales
%   Rtotal     Matriz que contiene las matrices de autocovarianzas obtenidas con cada filtro 
%   Rx         Matriz de autocovarianzas de las observaciones
%   Rxh        Matriz de autocavarianzas filtrada de las observaciones
%   W          Matriz de separacion
%   x          Vector que contiene las observaciones recibidas en los receptores
%   y          Estimaciones de las fuentes originales a partir de las mezclas
%
%
%---------------------------------------------------------------------------------------

close all;
clear all;

%Lista de parametros de la simulacion

T=90000;        %Numero de muestras que vamos a tomar
Fs=48000;       %Frecuencia de muestreo
Rp=0.5;         %Atenuacion en la banda de paso del filtro
Rs=30;          %Atenuacion en la banda de rechazo del filtro
wp=2000/Fs;     %Frecuencia de corte del filtro LP que usamos para filtrar las señales
ws=2080/Fs;     %Frecuencia de salida del filtro LP que usamos para filtrar las señales
MR=14;           %Numero de recptores empleados (un nº adecuado es 3*N+2)
N=3;            %Numero de fuentes que vamos a mezclar para despues separarlas
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
A_ruido=0.001;  %Factor que controla la cantidad de ruido añadida

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

%Matriz que contiene los filtros necesarios para modelar los retrasos desde cada fuente
%a los receptores. Dichos filtros son utilizados de forma aleatoria por la funcion 
%genera_mezclas para introducir los retardos en las fuentes antes de mezclarlas
%Con los parametros de la simulacion que aparecen por defecto, pueden introducirse
%retardos de hasta 2 muestras, dentro de los limites teoricos del algoritmo
                              
B=[1 0 0 ;0 1 0 ;0 1 0; 1 0 0]; 

%Cargamos las señales de voz que se van a utilizar en la simulacion. Si se desea añadir 
%alguna fuente, solo hay que cargarla añadiendo otra linea como las que siguen y aumentar 
%en una unidad la constante que indica el numero de fuentes a separar

F=wavread('voz1',T); 
G=wavread('voz2',T);
O=wavread('voz4',T);

%Almacenamos todas las fuentes en la matriz Fuentes
%Si la energia de alguna de ellas difiere mucho del resto, la ajustamos multiplicando por
%un factor adecuado


Fuentes=[F';G';5*O'];

%Una vez obtenidas las fuentes, las filtramos para poder cumplir la condicion de pequeños
%retrasos

[M,wn]=ellipord(wp,ws,Rp,Rs);   %Generamos un filtro eliptico
[b,a]=ellip(M,Rp,Rs,wn);        %Obtenemos los coeficientes del filtro





%Ahora realizamos el filtrado con el filtro diseñado previamente

for i=1:N
    
    Fuentes_fil(i,:)=filter(b,a,Fuentes(i,:));
    
end

for i=1:N       %Representamos las fuentes de voz originales que queremos separar
    subplot(N,1,i);
    plot(Fuentes_fil(i,:));
    
    if i==1
        
        title('Señales originales');
        
    end
end

xlabel('Indice Temporal');

%La siguiente funcion genera las mezclas de señales con retrasos recibidas en los 
%receptores, que se almacenan en le vector de mezclas x


ruido=A_ruido*[rand(MR,T)-0.5]; %Ruido que vamos a introducir en las mezclas

%Añadimos ruido a las señales recibidas en cada sensor.Si se desea variar el nivel de ruido
%solo tenemos que variar la constante A_ruido

x=genera_mezcla(Fuentes_fil,B,MR,N);

x=x+ruido;

%A continuacion, se van extrayendo las fuentes una a una. Por defecto se extraen las 2N
%fuentes. Si solo queremos las fuentes de interes, bastara con extraer las N primeras

   
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
%Representamos las fuentes que hemos separado y la estimacion de sus derivadas

figure;

for i=1:2*N
    
    subplot(2*N,1,i);
    plot(y(i,1:T));
    
    if i==1
        
        title('Estimacion de las fuentes originales y de sus derivadas');
        
    end
    
    
end

xlabel('Indice temporal');
