% == Nov. 2009 ==
% By Chien-Wei Hsin, Tsung-Hui Chang, Wing-Kin Ma
%
% This program presents a demo of the specialized interior-point algorithm  (IPA) developed 
% in the following paper
%
% Tsung-Hui Chang, Chien-Wei Hsin, Wing-Kin Ma and Chong-Yung Chi "A linear fractional semidefinite relaxation 
% approach to maximum-likelihood detection of higher-order QAM OSTBC in
% unknown channels", to appear in IEEE Trans. Signal Processing,
%
% for blind ML detection of higher-order QAM OSTBC.
% (All the equation numbers and paper references mentioned in the demo program are based on the above paper.) 
% 
% In this demo we consider the following MIMO OSTBC system model
%      Yp=HC(sp)+Wp  (Nr X T), p=1,....,P.    (1)
% Here, the notations are defined as
% Nt                        # of transmit antennas
% Nr                        # of receive antennas
% T                         OSTBC code length (so OSTBC code C(.) has dimension Nt by T)
% P                         OSTBC block numbers (so the coherence time is assumed to be PT samples)
% Yp (Nr X T)         received code matrix at block p
% H (Nr X Nt)         MIMO channel matrix, with each coefficient being
%                             i.i.d. complex Gaussian with mean 0 and variance 1
% sp (K X 1)          transmitted QAM symbol vector with K being the
%                             symbols per block
% C(.)  (Nt X T)      OSTBC function which maps up to an OSTBC block,
% and simulate the performance of the LFSDR-based blind ML OSTBC detector. 
%
% For simplicity, the 3 by 4 OSTBC (Nt=3, T=4, K=3) is considered [7].
%        | u1 -u2' -u3' 0   |
%        | u2  u1'  0  -u3' |  
%        | u3   0   u1' u2' |
% where up=[u1, u2, u3] and superscrip " ' " denotes the complex conjugate
%
% INPUT:: The parameters that the user can input are
%  -Nr is the # of receive antennas
%  -P  is the OSTBC block numbers
%  -trial is the # of simulation trials
%  -rand_num is the # of randomizations (see Sec. III-C))
%  -SNRt  is the signal-to-noise ratio range
%  -QAM  is the signal constellation of sp. This program provides two choices:
%            16-QAM (QAM=16) and 64-QAM (QAM=64)
%
% OUTPUT: When complete, the program will output the averaged symbol error
% rates of the proposed LFSDR detectors using both principle eigenvector
%  and randomization based solution approximation methods.
% The SER of the coherent ML detector (which has the perfect knowledge of
% H) will also be output. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all;
%------------------------- parameter settings -----------------------------
Nr=4;    % number of receive antenna
P=4;     % number of code block   %{{when run SD, set P=2}}%
trial=5e3;     % number of independent runs
rand_num=30;  % length of randomization
SNRt=[1:3:16];  % signal-to-noise ratio range
QAM=16;  %signal constellation  (16 for 16-QAM and 64 for 64-QAM)

%------------------------fixed parameters ----------------------------------
Nt=3;    % number of transmit antenna
T=4;     % OSTBC code length
K=2*3;   % no. of PAM symbols per code

% signal constellation
if QAM==16
PAMQ=[-3,-1,1,3];
eta=9;
else 
PAMQ=[-7,-5,-3,-1,1,3,5,7];
eta=49;
end

% basis matrices for complex 3x4 ostbc code 
X1 = [1 0 0 0;0 1 0 0;0 0 1 0];
X2 = [0 -1 0 0;1 0 0 0;0 0 0 1];
X3 = [0 0 -1 0;0 0 0 -1;1 0 0 0];
X4 = [j 0 0 0;0 -j 0 0;0 0 -j 0];
X5 = [0 j 0 0;j 0 0 0;0 0 0 -j];
X6 = [0 0 j 0;0 0 0 j;j 0 0 0];

%---------------------simulation loop--------------------------------------

for TT=1:trial
    for LL=1:length(SNRt)


%--------------signal model generation-------------------------
SNR=SNRt(LL); 
[G,F,Yp,H, NPW, sp,up] = model_generation(Nt,Nr,T,P,K,SNR,PAMQ,TT,LL,QAM);     
% the model_generation function outputs the matrix G in (20), which will be
% used in the following blind ML detection

%-------------proposed LFSDR based blind ML detection (the developed specialied IPA will be called)-----------
% build the initial feasible point for the specialized IPA (see Section
% III-B and Table I for the details)
n = length(up);
A = [ones(1,n-1) up(1,1)^2];
X = [5*eye(n-1) zeros(n-1,1);zeros(1,n-1) 1];
X = X/(A*diag(X)); 
z = (1.1)*abs(G)*ones(n,1);
t2 = 2*z(1:n-1,1);
t1 = z(1:n-1,1);
ld = [t1;t2];
v = (z(n,1)+ones(1,n-1)*(eta*t2-t1))/(up(1,1)^2);
vv=[v.' ld.']';   

% parameter setting
tol = 1e-3;  % stopping accuracy 
max_iter =100; % maximum # of iterations

% call the specialized IPA 
[X_opt,v_hat,psi,iter]= specialized_ipa(G,up(1,1)^2,eta,X,vv,tol,max_iter); % call the specialized IPA
X_hat=X_opt/X_opt(n,n); %(see (25))

%------------solution approximation of LFSDR (see Section III-C)-----------
% principal eigenvector based solution approxmation
[up_eigen] =  eigen_decomp(X_hat,up,n,PAMQ); % extract the principal eigenvector   
sp_eigen=qam_construct(up_eigen, P); % quantization to QAM signals

% randomization based solution approximation
[up_random] = randomization(X_hat,G,n,up,rand_num,PAMQ); % Gaussian randomization
sp_random=qam_construct(up_random, P);
       
%------------coherent ML detection [7]-----------
up_ml=ML_CSI(up,Yp,H,Nt,T,K,P,PAMQ);
sp_ml=qam_construct(up_ml, P);

%------------error counting-----------
error_eigen(TT,LL) = sum(sp~=sp_eigen);
error_random(TT,LL) = sum(sp~=sp_random);
error_ml(TT,LL) = sum(sp~=sp_ml);

    end
end


%-------------------------output average SER results-------------------------
% list averaged SERs
SER_ml = sum(error_ml,1)/trial./P/3
SER_eigen = sum(error_eigen,1)/trial./P/3
SER_random = sum(error_random,1)/trial./P/3
% plot figure
semilogy(SNRt,SER_eigen,'o-g',...
    SNRt,SER_random,'o-c',SNRt,SER_ml,'-')
xlabel('SNR (dB)');ylabel('SER');
legend('LF-SDR blind ML (principal eigenvector)','LF-SDR blind ML (randomization)','ML with CSI')























