frmLen = 100;              % frame length
numPackets = 1000;         % number of packets
EbNo = 0:2:12;             % Eb/No varying to 20 dB
N = 2;                     % maximum number of Tx antennas
M = 2;                     % maximum number of Rx antennas

%  and set up the simulation

% Create BPSK mod-demod objects
P = 2;                % modulation order

% Pre-allocate variables for speed
% Pre-allocate variables for speed
tx2 = zeros(frmLen, N);        %發(fā)射：100*N 
H   =  zeros(frmLen, N, M);       %信道：100*N*M

r11 = zeros(frmLen, 1);        %  uncoded 1x1 system

r12 = zeros(frmLen, M);        %  % combined 1xM system

r21 = zeros(frmLen, N);        % coded Nx1 system

r22=zeros(frmLen, M);


z21 = zeros(frmLen, 1); 
z21_1 = zeros(frmLen/N, 1); 
z21_2 = z21_1;
z12 = zeros(frmLen, M);

z_e = zeros(frmLen, M); 
z = z_e; 
z1 = zeros(frmLen/N, M);; 
z2 = zeros(frmLen/N, M);;

error11 = zeros(1, numPackets);
BER11 = zeros(1, length(EbNo));

error21 = error11; 
BER21 = BER11; 

error12 = error11; 
BER12 = BER11;

error22=error11;
BER22 = BER11;

% Set up a figure for visualizing BER results
h = gcf; 
clf(h); grid on; hold on;
set(gca,'yscale','log','xlim',[EbNo(1), EbNo(end)],'ylim',[1e-5 1]);
xlabel('Eb/No (dB)'); ylabel('BER'); set(h,'NumberTitle','off');
set(h,'Name','Orthogonal Space-Time Block Coding' );
set(h, 'renderer', 'zbuffer'); title('G2-coded 2x2 System');


for idx = 1:length(EbNo)   
   % Loop over the number of packets
          for packetIdx = 1:numPackets
              data = randint(frmLen, 1, P);  % data vector per user per channel,frmLen*1
             
              tx = pskmod(data,P);    % BPSK modulation
             % Alamouti Space-Time Block Encoder, G2, full rate
                    %  |    Antenna 1          2     
                    %  |  
                    %  |Time 0     s1          s2
                    %  | 
                    %  |Time 1     -s2*        s1*
                    
             s1 = tx(1:N:end);      s2 = tx(2:N:end); 
             
             tx2(1:N:end, :) = [ s1       s2];       
             tx2(2:N:end, :) = [-conj(s2) conj(s1)];  
             
             % Create the Rayleigh distributed channel response matrix for two transmit and two receive antennas
             H(1:N:end, :, :) = (randn(frmLen/2, N, M) +j*randn(frmLen/2, N, M))/sqrt(2);
             % assume held constant for 2 symbol periods
             H(2:N:end, :, :) = H(1:N:end, :, :);
             
         %----------------Received signals----------------------------
             % Received signals for uncoded 1x1 system
             % Y = AWGN(X,SNR) adds white Gaussian noise to X.  The SNR is in dB.  
               r11 = awgn(H(:, 1, 1).*tx, EbNo(idx));
                
             % for G2-coded 2x1 system - with normalized Tx power, i.e.,
             % the total transmitted power is assumed constant
               r21 = awgn(sum(H(:, :, 1).*tx2, 2)/sqrt(N), EbNo(idx));
               
             % for Maximal-ratio combined 1x2 system
            for i = 1:M
               r12(:, i) = awgn(H(:, 1, i).*tx, EbNo(idx));
            end
             % for G2-coded 2x2 system - with normalized Tx power, i.e.,
           
            for i = 1:M
            % with normalized Tx power
              r22(:, i) = awgn(sum(H(:, :, i).*tx2, 2)/sqrt(N), EbNo(idx));
           end
            
            % Front-end Combiners - assume channel response known at Rx
            % for G2-coded 2x1 system
            hidx = 1:N:length(H);
                z21_1 = r21(1:N:end).* conj(H(hidx, 1, 1)) + ...
                              conj(r21(2:N:end)).* H(hidx, 2, 1);
                z21_2 = r21(1:N:end).* conj(H(hidx, 2, 1)) - ...
                              conj(r21(2:N:end)).* H(hidx, 1, 1);
                z21(1:N:end) = z21_1;
                z21(2:N:end) = z21_2;
            
            % for Maximal-ratio combined 1x2 system
            for i = 1:M
                z12(:, i) = r12(:, i).* conj(H(:, 1, i));
            end

            
              % Combiner using known channel
            hidx = 1:N:length(H);
            for i = 1:M
                z1(:, i) = r22(1:N:end, i).* conj(H(hidx, 1, i)) + ...
                           conj(r22(2:N:end, i)).* H(hidx, 2, i);
                z2(:, i) = r22(1:N:end, i).* conj(H(hidx, 2, i)) - ...
                           conj(r22(2:N:end, i)).* H(hidx, 1, i);
            end
            z(1:N:end, :) = z1;
            z(2:N:end, :) = z2;
            
            
            % ML Detector (minimum Euclidean distance)
            demod11 = pskdemod(r11.*conj(H(:, 1, 1)),P);
            demod21 = pskdemod(z21,P);
            demod12 = pskdemod(sum(z12, 2),P);
            demod22=  pskdemod(sum(z, 2),P); % estimated
            % Determine errors
            error11(packetIdx) = biterr(demod11, data);
            error21(packetIdx) = biterr(demod21, data);
            error12(packetIdx) = biterr(demod12, data);
            error22(packetIdx) = biterr(demod22, data);
     end % end of FOR loop for numPackets
     
% Calculate BER for current idx
% for uncoded 1x1 system
           BER11(idx) = sum(error11)/(numPackets*frmLen);
% for G2 coded 2x1 system
           BER21(idx) = sum(error21)/(numPackets*frmLen);
% for Maximal-ratio combined 1x2 system
           BER12(idx) = sum(error12)/(numPackets*frmLen);
% for  combined 2x2 system
           BER22(idx) = sum(error22)/(numPackets*frmLen);
% Plot results
        semilogy(EbNo(1:idx), BER11(1:idx), 'r*', ...
        EbNo(1:idx), BER21(1:idx), 'go',...
        EbNo(1:idx), BER12(1:idx), 'bs',...
        EbNo(1:idx), BER22(1:idx),'kd');
        legend('No Diversity (1Tx, 1Rx)', 'Alamouti (2Tx, 1Rx)',...
        'Maximal-Ratio Combining (1Tx, 2Rx)','Alamouti (2Tx, 2Rx)' );
        drawnow;
end % end of for loop for EbNo
        % Perform curve fitting and replot the results
        fitBER11 = berfit(EbNo, BER11);
        fitBER21 = berfit(EbNo, BER21);
        fitBER12 = berfit(EbNo, BER12);
        fitBER22 = berfit(EbNo, BER22);
        semilogy(EbNo, fitBER11, 'r', EbNo, fitBER21, 'g', EbNo, fitBER12, 'b',EbNo, fitBER22, 'k');
 hold off;