%function BER2m_e = ostbc2m_e(M, frLen, numPackets, EbNo, pLen)
%OSTBC2M_E  Orthogonal space-time block coding with channel estimation for
%   2xM antenna configurations. 
%   Suggested parameter values:
%       M = 1 or 2; FRLEN = 100; NUMPACKETS = 1000; EBNOVEC = 0:2:20, PLEN = 8;
%
%   Example:
%       ber22_e = ostbc2m_e(2, 100, 1000, 0:2:20, 8);

%% Simulation parameters
clear all;
close all;
N = 2;              % Number of transmit antennas
rate = 1; inc = N/rate;
M = 1; frLen = 64; numPackets = 1500; EbNo = 0:2:24, pLen = 64; cp = 16; L = 18;%L-taps of the channel

h = waitbar(0, 'Percentage Completed');
set(h, 'name', 'Please wait...');
wb = 100/(length(EbNo)*4);

% Create  mod-demod objects
for P = [2 4 16 64]
%P = 16; % modulation order
spmod = modem.qammod('M', P, 'SymbolOrder', 'Gray', 'InputType', 'Integer');
spdemod = modem.qamdemod(spmod);
% spmod = modem.qammod(P);
% spdemod = modem.qamdemod(spmod);


% Pilot sequences - orthogonal set over N
%function [t_preamble_cp,f_preamble] = preamble_gen(GuardInterval,N_STS)
[pilots_a pilots pilots_f N_LTF] =  preamble_gen(2);

%%  Pre-allocate variables for speed
txEnc(1:52, 1:2) = zeros(52, 2);
r = zeros(pLen + frLen/rate, M);
error2m_e = zeros(1, numPackets); 
z_e = zeros(1,52);
%BER2m_e = zeros(1, length(EbNo)); 
FFT_64 = fft(eye(64))/sqrt(N);
FFT_L = FFT_64(:, 1:L);

%% Loop over EbNo points
for idx = 1:length(EbNo)
    % Loop over the number of packets
    for packetIdx = 1:numPackets
        data = randint(frLen-12, 1, P);         % data vector per user/channel
        %data = ones(frLen-12, 1);              %Just for Test
        tx = modulate(spmod, data);        % M-PSK modulation
        
        % Alamouti ST-Block Encoder, G2, full rate
        % G2 = [s1 s2; -s2* s1*]
        s1 = tx(1:N:end); s2 = tx(2:N:end);
        x1 = [s1 s2];
        txEnc(1:2:51, 1:2) = [s1 s2];
        txEnc(2:2:52, 1:2) = [-conj(s2) conj(s1)];
        txEnc = txEnc.'; %'
        
        %Add pilots
        txEnc_p(1, 1:56) = [ txEnc(1, 1:7) 1 txEnc(1, 8:20) 1 txEnc(1, 21:32) 1 txEnc(1, 33:45) -1 txEnc(1, 46:52) ];
        txEnc_p(2, 1:56) = [ txEnc(2, 1:7) 1 txEnc(2, 8:20) -1 txEnc(2, 21:32) -1 txEnc(2, 33:45) 1 txEnc(2, 46:52) ];
        
        %IFFT
        txEnc_t(1, :) = ifft([0 txEnc_p(1,29:56) zeros(1, 7),  txEnc_p(1,1:28)], frLen)*sqrt(frLen);
        txEnc_t(2, :) = ifft([0 txEnc_p(2,29:56) zeros(1, 7),  txEnc_p(2,1:28)], frLen)*sqrt(frLen);
        txEnc_t_c = [ txEnc_t(:, (49:64)) txEnc_t(:,:) ];
        
        % Prepend pilot symbols for each frame
        transmit = [pilots txEnc_t_c];

        % Create the MultiPath channel response matrix
        % Received signal for each Rx antenna with pilot

		[r_t_sum_p,r_f_sum]=channel(transmit, 3, 2, 2,64,64,(EbNo(idx)+10*log10(log2(P))));

		r_f_sum_p = r_f_sum(1:128);
		r_f = r_f_sum(129:192);
			
           
        % Channel Estimation(LS)
        %   For each link => N*M estimates
        %         for n = 1:N
        %             H_e(1, n, :) = (r(1:pLen, :).' * pilots(:, n))./pLen;
        %         end
        %         %   held constant for the whole frame
        %         H_e = H_e(ones(frLen/rate, 1), :, :);
        %[h_t mse mmse] = channel_ls(t_preamble_cp,f_preamble,n_tx_antenna,switch_awgn_rayleigh,snr,M,L,fd)
        %ls估計
        %Generate N-FFT Matrix
                %Generate A Matrix
        A_1 = [ diag(pilots_f(1,1:64))*FFT_L     diag(pilots_f(2,1:64))*FFT_L  ];
        AH_1 = conj(A_1.'); 
        h1_1 = inv(AH_1*A_1);
        h2_1 = h1_1*AH_1;
        r_f_sum_1 = r_f_sum_p(1:64);
         h_t_1 = h2_1*r_f_sum_1.';
    
        A_2 = [ diag(pilots_f(1,65:128))*FFT_L   diag(pilots_f(2,65:128))*FFT_L  ];
        AH_2 = conj(A_2.'); 
        h1_2 = inv(AH_2*A_2);
        h2_2 = h1_2*AH_2;
        r_f_sum_2 = r_f_sum_p(65:128);
        h_t_2 = h2_2*r_f_sum_2.';
    
        h_t = (h_t_1 + h_t_2)/2;
   
        h_t_padding0(1,:) = [h_t(1    :L  ).',zeros(1,frLen-L)];
        h_t_padding0(2,:) = [h_t(L+1  :2*L).',zeros(1,frLen-L)];
        stem(abs(h_t_padding0(1,:)));
        
        for index_tx_antenna = 1:N,  %receive index_tx_antenna'th antenna frenquency response
            h_f_padding0(index_tx_antenna,:) = fft(h_t_padding0(index_tx_antenna,:),frLen)/sqrt(frLen);
        end

        %Adjust the frenquency response to make them comply the order of
        %the S-T decoder
        h_f_padding0_order(1, :) = [h_f_padding0(1,37:64) h_f_padding0(1,2:29)];
        h_f_padding0_order(2, :) = [h_f_padding0(2,37:64) h_f_padding0(2,2:29)];
        
        %Adjust the frenquency response to make them comply the order after
        %remove pilots
        h_f_padding0_order_rp(1, :) = [h_f_padding0_order(1, 1:7) h_f_padding0_order(1, 9:21) h_f_padding0_order(1, 23:34) h_f_padding0_order(1, 36:48) h_f_padding0_order(1, 50:56)];
        h_f_padding0_order_rp(2, :) = [h_f_padding0_order(2, 1:7) h_f_padding0_order(2, 9:21) h_f_padding0_order(2, 23:34) h_f_padding0_order(2, 36:48) h_f_padding0_order(2, 50:56)];
        
       
        r_f_order = [r_f(37:64) r_f(2:29)];
        r_f_order_rp = [r_f_order(1:7) r_f_order(9:21) r_f_order(23:34) r_f_order(36:48) r_f_order(50:56)];
        
        %Alamouti Space-Time Decoder
        z1_e = (r_f_order_rp(1:N:end).* conj(h_f_padding0_order_rp(1,1:N:end)) + conj(r_f_order_rp(2:N:end)).* h_f_padding0_order_rp(2,2:N:end))./(abs(h_f_padding0_order_rp(1,1:N:end).^2)+(abs(h_f_padding0_order_rp(2,2:N:end)).^2));          
        z2_e = (r_f_order_rp(1:N:end).* conj(h_f_padding0_order_rp(2,1:N:end)) - conj(r_f_order_rp(2:N:end)).* h_f_padding0_order_rp(1,2:N:end))./(abs(h_f_padding0_order_rp(2,1:N:end).^2)+(abs(h_f_padding0_order_rp(1,2:N:end)).^2));
        
        z_e(1:N:end) = z1_e/8*sqrt(2); 
        z_e(2:N:end) = z2_e/8*sqrt(2);
        % ML Detector (minimum Euclidean distance)
        demod2m_e = demodulate(spdemod, z_e); % estimated

        % Determine bit errors
        error2m_e(packetIdx) = biterr(demod2m_e.', data);
        error2m_e_s(packetIdx) = symerr(demod2m_e.', data);
    end % end of FOR loop for numPackets

    % Calculate BER for current idx
    BER2m_e(P,idx) = sum(error2m_e)/(numPackets*(frLen-12)*log2(P)); % G2 coded
    SER2m_e(P,idx) = sum(error2m_e_s)/(numPackets*(frLen-12));

    str_bar = [num2str(wb) '% Completed'];
    waitbar(wb/100, h, str_bar);
    wb = wb + 100/(length(EbNo)*4);
end  % end of for loop for EbNo
end
close(h);
semilogy(EbNo, BER2m_e(2,:), '-*', EbNo, BER2m_e(4,:), '-+', EbNo, BER2m_e(16,:), '-.', EbNo, BER2m_e(64,:), '-*');
title('QAM 2*1 Alamouti STBC LS-CE-P BER performance');
xlabel('SNR(dB)'); ylabel('BER');
legend('2PSK', 'QPSK', '16-QAM', '64-QAM');
grid on;
eval(['print -djpeg -r300 ', 'result_ber_perfect_18.jpg;']);   
close;
semilogy(EbNo, SER2m_e(2,:), '-*', EbNo, SER2m_e(4,:), '-+', EbNo, SER2m_e(16,:), '-.', EbNo, SER2m_e(64,:), '-*');
title('QAM 2*1 Alamouti STBC LS-CE-P SER performance');
xlabel('SNR(dB)'); ylabel('SER');
legend('2PSK', 'QPSK', '16-QAM', '64-QAM');
grid on;
eval(['print -djpeg -r300 ', 'result_ser_perfect_18.jpg;']);   
close;
% [EOF]