% STBC 2*2
% wiht estimation and without estimation
frmLen = 100; % frame length
maxNumErrs = 300; % maximum number of errors
maxNumPackets = 3000; % maximum number of packets
EbNo = 0:2:12; % Eb/No varying to 12 dB
N = 2;            % number of Tx antennas
M = 2;            % number of Rx antennas
pLen = 20;         % number of pilot symbols per frame
W = hadamard(pLen);
pilots = W(:,1:N); % orthogonal set per transmit antenna

P = 2; % modulation order

% Pre-allocate variables for speed
tx2 = zeros(frmLen, N);
r = zeros(pLen + frmLen, M);
H = zeros(pLen + frmLen, N, M);
H_e = zeros(frmLen,N, M);

z1_e = zeros(frmLen/N, M); 
z2_e = z1_e;

z_e = zeros(frmLen, M); 
z = z_e; 
z1 = z1_e; 
z2 = z2_e;
BER22_e = zeros(1, length(EbNo));
BER22 = BER22_e;

% 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');


% Loop over several EbNo points
for idx = 1:length(EbNo)
    numPackets = 0;
    totNumErr22 = 0; 
    totNumErr22_e = 0;
    
    % Loop till the number of errors exceed 'maxNumErrs'
    % or the maximum number of packets have been simulated
    
    while (totNumErr22 < maxNumErrs) && (totNumErr22_e < maxNumErrs) && (numPackets < maxNumPackets)
        
        data = randint(frmLen, 1, P); % data vector per user per channel
        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)];
        
        % Prepend pilot symbols for each frame
        transmit = [pilots; tx2];
        
        % Create the Rayleigh distributed channel response matrix
        H(1, :, :) = (randn(N, M) + j*randn(N, M))/sqrt(2);
        % assume held constant for the whole frame and pilot symbols
        H = H(ones(pLen + frmLen, 1), :, :);
      
        % Received signal for each Rx antenna
        % with pilot symbols transmitted 
        for i = 1:M
        % with normalized Tx power
           r(:, i) = awgn(sum(H(:, :, i).*transmit, 2)/sqrt(N), EbNo(idx));
        end
        
        
        % Channel Estimation
        % For each link => N*M estimates
        for n = 1:N
           H_e(1, n, :) = (r(1:pLen, :).' * pilots(:, n))./pLen;
        end
       
        % assume held constant for the whole frame
        H_e = H_e(ones(frmLen, 1), :, :);

       % Combiner using estimated channel
        heidx = 1:N:length(H_e);
        for i = 1:M
                z1_e(:, i) = r(pLen+1:N:end, i).* conj(H_e(heidx, 1, i)) + ...
                conj(r(pLen+2:N:end, i)).* H_e(heidx, 2, i);
                z2_e(:, i) = r(pLen+1:N:end, i).* conj(H_e(heidx, 2, i)) - ...
                conj(r(pLen+2:N:end, i)).* H_e(heidx, 1, i);
        end
       z_e(1:N:end, :) = z1_e; 
       z_e(2:N:end, :) = z2_e;

       
   % Combiner using known channel
    hidx = pLen+1:N:length(H);
    for i = 1:M
        z1(:, i) = r(pLen+1:N:end, i).* conj(H(hidx, 1, i)) + ...
                   conj(r(pLen+2:N:end, i)).* H(hidx, 2, i);
        z2(:, i) = r(pLen+1:N:end, i).* conj(H(hidx, 2, i)) - ...
                   conj(r(pLen+2:N:end, i)).* H(hidx, 1, i);
    end
    z(1:N:end, :) = z1;
    z(2:N:end, :) = z2;
    
    % ML Detector (minimum Euclidean distance)
   demod22_e = pskdemod(sum(z_e, 2),2); % estimated
   demod22 = pskdemod(sum(z, 2),2); % known
   
   % Determine errors
   numPackets = numPackets + 1;
   totNumErr22_e = totNumErr22_e + biterr(demod22_e, data);
   totNumErr22 = totNumErr22 + biterr(demod22, data);
   
   end % end of FOR loop for numPackets

    % Calculate BER for current idx
    % for estimated channel
      BER22_e(idx) = totNumErr22_e/(numPackets*frmLen);
    % for known channel
      BER22(idx) = totNumErr22/(numPackets*frmLen);
      
      % Plot results
    semilogy(EbNo(1:idx), BER22_e(1:idx), 'ro');
    semilogy(EbNo(1:idx), BER22(1:idx), 'g*');
    legend(['Channel estimated with ' num2str(pLen) ' pilot symbols/frame'],...
    'Known channel');
    drawnow;
end % end of for loop for EbNo

% Perform curve fitting and replot the results
fitBER22_e = berfit(EbNo, BER22_e);
fitBER22 = berfit(EbNo, BER22);
semilogy(EbNo, fitBER22_e, 'r', EbNo, fitBER22, 'g'); hold off;