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

randn('state',0);
rand('state',0);

itN = 50;
R1 = zeros(itN,3);
R2 = zeros(itN,3);
R3 = zeros(itN,3);
R4 = zeros(itN,3);
T = zeros(1,itN);
%mi=zeros(2,4,itN);
save DCTtest_temp R1 R2 R3 R4 T;

for itn = 1:itN
  
data.signalSize = 2048;       
block = data.signalSize;
t0 = clock;   
    
% randomly produce mixing matrix and normalize it.
while 1
        
  mixing = rand(2,4);
  mm = sqrt(sum(mixing.*mixing,1));
  mixing = mixing .* repmat(1./mm,2,1);

  aa=sum(mixing);
  aaa=(aa>ones(1,4));
  cc=cond(mixing);
  bb=sum((mixing(1,:)-mixing(2,:))>0);  
  tttt= (sum(abs(mixing(:,1)-mixing(:,2))>0.1))+(sum(abs(mixing(:,2)-mixing(:,3))>0.1))+(sum(abs(mixing(:,3)-mixing(:,4))>0.1))+(sum(abs(mixing(:,4)-mixing(:,1))>0.1))+(sum(abs(mixing(:,1)-mixing(:,3))>0.1))+(sum(abs(mixing(:,2)-mixing(:,4))>0.1));
  ttt=(abs(mixing(1,:)-mixing(2,:))>0.1);
  tt=(mixing>0.1);
  if sum(sum(tt))==8 && cc<4 && sum(ttt)==4 && tttt==12 && sum(aaa)==4 && bb==2
      break;
  end
  
end
       
s = zeros(50000,12); SN = 4;
s(:,1) = wavread('ukma.wav');
s(:,2) = wavread('frma.wav');
s(:,3) = wavread('itfe.wav');
s(:,4) = wavread('cnfe.wav');
s(:,5) = wavread('rufe.wav');
s(:,6) = wavread('gema.wav');
s(:,7) = wavread('nlma.wav');
s(:,8) = wavread('jpfe.wav');
s(:,9) = wavread('brfe.wav');
s(:,10) = wavread('esma.wav');
s(:,11) = wavread('dkma.wav');
s(:,12) = wavread('ukfe.wav');

% randomly pick SN source signal from the pool of size(s,2) source signals.

NS = 1:size(s,2);
for i = 1:SN
    rnd = ceil(rand(1)*(size(s,2)-i+1));
    Ns(i) = NS(rnd);
    NS(rnd) = [];% actually we pick by removing size(s,2)-N sources from the pool
end
   
S = s(:,Ns);

s1 = S(:,1);
s2 = S(:,2);
s3 = S(:,3);
s4 = S(:,4);

mixture = [s1 s2 s3 s4] * mixing';

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

DCT=opFFT(512);

D = DCT;

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


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

for N=1:ns

mixture_block(:,1) = mixture(seg+1:seg+block,1); mixture_block(:,2) = mixture(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);

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
 
time = etime(clock, t0);

s1=s1(1:block*ns);
s2=s2(1:block*ns);
s3=s3(1:block*ns);
s4=s4(1:block*ns);

es = [es1',es2',es3',es4'];
Y = [s1,s2,s3,s4]; 
Y= Y';

for i=1:4
    a=sum((es1-Y(i,:)).^2);b=sum((es2-Y(i,:)).^2);c=sum((es3-Y(i,:)).^2);d=sum((es4-Y(i,:)).^2);
    C=[a,b,c,d];
    [max,maxindex]=min(C,[],2);
    ees(:,i)=es(:,maxindex);
end

es1=ees(:,1)';
es2=ees(:,2)';
es3=ees(:,3)';
es4=ees(:,4)';


[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];

load DCTtest_temp R1 R2 R3 R4 T;

R1(itn,:) = result1;
R2(itn,:) = result2;
R3(itn,:) = result3;
R4(itn,:) = result4;

T(itn)=time;
%mi(:,:,itn)=mixing;

save DCTtest_temp R1 R2 R3 R4 T;
end

save DCTtest R1 R2 R3 R4 T;

DCTtest = load('DCTtest.mat');

avg_result1 = sum(MDCTtest.R1,1)/itN
avg_result2 = sum(MDCTtest.R2,1)/itN
avg_result3 = sum(MDCTtest.R3,1)/itN
avg_result4 = sum(MDCTtest.R4,1)/itN
avg_time = sum(DCTtest.T)/itN



