clear all;
close all;

randn('state',0);
rand('state',0);
t0 = cputime;
nbrloops = 1;

loopbegin = cputime - t0

%L = 2; 			%% # of signals
N = 8; 			%% # of sensors
T = 1000;		%% # of samples of X (snapshots)
	     		%% The matrix C = U * U' with U = W * A looks better when T > 1000 
K = 10; 		%% # of sample covariances matrices of Z, # of delays
thold = 1.0e-12; 	%% Just a threshold for the joint approximate diagonalization 
x1 = 550; 		%% Start of the interval for the plot
x2 = 650;		%% End of the interval for the plot
SNR = 10;


L_max = 7;
L_min = 2;
SINR_SNR = zeros(1,L_max);
L_set = L_min:1:L_max;

for l=0:(nbrloops-1),
	for k=0:size(L_set,2)-1,
		%loopbegin = cputime - t0
		L = L_set(k+1);

		[X, S, A, Ps, Pw, Noise] = signal(T, N, L, SNR);

		[Rx, Lambda, H, Lambda_small] = cov_sample(X, T, N, L);

		[sigma, W] = whitening(Lambda_small, Lambda, H, N, L);

		[U, C] = check_w(W, A);

		[Z] = matrix_z(W, X);

		[R] = estimate_R(Z, T, K, L);

		[V, D] = joint_diag(R,thold);

		[SE, AE, QE, IW] = final_estimation(V, W, X);

		[xcorr, SEF, AEF, QEF] = max_corr(S, SE, AE, QE, N, L, T);

		[SINR, R_s, R_in] = sinr(S, A, QEF, Ps, Pw, L, N, T);

		SINR_L(k+1) = SINR(1);
		%loopend = cputime - t0
	end;

	plotbegin = cputime - t0
	SINR_L_f(l+1,:) = SINR_L
end;

for k=0:size(L_set,2)-1,
	SINR_L_m(k+1) = mean(SINR_L_f(:,k+1));
end;

SINR_L_m

figure(1)
plot((L_min:10:L_max),SINR_L_m,'r')
%hold on
%plot((T_min:10:T_max),SINR_T_m(99)*ones(1,size(T_min:10:T_max,2)),'b') 
xlabel('Number of signals L');
ylabel('SINR of signal 1 in dB');
title(['Influence of L for SNR = 10dB, T = 1000, K = 10 and for nbrloops = ',num2str(nbrloops)]);