function [G,F,Yp,H, NPW, sp,up] =model_generation(Nt, Nr,T,P,K,SNR,PAMQ,TT,LL,QAM)
% == Nov. 2009 ==
% This function generates 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
%
% OUTPUT: This function returns the following items
% -G   is the matrix in (20)  
% -F   is the matrix in (14)
% -Yp is the received code matrix
% -H   is the MIMO channel matrix
% -NPW is the AWGN power
% -sp  is the transmitted QAM symbol vector
% -up  is transmitted PAM symbol vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----QAM signal generation----
rand('seed',123*TT+98); 
[up] = quan_PAM_gen(Nt,K,P,PAMQ); % 8-PAM generation (data source)
sp = qam_construct(up, P);
        
%----MIMO channel generation----
randn('seed',7621*TT+43); 
H = sqrt(1/2)*(randn(Nr,Nt)+j*randn(Nr,Nt));  % Rayleigh fading channel

%----AWGN generation----
if QAM==16
        NPW = (10*Nt*K/T)*10^(-SNR/10);  % noise power
else
        NPW = (42*Nt*K/T)*10^(-SNR/10);  % noise power
end
randn('seed',7621*TT+4399*LL+674); 
N = sqrt(NPW/2)*(randn(Nr,P*T)+j*randn(Nr,P*T));  % AWGN

%----OSTBC generation (3 by 4 OSTBC (Nt=3, T=4, K=3))----
C = zeros(Nt,P*T);
 for ii=1:P
       C(:,4*ii-3:4*ii)=[sp(3*ii-2,1) -sp(3*ii-1,1)' -sp(3*ii,1)'      0     ;
        sp(3*ii-1,1)  sp(3*ii-2,1)'      0     -sp(3*ii,1)' ;
        sp(3*ii,1)          0       sp(3*ii-2,1)' sp(3*ii-1,1)'];
 end % end of OSTBC mapping
        
%----Received signal----
Yp=H*C+N; 
        
%----Construction of matrix F in (14)----        
       %--- 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];
        %
        F=zeros(P*T,P*T);   
        for f=1:P
            for k=1:P
                F(6*f-5:6*f,6*k-5:6*k)=[real(trace(Yp(:,4*f-3:4*f)*X1'*X1*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X1'*X2*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X1'*X3*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X1'*X4*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X1'*X5*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X1'*X6*Yp(:,4*k-3:4*k)'));
                    real(trace(Yp(:,4*f-3:4*f)*X2'*X1*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X2'*X2*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X2'*X3*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X2'*X4*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X2'*X5*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X2'*X6*Yp(:,4*k-3:4*k)'));
                    real(trace(Yp(:,4*f-3:4*f)*X3'*X1*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X3'*X2*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X3'*X3*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X3'*X4*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X3'*X5*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X3'*X6*Yp(:,4*k-3:4*k)'));
                    real(trace(Yp(:,4*f-3:4*f)*X4'*X1*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X4'*X2*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X4'*X3*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X4'*X4*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X4'*X5*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X4'*X6*Yp(:,4*k-3:4*k)'));
                    real(trace(Yp(:,4*f-3:4*f)*X5'*X1*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X5'*X2*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X5'*X3*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X5'*X4*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X5'*X5*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X5'*X6*Yp(:,4*k-3:4*k)'));
                    real(trace(Yp(:,4*f-3:4*f)*X6'*X1*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X6'*X2*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X6'*X3*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X6'*X4*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X6'*X5*Yp(:,4*k-3:4*k)')) real(trace(Yp(:,4*f-3:4*f)*X6'*X6*Yp(:,4*k-3:4*k)'))];
            end
        end
        
%----Construction of matrix G in (20)----        
        R = F(2:P*K,2:P*K);p=F(1,1);  
        q = F(2:P*K,1);
        G = [R up(1,1)*q;up(1,1)*q' up(1,1)^2*p]; % assume up_1 is known
        
        
        