
%----------------Algoritmo para la separacion de fuentes con retardos-----------------
%
%   Implementacion de la version basica 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
%   x           Vector que contiene las observaciones recibidas en los receptores
%   Rx          Matriz de autocovarianzas de las observaciones
%   Rxh         Matriz de autocavarianzas filtrada de las observaciones
%   W           Matriz de separacion
%   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=11;          %Numero de receptores empleados
N=3;            %Numero de fuentes que vamos a mezclar para despues separarlas
A_ruido=0.001;  %Factor que controla la cantidad de ruido añadida
filtro=[0, ones(1,10)]; %Filtro con respuesta en cero nula utilizado para filtrar Rx

%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 1; 1 0;0 1];  
                              

%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';10*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 de las fuentes 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

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


%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=x+A_ruido*[rand(MR,T)-0.5];  
    
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(filtro,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;

%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('Estimaciones de las fuentes transmitidas y de sus derivadas')
        
    end
    
    
end


xlabel('Indice temporal');


