mimo.m

THIS matlab code is STTD

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


% Seed states for repeatability
seed = [98765 12345]; randn('state', seed(1)); rand('state', seed(2));

% Create BPSK mod-demod objects
P = 2;                % modulation order
bpskmod = modem.pskmod('M', P, 'SymbolOrder', 'Gray');
bpskdemod = modem.pskdemod(bpskmod);

% Pre-allocate variables for speed
tx2 = zeros(frmLen, N); H  = zeros(frmLen, N, M);
r21 = zeros(frmLen, 1); r12  = zeros(frmLen, 2);
z21 = zeros(frmLen, 1); z21_1 = zeros(frmLen/N, 1); z21_2 = z21_1;
z12 = zeros(frmLen, M);
error11 = zeros(1, numPackets); BER11 = zeros(1, length(EbNo));
error21 = error11; BER21 = BER11; error12 = error11; BER12 = BER11;

% Set up a figure for visualizing BER results
h = gcf; 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, 'renderer', 'zbuffer');
set(h,'Name','Transmit vs. Receive Diversity');
title('Transmit vs. Receive Diversity');

% Loop over several EbNo points
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
        tx = modulate(bpskmod, data);   % BPSK modulation

        % Alamouti Space-Time Block Encoder, G2, full rate
        %   G2 = [s1 s2; -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
        %   for uncoded 1x1 system
        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

        % 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

        % ML Detector (minimum Euclidean distance)
        demod11 = demodulate(bpskdemod, r11.*conj(H(:, 1, 1)));
        demod21 = demodulate(bpskdemod, z21);
        demod12 = demodulate(bpskdemod, sum(z12, 2));

        % Determine errors
        error11(packetIdx) = biterr(demod11, data);
        error21(packetIdx) = biterr(demod21, data);
        error12(packetIdx) = biterr(demod12, 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);

    % Plot results
    semilogy(EbNo(1:idx), BER11(1:idx), 'r*', ...
             EbNo(1:idx), BER21(1:idx), 'go',...
             EbNo(1:idx), BER12(1:idx), 'bs');
    legend('No Diversity (1Tx, 1Rx)', 'Alamouti (2Tx, 1Rx)',...
           'Maximal-Ratio Combining (1Tx, 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);
semilogy(EbNo, fitBER11, 'r', EbNo, fitBER21, 'g', EbNo, fitBER12, 'b');
hold off;