clear all;
close all;

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

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


T_max = 1000;
T_min = 1;
SINR_SNR = zeros(1,(T_max/10)+1);
T_set = T_min:10:T_max;

for l=0:(nbrloops-1),
	for k=1:size(T_set,1),
		%loopbegin = cputime - t0
		T = T_set(k);

		[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_SNR(k) = SINR(1);
		%loopend = cputime - t0
	end;

	%plotbegin = cputime - t0
	SINR_SNR_f(l+1,:) = SINR_SNR;
end;

for k=1:size(T_set,1),
	SINR_SNR_m(k) = mean(SINR_SNR_f(:,k));
end;

SINR_SNR_m

%figure(1)
%plot((SNR_min:2:SNR_max),SINR_SNR_m,'r')
%hold on
%plot((SNR_min:2:SNR_max),(SNR_min+10*log10(N):2:SNR_max+10*log10(N)),'b')
%xlabel('SNR in dB');
%ylabel('SINR of signal 1 in dB');
%title(['Influence of SNR for T = 1000 and for nbrloops = ',num2str(nbrloops)]);