clear all;
close all;
clc;


addpath('C:\Users\TT\Desktop\academics\FinalCode_UBSS\signals');
addpath('C:\Users\TT\Desktop\academics\FinalCode_UBSS\sparco-1.2');
addpath('C:\Users\TT\Desktop\academics\FinalCode_UBSS\l1_ls_matlab');
addpath('C:\Users\TT\Desktop\academics\FinalCode_UBSS\Solvers');
addpath('C:\Users\TT\Desktop\academics\FinalCode_UBSS\bss_eval');
addpath('C:\Users\TT\Desktop\academics\FinalCode_UBSS\spgl1-1.7');
addpath('C:\Users\TT\Desktop\academics\FinalCode_UBSS\omp');


t0 = clock;

%female
%mixing = [0.3754,  0.5235,  0.7660,  0.8151;...
%         1.0600   0.6129   0.5264   0.2967];

%DCol2 = sqrt(sum(mixing.*mixing,1));
%mixing = mixing .* repmat(1./DCol2,2,1);   
%male
%mixing = [ 0.3865    0.5391    0.7887    0.8394;...
%           1.0915    0.6311    0.5421    0.3055];

mm=wavread('male_inst_mix.wav');
mixing = kmeans(mm);

%%%%DCol2 = sqrt(sum(mixing.*mixing,1));
%%%%mixing = mixing .* repmat(1./DCol2,2,1);

data.signalSize=2048;   
      
block = data.signalSize;

data.M = opMatrix(kron(mixing,speye(data.signalSize,data.signalSize)), ...
              'Mixing matrix: kron(mixing,speye)');

 DCT = opFFT(512);
 D = DCT;

% m1=mm(:,1);m2=mm(:,2);
% mm1=buffer(m1,512,450,'nodelay');
% mm2=buffer(m2,512,450,'nodelay'); 
% TSM=[mm1,mm2];               
% DTr = SimCOTr(TSM);
% D = opMatrix(DTr);

WindowedD = opWindowedOp(data.signalSize, D, ones(512,1), 450);
data.B = opBlockDiag([1,1,1,1], WindowedD);  

seg=0;
ns = fix(160000/block);

for N=1:ns

mixture_block(:,1) = mm(seg+1:seg+block,1); mixture_block(:,2) = mm(seg+1:seg+block,2);

data.b = mixture_block(:);

data = completeOps(data);

opts = spgSetParms('optTol',1e-2);

tau = pi;    % Initial one-norm of solution

sigma = 0;  % Go for a basis pursuit solution

y = spgl1(data.A, data.b, tau, sigma, [], opts);
%[y, i2, ah2] = SolveMP(data.A, data.b, data.sizeA(2), 1000, 0, 0, 1e-2);
%At= opTranspose(data.A);
%lambda = 0.01;      % regularization parameter
%rel_tol = 1;     % relative target duality gap  opt 0.05
%m=data.sizeA(1);
%n=data.sizeA(2);
%w=2*m;
%[y,status,history]=l1_ls(data.A,At,w,n,data.b,lambda,rel_tol);
f = data.B(y,1);

es1(seg+1:seg+block)=f(1:block);
es2(seg+1:seg+block)=f(block+1:2*block);
es3(seg+1:seg+block)=f(2*block+1:3*block);
es4(seg+1:seg+block)=f(3*block+1:4*block);

seg=seg+block;

end

save test2007malefft es1 es2 es3 es4;

time = etime(clock, t0);

s1 = wavread('male_inst_sim_1.wav');
s2 = wavread('male_inst_sim_2.wav');
s3 = wavread('male_inst_sim_3.wav');
s4 = wavread('male_inst_sim_4.wav');

s1=s1';
s2=s2';
s3=s3';
s4=s4';

a1=s1(1,1:block*ns);
b1=s1(2,1:block*ns);
a2=s2(1,1:block*ns);
b2=s2(2,1:block*ns);
a3=s3(1,1:block*ns);
b3=s3(2,1:block*ns);
a4=s4(1,1:block*ns);
b4=s4(2,1:block*ns);

Y(1,:,1)=a1;
Y(2,:,1)=b1;
Y(1,:,2)=a2;
Y(2,:,2)=b2;
Y(1,:,3)=a3;
Y(2,:,3)=b3;
Y(1,:,4)=a4;
Y(2,:,4)=b4;

se1(1,:)=mixing(1,1)*es1;
se1(2,:)=mixing(2,1)*es1;

se2(1,:)=mixing(1,2)*es2;
se2(2,:)=mixing(2,2)*es2;

se3(1,:)=mixing(1,3)*es3;
se3(2,:)=mixing(2,3)*es3;

se4(1,:)=mixing(1,4)*es4;
se4(2,:)=mixing(2,4)*es4;

se(:,:,1) = se1;
se(:,:,2) = se2;
se(:,:,3) = se3;
se(:,:,4) = se4;

for i=1:4
    a=sum((se1(1,:)-Y(1,:,i)).^2)+sum((se1(2,:)-Y(2,:,i)).^2);b=sum((se2(1,:)-Y(1,:,i)).^2)+sum((se2(2,:)-Y(2,:,i)).^2);c=sum((se3(1,:)-Y(1,:,i)).^2)+sum((se3(2,:)-Y(2,:,i)).^2);d=sum((se4(1,:)-Y(1,:,i)).^2)+sum((se4(2,:)-Y(2,:,i)).^2);
    C=[a,b,c,d];
    [max,maxindex]=min(C,[],2);
    ees(:,:,i)=se(:,:,maxindex);
end

se1=ees(:,:,1);
se2=ees(:,:,2);
se3=ees(:,:,3);
se4=ees(:,:,4);

[s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se1,Y,1,512);
[SDR,ISR,SIR,SAR]=bss_mcrit(s_true,e_spat,e_interf,e_artif);
result1=[SDR,ISR,SIR,SAR];

[s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se2,Y,2,512);
[SDR,ISR,SIR,SAR]=bss_mcrit(s_true,e_spat,e_interf,e_artif);
result2=[SDR,ISR,SIR,SAR];

[s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se3,Y,3,512);
[SDR,ISR,SIR,SAR]=bss_mcrit(s_true,e_spat,e_interf,e_artif);
result3=[SDR,ISR,SIR,SAR];

[s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se4,Y,4,512);
[SDR,ISR,SIR,SAR]=bss_mcrit(s_true,e_spat,e_interf,e_artif);
result4=[SDR,ISR,SIR,SAR];


% wavwrite(es1,16000,'es1.wav');
% wavwrite(es2,16000,'es2.wav');
% wavwrite(es3,16000,'es3.wav');
% wavwrite(es4,16000,'es4.wav');

% [s_target,e_interf,e_artif] = bss_decomp_gain(es1,1,Y);
% [SDR,SIR,SAR] = bss_crit(s_target,e_interf,e_artif);
% result1=[SDR,SIR,SAR];
% 
% [s_target,e_interf,e_artif] = bss_decomp_gain(es2,2,Y);
% [SDR,SIR, SAR] = bss_crit(s_target,e_interf,e_artif);
% result2=[SDR,SIR,SAR];
% 
% [s_target,e_interf,e_artif] = bss_decomp_gain(es3,3,Y);
% [SDR,SIR,SAR] = bss_crit(s_target,e_interf,e_artif);
% result3=[SDR,SIR,SAR];
% 
% [s_target,e_interf,e_artif] = bss_decomp_gain(es4,4,Y);
% [SDR,SIR,SAR] = bss_crit(s_target,e_interf,e_artif);
% result4=[SDR,SIR,SAR];

result1
result2
result3
result4