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);
time_start = clock; 

m = 16; %the dimension of the training signal
d = 32; %the number of atoms in training dictionary
c = 1; %the number of nonzero elements in coefficient vector
T = 5120;

A = rand(2,4);
S = rand(4,T);
X = A*S;

%generate a dictionary and normalize it
D = randn(m,d);
DCol2 = sqrt(sum(D.*D,1));
D = D .* repmat(1./DCol2,m,1);

SR = pinv(A)*X;
mixing = A;

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

mixture = X;

data.signalSize = 5120;   
      
block = data.signalSize;

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

DCT = opDCT(16);

WindowedD = opWindowedOp(data.signalSize, DCT, ones(16,1), 8);
data.B = opBlockDiag([1,1,1,1], WindowedD);  

data.b = mixture(:);
 
data = completeOps(data);
 
%opts = spgSetParms('optTol',1e-4,'iterations',100);

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

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=f(1:block);
es2=f(block+1:2*block);
es3=f(2*block+1:3*block);
es4=f(3*block+1:4*block);
 
%time = etime(clock, t0);

es = [es1 es2 es3 es4];

Xr = X - A*es';
ftt = sum(sum(Xr.*Xr));

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

result1
result2
result3
result4

