%---------------------------------------------------------------
%                20070101 simulation para records 
%---------------------------------------------------------------
%function [tag_opt,tag_decode,Tag_diff] = run_sim(SNR_dB);

% -------------------- simulation parameter--------------------
m_transmitter = 4;  n_receiver = 8;         
fir_len = 1;        QAM_bit = 2;
N_point_FFT = 4;    symbol_block_K = 25;
SNR_dB = 30;

%------------------------trans symbols--------------------------
% S = transmitter(m_transmitter, QAM_bit, N_point_FFT, symbol_block_K);
load S;
%----------------------------channel----------------------------
chan_h = firChannel(m_transmitter, n_receiver, fir_len);
channel_H = gen_freq_H(chan_h, N_point_FFT);
R_HH = channel_H * channel_H';

%----------------------------step 1----------------------------
Y_opt = S * channel_H;      
Y = awgn(Y_opt,(SNR_dB + 10*log10(QAM_bit)),'measured');          %   received Signals Y

A = diag(gen_tag(m_transmitter, N_point_FFT));
Z_opt = S * A * channel_H;  
Z = awgn(Z_opt,(SNR_dB + 10*log10(QAM_bit)),'measured');          %   received Signals Z

R_YY = Y * Y';                      % received Signals Y s self-relations
R_YY_opt = Y_opt * Y_opt';          % Optimal received Signals Y's self-relations 
%R_YY0 = S * R_HH * S'; 

R_ZY = Z * Y';                      % received Signals Y & Z corelations
R_ZY_opt = Z_opt * Y_opt';          % Optimal received Signals corelations of Z & Y
%R_ZY0 = S*A*R_HH*S';

%----------------------------step 2-------------------------------
[VV, DD] = eig(R_YY); % R_YY * V = V * D
[VV_opt,DD_opt] = eig(R_YY_opt);
mu = diag(DD)';
mu_opt = diag(DD_opt)';

%----------------------------step 3------------------------------
sigma2 = sum(mu(m_transmitter * N_point_FFT + 1 : symbol_block_K))/...
    (symbol_block_K - m_transmitter * N_point_FFT);

%----------------------------step 4------------------------------
R_YY0_P = 0;
for index = 1 : (m_transmitter * N_point_FFT)
    temp = (mu(index) - sigma2)^(-1) * VV(:,index) * VV(:,index)';
    R_YY0_P = R_YY0_P + temp;
end;

R_YY0_P_opt = 0;
for index = 1 : (m_transmitter * N_point_FFT)
    temp = mu_opt(index)^(-1) * VV_opt(:,index) * VV_opt(:,index)';
    R_YY0_P_opt = R_YY0_P_opt + temp;
end;

R = R_ZY * R_YY0_P;
R_opt = R_ZY_opt * R_YY0_P_opt;     %optimal recieved signals
[V_R, D_R] = eig(R);
mu_R = diag(D_R);
[V_R_opt, D_R_opt] = eig(R_opt);
mu_R_opt = diag(D_R_opt);

% verify1:
result = R * V_R(:,1) - V_R(:,1) * D_R(1,1)
result_opt = R_opt * V_R_opt(:,1) - V_R_opt(:,1) * D_R_opt(1,1)

%----zhanGs new added 20070108
[D_R_opt_sort, V_R_opt_sort] = sort_D_V_opt(A, D_R_opt, V_R_opt);
test_opt = angle(A)/pi - angle(D_R_opt_sort(1:16,1:16))/pi;

[D_R_sort, V_R_sort] = sort_D_V(A, D_R, V_R);
test = angle(A)/pi - angle(D_R_sort(1:16,1:16))/pi;

% to do the "men xian pan jue"
%----decision
V_R_sort_decision = decision(V_R_opt_sort, V_R_sort, N_point_FFT);
test_div = V_R_sort_decision(:,:) ./ V_R_opt_sort(:,:)

% verify2:
result_opt_sort = R_opt * V_R_opt_sort - V_R_opt_sort * D_R_opt_sort;
result_sort = R * V_R_sort - V_R_sort * D_R_sort;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%