% Demonstration program of the C-FICA algorithm used with artificial
% colored sources and real filters

clear all; clc; close all;

% Number of samples
T=5e4;

% Choose the sources type
% 1: MA with uniforms innovation processes
% 2: MA with Laplacian innovation processes
% 3: AR with uniforms innovation processes
% 4: AR with Laplacian innovation processes
choice=2;

switch choice
    case 1
    s(1,:)=filter(rand(1,10),1,rand(1,T)-.5);
    s(2,:)=filter(rand(1,10),1,rand(1,T)-.5);
    case 2
    s(1,:)=filter(rand(1,10),1,laplace(1,T));
    s(2,:)=filter(rand(1,10),1,laplace(1,T));
    case 3
    s(1,:)=filter(1,[1 -.9],rand(1,T)-.5);
    s(2,:)=filter(1,[1 -.9],rand(1,T)-.5);
    case 4
    s(1,:)=filter(1,[1 -.9],laplace(1,T));
    s(2,:)=filter(1,[1 -.9],laplace(1,T));
end

% Definition of the filters by using the following data base:
% B. Gardner and K. Martin. Head Related Transfer Functions of a Dummy Head
% http://sound.media.mit.edu/ica-bench/.
A(:,1,:)=hrtf(30);
A(:,2,:)=hrtf(-80);

% Convolutive mixing process
x=conv_mix(A,s);

% C-FICA algorithm
[Sc,sc,Isc]=cfica(x,1e-16,80,400,'v');

% Computation of the true contributions
xc=contributions(A,s);

% Identification of the source permutations
perm=permutations(xc,sc);
sc=sc(:,perm,:);Isc=Isc(:,perm);

% SIRs computation
SIR1=sir(xc(Isc(1),1,:),sc(Isc(1),1,:))
SIR2=sir(xc(Isc(2),2,:),sc(Isc(2),2,:))

% Plotting of the true (blue) and estimated (red) contributions
figure
hold on
plot(squeeze(xc(Isc(1),1,:)));
plot(squeeze(sc(Isc(1),1,:)),'r');

figure
hold on
plot(squeeze(xc(Isc(2),2,:)));
plot(squeeze(sc(Isc(2),2,:)),'r');