%------------------------------------------------------------------------
%Author: Hoang Nguyen                        hoangnmt@yahoo.com
%------------------------------------------------------------------------

%STBC(Num_Frame,Num_TxAnt,Num_RxAnt,Path_Delay,Path_Gain,v,M);
function STBC(Num_Frame,Nt,Nr,M)

if Nt == 2
    switch M
        case 1
            O=[1 2;-2+j 1+j]; % Choose for real digital modulator
        case 2
            O=[1 2;-2 1];     % Choose for complex digital modulator
    end
else
    O=[1 -2 -3;2+j 1+j 0;3+j 0 1+j;0 -3+j 2+j];  %Complex or Real Orthogonal Matrix **define this**  
end
Nt=size(O,2);                                                        %Number of Transmit antennas
co_time=size(O,1);                                                   %Block time length
Nr=1;                                                                %Number of Receive antennas        **define this**
Nit=Num_Frame*10*9*log2(M);                                          %Number of repeates for each snr   **define this**
M_psk=M;                                                             %M-PSK  constellation,M_psk=2^k    
snr_min=1;                                                           %Min snr range for simulation      
snr_max=20;
%SNR=snr_min:snr_max;
%Max snr rande for simulation      
graph_inf_bit=zeros(snr_max-snr_min+1,2);                            %Plot information
graph_inf_sym=zeros(snr_max-snr_min+1,2);                            %Plot information
num_X=1;
num_bit_per_sym=log2(M_psk);
for cc_ro=1:co_time
    for cc_co=1:Nt
        num_X=max(num_X,abs(real(O(cc_ro,cc_co))));
    end
end
co_x=zeros(num_X,1);
for con_ro=1:co_time                                                 %Compute delta,epsilon,eta and conj matrices
    for con_co=1:Nt
        if abs(real(O(con_ro,con_co)))~=0
            delta(con_ro,abs(real(O(con_ro,con_co))))=sign(real(O(con_ro,con_co)));
            epsilon(con_ro,abs(real(O(con_ro,con_co))))=con_co;
            co_x(abs(real(O(con_ro,con_co))),1)=co_x(abs(real(O(con_ro,con_co))),1)+1;
            eta(abs(real(O(con_ro,con_co))),co_x(abs(real(O(con_ro,con_co))),1))=con_ro;
            coj_mt(con_ro,abs(real(O(con_ro,con_co))))=imag(O(con_ro,con_co));
        end
    end
end
eta=eta.';                                                           %Sort is not necessary
eta=sort(eta);
eta=eta.';
for SNR=snr_min:snr_max                                              %Start simulation
    clc
    %disp('Wait until SNR=');disp(snr_max);
    %SNR
    n_err_sym=0;
    n_err_bit=0;
    graph_inf_sym(SNR-snr_min+1,1)=SNR;
    graph_inf_bit(SNR-snr_min+1,1)=SNR;
    for con_sym=1:Nit
        bi_data=randint(num_X,num_bit_per_sym);                      %Random binary data
        de_data=bi2de(bi_data);                                      %Convert binary data to decimal for use in M-PSK mod.
        data=pskmod(de_data,M_psk,0);
        H=randn(Nt,Nr)+j*randn(Nt,Nr);                               %Path gains matrix 
        XX=zeros(co_time,Nt);
        for con_r=1:co_time                                          %Start space time coding
            for con_c=1:Nt
                if abs(real(O(con_r,con_c)))~=0
                    if imag(O(con_r,con_c))==0
                        XX(con_r,con_c)=data(abs(real(O(con_r,con_c))),1)*sign(real(O(con_r,con_c)));
                    else
                        XX(con_r,con_c)=conj(data(abs(real(O(con_r,con_c))),1))*sign(real(O(con_r,con_c)));
                    end
                end
            end
        end                                                          %End space time coding
        H=H.';
        XX=XX.';
        snr=10^(SNR/10);
        Noise=(randn(Nr,co_time)+j*randn(Nr,co_time));               %Generate Noise
        Y=(sqrt(snr/Nt)*H*XX+Noise).';                               %Received signal
        H=H.';                               
        %Start decoding with perfect channel estimation
        for co_ii=1:num_X
            for co_tt=1:size(eta,2)
                if eta(co_ii,co_tt)~=0
                    if coj_mt(eta(co_ii,co_tt),co_ii)==0
                        r_til(eta(co_ii,co_tt),:,co_ii)=Y(eta(co_ii,co_tt),:);
                        a_til(eta(co_ii,co_tt),:,co_ii)=conj(H(epsilon(eta(co_ii,co_tt),co_ii),:));
                    else
                        r_til(eta(co_ii,co_tt),:,co_ii)=conj(Y(eta(co_ii,co_tt),:));
                        a_til(eta(co_ii,co_tt),:,co_ii)=H(epsilon(eta(co_ii,co_tt),co_ii),:);
                    end
                end
            end
        end
        RR=zeros(num_X,1);
        for ii=1:num_X                                               %Generate decision statistics for the transmitted signal "xi"
            for tt=1:size(eta,2)
                for jj=1:Nr
                    if eta(ii,tt)~=0
                        RR(ii,1)=RR(ii,1)+r_til(eta(ii,tt),jj,ii)*a_til(eta(ii,tt),jj,ii)*delta(eta(ii,tt),ii);
                    end
                end
            end
        end
        re_met_sym=pskdemod(RR,M_psk,0);                      % = ML decision for M-PSK
        re_met_bit=de2bi(re_met_sym);
        re_met_bit(1,num_bit_per_sym+1)=0;                           %For correct demension of "re_met_bit"
        for con_dec_ro=1:num_X                                             
            if re_met_sym(con_dec_ro,1)~=de_data(con_dec_ro,1)
                n_err_sym=n_err_sym+1;
                for con_dec_co=1:num_bit_per_sym
                    if re_met_bit(con_dec_ro,con_dec_co)~=bi_data(con_dec_ro,con_dec_co)
                        n_err_bit=n_err_bit+1;
                    end
                end
            end
        end
    end
    Perr_sym=n_err_sym/(num_X*Nit);                                  %Count number of error bits and symbols 
    graph_inf_sym(SNR-snr_min+1,2)=Perr_sym;
    Perr_bit=n_err_bit/(num_X*Nit*num_bit_per_sym);
    graph_inf_bit(SNR-snr_min+1,2)=Perr_bit;
end

%Create plot
axes(handles.axes1);
%grid on
%x_sym=graph_inf_sym(:,1);                                            %Generate plot
%y_sym=graph_inf_sym(:,2);
%subplot(2,1,1);
%semilogy(x_sym,y_sym,'k-v');
%xlabel('SNR, [dB]');
%ylabel('Symbol Error Probability');
grid on
x_bit=graph_inf_bit(:,1);
y_bit=graph_inf_bit(:,2);
%subplot(2,1,2);
semilogy(x_bit,y_bit,'r-*');
xlabel('SNR, [dB]');
ylabel('Bit Error Probability');
grid on
set(handles.axes1,'XMinorTick','on')