stbcxindaoguji.m

STBC_ofdm的最大似然信道估计算法实现（ML）

%for ML channel 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 = 1;                  % number of Rx antennas
pLen = 8;               % number of pilot symbols per frame
W = hadamard(pLen);
pilots = W(:, 1:N);     % orthogonal set per transmit antenna
 %Seed states for repeatability
seed = [98765 12345];
randn('state', seed(1)); 
rand('state', seed(2));

% 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);
z_e = zeros(frmLen, M); 
z1_e = zeros(frmLen/N, M); 
z2_e = z1_e;
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
%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('STBC 2x1 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 = 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)];

        % 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 = demodulate(bpskdemod, sum(z_e, 2)); % estimated
        demod22   = demodulate(bpskdemod, sum(z, 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;