RENEW Project Documentation

Version 1.0

Reconfigurable Ecosystem for Next-generation End-to-end Wireless

Basic MIMO Tutorial

In this tutorial, we go over a simple Matlab MIMO script that implements a single-shot transmission from multiple clients (also known as User Equipment or UE) to multiple base station radios. We define two different modes: a) OTA (Over-the-air) and b) SIM_MOD (simulation).

If in simulation mode we simply use a Rayleigh channel whereas the OTA mode relies on the Iris hardware for transmission and reception. In both cases the clients transmit an OFDM signal that resembles a typical 802.11 WLAN waveform. If the transmission is OTA, then the user specifies a TDD schedule that tells all clients when to transmit their frame. The base station initiates the schedule by sending a beacon signal that synchronizes clients. After that, all clients will transmit simultaneously. We implement a frame structure that allows the base station to capture clean (non-overlaping/staggered) training sequences for equalization and demultiplexing of the concurrent data streams. In the figure below we observe the structure of the transmit signals for each client. Notice this is for a four-client topology. You will notice that the training sequences of the multiple clients are staggered. This allows the base station to estimate the channel from each of its antennas to each of the clients in a clean manner. After all the training sequences, all clients send their data frames simultaneously. On the receiver side, the base station will perform channel estimation, equalization, and use zero-forcing or conjugate beamforming to demultiplex the different data streams.

NOTE: In this tutorial we assume that the user is familiar with Matlab.


close all;

% For python bindings
[version, executable, isloaded] = pyversion;
if ~isloaded
    pyversion /usr/bin/python
    py.print() %weird bug where py isn't loaded in an external script

Configure basic experiment parameters

If simulation mode, specify Rayleigh channel. If over-the-air, specify serial numbers of base station nodes and client nodes. Here we also configure the frequency, transmit and receive gains, and sampling rate.

% Params:
N_BS_NODE               = 4;
N_UE                    = 2;
WRITE_PNG_FILES         = 0;           % Enable writing plots to PNG
SIM_MOD                 = 1;
DEBUG                   = 0;
PLOT                    = 1;
    chan_type               = "rayleigh"; % Will use only Rayleigh for simulation
    sim_SNR_db              = 15;
    TX_SCALE                = 1;          % Scale for Tx waveform ([0:1])
    bs_ids                  = ones(1, N_BS_NODE);
    ue_ids                  = ones(1, N_UE);

    %Iris params:
    TX_SCALE                = 0.5;         % Scale for Tx waveform ([0:1])
    chan_type               = "iris";
    USE_HUB                 = 0;
    TX_FRQ                  = 2.5e9;
    RX_FRQ                  = TX_FRQ;
    TX_GN                   = 42;
    TX_GN_ue                = 42;
    RX_GN                   = 20;
    SMPL_RT                 = 5e6;
    N_FRM                   = 10;
    bs_ids                  = string.empty();
    ue_ids                  = string.empty();
    ue_scheds               = string.empty();
    if USE_HUB
        % Using chains of different size requires some internal
        % calibration on the BS. This functionality will be added later.
        % For now, we use only the 4-node chains:
        bs_ids = ["RF3E000134", "RF3E000191", "RF3E000171", "RF3E000105",...
            "RF3E000053", "RF3E000177", "RF3E000192", "RF3E000117",...
            "RF3E000183", "RF3E000152", "RF3E000123", "RF3E000178", "RF3E000113", "RF3E000176", "RF3E000132", "RF3E000108", ...
            "RF3E000143", "RF3E000160", "RF3E000025", "RF3E000034",...
            "RF3E000189", "RF3E000024", "RF3E000139", "RF3E000032", "RF3E000154", "RF3E000182", "RF3E000038", "RF3E000137", ...
            "RF3E000103", "RF3E000180", "RF3E000181", "RF3E000188"];
        hub_id = "FH4A000001";
        bs_ids = ["RF3E000189", "RF3E000024", "RF3E000139", "RF3E000032", "RF3E000154", "RF3E000182", "RF3E000038", "RF3E000137"];
    ue_ids= ["RF3E000060", "RF3E000157"];
    N_BS_NODE               = length(bs_ids);           % Number of nodes/antennas at the BS
    N_UE                    = length(ue_ids);           % Number of UE nodes

Configure signal parameters (OFDM)

First, select the receiver beamforming algorithm to be used. Currently we support both Zero-Forcing and Conjugate Beamforming. In addition, we specify all the necessary OFDM parameters such as number of subcarriers, numbre of OFDM symbols, modulation order, and index of pilot and data subcarriers. Notice we follow the basic parameters from a 20 MHz 802.11a waveform.

MIMO_ALG                = 'ZF';      % MIMO ALGORITHM: ZF or Conjugate 
fprintf("MIMO algorithm: %s \n",MIMO_ALG);

% Waveform params
N_OFDM_SYM              = 44;         % Number of OFDM symbols for burst, it needs to be less than 47
MOD_ORDER               = 16;          % Modulation order (2/4/16/64 = BSPK/QPSK/16-QAM/64-QAM)

% OFDM params
SC_IND_PILOTS           = [8 22 44 58];                           % Pilot subcarrier indices
SC_IND_DATA             = [2:7 9:21 23:27 39:43 45:57 59:64];     % Data subcarrier indices
SC_IND_DATA_PILOT       = [2:27 39:64]';
N_SC                    = 64;                                     % Number of subcarriers
CP_LEN                  = 16;                                     % Cyclic prefix length
N_DATA_SYMS             = N_OFDM_SYM * length(SC_IND_DATA);       % Number of data symbols (one per data-bearing subcarrier per OFDM symbol) per UE
N_LTS_SYM               = 2;                                      % Number of 
N_SYM_SAMP              = N_SC + CP_LEN;                          % Number of samples that will go over the air
N_ZPAD_PRE              = 90;                                     % Zero-padding prefix for Iris
N_ZPAD_POST             = N_ZPAD_PRE -14;                         % Zero-padding postfix for Iris

% Rx processing params
FFT_OFFSET                    = 16;          % Number of CP samples to use in FFT (on average)
DO_APPLY_PHASE_ERR_CORRECTION = 1;           % Enable Residual CFO estimation/correction

Define frame preamble sequence

We follow the 802.11 convention and define the preamble as frequency domain long-term training sequences (LTS or LTF) and convert them to time domain samples through an IFFT. Then, we generate the staggered (time-orthogonal) preambles for the different clients. These preambles are the ones that will allow the receiver to cleanly estimate the channel to each client in order to perform channel equalization.

% LTS for fine CFO and channel estimation
lts_f = [0 1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 ...
    1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1].';
lts_t = ifft(lts_f, N_SC); %time domain

% Arrange time-orthogonal pilots 
preamble_common = [lts_t(33:64); repmat(lts_t,N_LTS_SYM,1)];
l_pre = length(preamble_common);
pre_z = zeros(size(preamble_common));
preamble = zeros(N_UE * l_pre, N_UE);
for jp = 1:N_UE
    preamble((jp-1)*l_pre + 1: (jp-1)*l_pre+l_pre,jp) = preamble_common;

Generate full data signal

We generate multiple random sequences of symbols, one for each client. Then we build an array that will contain all the frequency-domain samples across subcarriers and OFDM symbols. This array combines both the pilots and data. At the end we simply construct the full time-domain OFDM waveform to be transmitted (one vector per client).

%% Generate a payload of random integers
tx_data = randi(MOD_ORDER, N_DATA_SYMS, N_UE) - 1;
tx_syms = mod_sym(tx_data, MOD_ORDER);

% Reshape the symbol vector to a matrix with one column per OFDM symbol
tx_syms_mat = reshape(tx_syms, length(SC_IND_DATA), N_OFDM_SYM, N_UE);

% Define the pilot tone values as BPSK symbols
pilots = [1 1 -1 1].';

% Repeat the pilots across all OFDM symbols
pilots_mat = repmat(pilots, 1, N_OFDM_SYM, N_UE);

% Construct the IFFT input matrix
ifft_in_mat = zeros(N_SC, N_OFDM_SYM, N_UE);

% Insert the data and pilot values; other subcarriers will remain at 0
ifft_in_mat(SC_IND_DATA, :, :)   = tx_syms_mat;
ifft_in_mat(SC_IND_PILOTS, :, :) = pilots_mat;

%Perform the IFFT
tx_payload_mat = ifft(ifft_in_mat, N_SC, 1);

% Insert the cyclic prefix
if(CP_LEN > 0)
    tx_cp = tx_payload_mat((end-CP_LEN+1 : end), :, :);
    tx_payload_mat = [tx_cp; tx_payload_mat];

% Reshape to a vector
tx_payload_vecs = reshape(tx_payload_mat, ceil(numel(tx_payload_mat)/N_UE), N_UE);

% Construct the full time-domain OFDM waveform
tx_vecs = [zeros(N_ZPAD_PRE, N_UE); preamble; tx_payload_vecs; zeros(N_ZPAD_POST, N_UE)];

% Leftover from zero padding:
tx_vecs_iris = tx_vecs;
% Scale the Tx vector to +/- 1
tx_vecs_iris = TX_SCALE .* tx_vecs_iris ./ max(abs(tx_vecs_iris));

Configure TDD schedule and create struct with all parameters

% Create BS Hub and UE objects. Note: BS object is a collection of Iris
% nodes.
bs_sched = ["BGGGGGRG"];           % BS schedule
ue_sched = ["GGGGGGPG"];           % UE schedule

%number of samples in a frame
n_samp = length(tx_vecs_iris);

% Iris nodes' parameters
bs_sdr_params = struct(...
    'id', bs_ids, ...
    'n_sdrs', N_BS_NODE, ...        % number of nodes chained together
    'txfreq', TX_FRQ, ...   
    'rxfreq', RX_FRQ, ...
    'txgain', TX_GN, ...
    'rxgain', RX_GN, ...
    'sample_rate', SMPL_RT, ...
    'n_samp', n_samp, ...          % number of samples per frame time.
    'n_frame', N_FRM, ...
    'tdd_sched', bs_sched, ...     % number of zero-paddes samples
    'n_zpad_samp', N_ZPAD_PRE ...

ue_sdr_params = bs_sdr_params; =  ue_ids;
ue_sdr_params.n_sdrs = N_UE;
ue_sdr_params.txgain = TX_GN_ue;
ue_sdr_params.tdd_sched = ue_sched;

Set up and trigger transmission and reception

If we are operating in simulation mode, we simply take the transmit vector and pass it through a Rayleigh channel. However, if we are actually dealing with hardware and transmitting signals over the air, the process is a bit more involved.

This happens inside the getRxVec() function. Below we wil go over the actions taking place inside the function.

if (SIM_MOD)
    rx_vec_iris = getRxVec(tx_vecs_iris, N_BS_NODE, N_UE, chan_type, sim_SNR_db, bs_sdr_params, ue_sdr_params, []);
    rx_vec_iris = rx_vec_iris.'; % just to agree with what the hardware spits out.
    if USE_HUB
        rx_vec_iris = getRxVec(tx_vecs_iris, N_BS_NODE, N_UE, chan_type, [], bs_sdr_params, ue_sdr_params, hub_id);
        rx_vec_iris = getRxVec(tx_vecs_iris, N_BS_NODE, N_UE, chan_type, [], bs_sdr_params, ue_sdr_params, []);

Inside getRxVec()

We first initialize both the BS nodes and UE nodes:

n_samp = bs_param.n_samp;
if ~isempty(hub_id)
    node_bs = iris_py(bs_param,hub_id);
    node_bs = iris_py(bs_param,[]);         % initialize BS
node_ue = iris_py(ue_param,[]);             % initialize UE

Then we set up the transmission and reception by setting up streams and synchronizing base station radios. The following lines (i) reset the correlator on the UE, hence, preparing it to synchronize based on the beacon it will have to detect from the base station, (ii) set up reading streams on both nodes before the transmission starts, and (iii) configure the schedules on the BS and the UE, respectively:

node_bs.sdrsync();                 % synchronize delays only for BS
node_ue.sdrrxsetup();              % set up reading stream
tdd_sched_index = 1; % for uplink only one frame schedule is sufficient
node_bs.set_tddconfig(1, bs_param.tdd_sched(tdd_sched_index)); % configure the BS: schedule, etc.
node_ue.set_tddconfig(0, ue_param.tdd_sched(tdd_sched_index)); % configure the UE: schedule, etc.

Next, we write the data to be transmitted from both the BS and UE. The BS will initiate the writing of the predefined beacon on the FPGA block ram for transmission during the ‘P’ slots of the BS schedule. The beacon sequence is hardcoded in the Python driver and the Matlab user should not be concerned with it. The UE in turn also writes its pilots onto the block ram so that they will be transmitted during the ‘P’ slots of its schedule:

node_bs.sdr_setupbeacon();                      % Burn beacon to the BS(1) RAM
for i=1:n_ue
    node_ue.sdrtx_single(tx_data(:,i), i);      % Burn data to the UE RAM

The next step is to set up reception by activating the RX stream on the base station and enabling the correlator on the UE:

node_bs.sdr_activate_rx();         % activate reading stream
node_ue.sdr_setcorr()              % activate correlator

The BS then reads and stores the received data from the receiver stream that we had set up and activated earlier:

[y, data0_len] = node_bs.sdrrx(n_samp); % read data

Finally we end the streams on both nodes:


Receiver Processing

Find start of frame and perform channel estimation

Upon reception, each base station node performs a cross-correlation in an attempt to find the LTS sent by each client. Since we know the structure of each frame, once we find the LTS we can determine the start of the payload (data).


%% Correlate for LTS
% Complex cross correlation of Rx waveform with time-domain LTS

a = 1;
unos = ones(size(preamble_common'));
lts_corr = zeros(N_BS_NODE, length(rx_vec_iris));
data_len = (N_OFDM_SYM)*(N_SC +CP_LEN);
rx_lts_mat = double.empty();
payload_ind = int32.empty();
payload_rx = zeros(N_BS_NODE, data_len);
lts_peaks = zeros(N_BS_NODE, N_UE);

for ibs =1:N_BS_NODE
        % Correlation through filtering
        v0 = filter(fliplr(preamble_common'),a,rx_vec_iris(ibs,:));
        v1 = filter(unos,a,abs(rx_vec_iris(ibs,:)).^2);
        lts_corr(ibs,:) = (abs(v0).^2)./v1; % normalized correlation
        % Sort the correlation values
        sort_corr = sort(lts_corr(ibs,:), 'descend');
        % Take the N_UE largest values
        rho_max = sort_corr(1:N_UE);
        % Get the indices of N_UE largest corr. values
        lts_peaks(ibs,:) = find(lts_corr(ibs,:) >= min(rho_max));
        % position of the last peak
        max_idx = max(lts_peaks(ibs,:));
        % In case of bad correlatons:
        if (max_idx + data_len) > length(rx_vec_iris) || (max_idx < 0) || (max_idx - length(preamble) < 0)
            fprintf('Bad correlation at antenna %d max_idx = %d \n', ibs, max_idx);
            % Real value doesn't matter since we have corrrupt data:
            max_idx = length(rx_vec_iris)-data_len -1;
        payload_ind(ibs) = max_idx +1;
        pream_ind_ibs = payload_ind(ibs) - length(preamble);
        pl_idx = payload_ind(ibs) : payload_ind(ibs) + data_len;
        rx_lts_mat(ibs,:) = rx_vec_iris(ibs, pream_ind_ibs: pream_ind_ibs + length(preamble) -1 );
        payload_rx(ibs,1:length(pl_idx) -1) = rx_vec_iris(ibs, payload_ind(ibs) : payload_ind(ibs) + length(pl_idx) -2);

Channel estimation

From the received preambles we can estimate the channel matrix between all TX and RX nodes. We convert both preambles and data to frequency domain and use our knowledge of the TX training sequence and the RX signal to do channel estimation, averaging over two LTS symbols.

% Construct a matrix from the received pilots
n_plt_samp = floor(length(preamble)/ N_UE);     % number of samples in a per-UE pilot 
Y_lts = zeros(N_BS_NODE,N_UE, n_plt_samp);
for iue = 1:N_UE
   plt_j_ix = (iue-1) * n_plt_samp +1:iue * n_plt_samp;
   Y_lts(:,iue,:) = rx_lts_mat(:,plt_j_ix);

% Take N_SC samples from each LTS
rx_lts_idx = CP_LEN +1 : N_LTS_SYM * N_SC +CP_LEN;
Y_lts = Y_lts(:,:,rx_lts_idx);

% Reshape the matix to have each lts pilot in a different dimension:
% N_BS_NODE x N_UE x 64 x 2
Y_lts = reshape(Y_lts, N_BS_NODE, N_UE, [], N_LTS_SYM); 

% Take FFT:
Y_lts_f = fft(Y_lts, N_SC,3);
% Construct known pilot matrix to use i next step:
lts_f_mat = repmat(lts_f, N_BS_NODE *N_UE * N_LTS_SYM,1);
lts_f_mat = reshape(lts_f_mat, [], N_LTS_SYM, N_UE, N_BS_NODE);
lts_f_mat = permute(lts_f_mat, [4 3 1 2]);

% Get the channel by dividing by the pilots
G_lts = Y_lts_f ./ lts_f_mat;
% Estimate the channel by averaging over the two LTS symbols:
H_hat = mean(G_lts, 4);

% Reshape the payload and take subcarriers without the CP
payload_rx = reshape(payload_rx,N_BS_NODE, (N_SC + CP_LEN),[]);
payload_noCP = payload_rx(:,CP_LEN-FFT_OFFSET+(1:N_SC),:);
% Take FFT
Y_data = fft(payload_noCP, N_SC, 2);

RX Beamforming

Compute beamforming weights (zero-forcing and conjugate supported), and apply them to separate data streams. We also compute the channel condition number.

% Apply ZF by multiplying the pseudo inverse of H_hat[N_BS_NODE x NUE] for each suubcarrier:
nz_sc = find(lts_f ~= 0); % non-zero subcarriers
syms_eq = zeros(N_UE,N_SC,N_OFDM_SYM);
channel_condition = double.empty();
channel_condition_db = double.empty();
for j=1:length(nz_sc)
       % Pseudo-inverse:(H'*H)^(-1)*H':
        HH_inv = inv((squeeze(H_hat(:,:, nz_sc(j) ) )' * squeeze(H_hat(:,:, nz_sc(j) ) ) ) ) * squeeze(H_hat(:,:, nz_sc(j) ) )';
        x = HH_inv*squeeze(Y_data(:,nz_sc(j),:));
        % Do yourselves: Conj BF: 
        % Normalization coeff:
        H_pow = diag(abs (H_hat(:,:, nz_sc(j) )' * H_hat(:,:, nz_sc(j) ) ));
        % Apply BF: 
        x = (H_hat(:,:, nz_sc(j) )') * squeeze(Y_data(:,nz_sc(j),:))./ repmat(H_pow, 1, N_OFDM_SYM);
    syms_eq(:,nz_sc(j),:) = x;
    channel_condition(nz_sc(j)) = cond(H_hat(:,:,nz_sc(j) ) );
    channel_condition_db(nz_sc(j)) = 10*log10(channel_condition(nz_sc(j)) );

Apply phase correction (if enabled)

Phase error is obtained from pilots across symbols. Once we measure the error, we apply it to entire frame.

    % Extract the pilots and calculate per-symbol phase error
    pilots_f_mat = syms_eq(:,SC_IND_PILOTS,:);
    pilots_f_mat_comp = pilots_f_mat.* permute(pilots_mat, [3 1 2]);
    pilot_phase_err = angle(mean(pilots_f_mat_comp,2));
	% Define an empty phase correction vector (used by plotting code below)
    pilot_phase_err = zeros(N_UE, 1, N_OFDM_SYM);

pilot_phase_err_corr = repmat(pilot_phase_err, 1, N_SC, 1);
pilot_phase_corr = exp(-1i*(pilot_phase_err_corr));

% Apply the pilot phase correction per symbol
syms_eq_pc = syms_eq.* pilot_phase_corr;

Take only data subcarriers and prepare the structure to be used for demodulation

syms_eq_pc = syms_eq_pc(:,SC_IND_DATA,:);

% Reshape the 3-D matrix to 2-D:
syms_eq_pc = reshape(syms_eq_pc, N_UE, [] );
syms_eq_pc = syms_eq_pc.';

Finally, call function to demodulate the RX signal

rx_data = demod_sym(syms_eq_pc ,MOD_ORDER);

Plot received samples

Finally, we plot a variety of results that can be used to understand the performance of the system:

Start with plotting setup

cf = 0;
fst_clr = [0, 0.4470, 0.7410];
sec_clr = [0.8500, 0.3250, 0.0980];

  1. The received signal from all clients at each base station antenna. You will notice the discrepacy in power from the different staggered training sequences (different clients at different locations). The signal that follows is simply the combination of all the concurrent data from the multiple clients.
% Rx signal
cf = cf + 1;
figure(cf); clf;
for sp = 1:N_BS_NODE
    subplot(N_BS_NODE,2,2*(sp -1) + 1 );
    axis([0 length(rx_vec_iris(sp,:)) -TX_SCALE TX_SCALE])
    grid on;
    title(sprintf('BS antenna %d Rx Waveform (I)', sp));

    plot(imag(rx_vec_iris(sp,:)), 'color' , sec_clr);
    axis([0 length(rx_vec_iris(sp,:)) -TX_SCALE TX_SCALE]);
    grid on;
    title(sprintf('BS antenna %d Rx Waveform (Q)', sp));
  1. Received symbol constellation for a 16-QAM constellation. We plot both the equalized and demultiplexed stream from each user, as well as the constellation seen by each base station antenna before any processing on the data.
%% Constellations
cf = cf+ 1;
figure(cf); clf;
if N_BS_NODE >=4
    sp_rows = ceil(N_BS_NODE/4)+1;
    sp_rows = ceil(N_BS_NODE/2)+1;
sp_cols = ceil(N_BS_NODE/(sp_rows -1));

for sp=1:N_UE
    subplot(sp_rows,sp_cols, sp);
    plot(syms_eq_pc(:,sp),'o','MarkerSize',1, 'color', sec_clr);
    axis square; axis(1.5*[-1 1 -1 1]);
    grid on;
    hold on;
    plot(tx_syms(:, sp),'*', 'MarkerSize',10, 'LineWidth',2, 'color', fst_clr);
    title(sprintf('Equalized Uplink Tx (blue) and Rx (red) symbols for stream %d', sp));
    % legend({'Rx','Tx'},'Location','EastOutside', 'fontsize', 12);

for sp=1:N_BS_NODE
    subplot(sp_rows,sp_cols, sp_cols+sp);
    axis square; axis(max(max(max( abs( Y_data)) ))*[-1 1 -1 1]);
    title(sprintf('Unequalized received symbols at BS ant. %d', sp));
    grid on;
    hold on;
  1. Channel magnitude obtained at each base station antenna from the client training sequences, as well as the channel 2-norm condition number for inversion. This is based on Matlab’s cond() function.
%% Channel Estimates
cf = cf + 1;
cond_clr = [0.8500, 0.250, 0.1980];

bw_span = (20/N_SC) * (-(N_SC/2):(N_SC/2 - 1)).';

figure(cf); clf;
sp = 0;
for ibs = 1:N_BS_NODE
    for iue = 1:N_UE
        sp = sp+1;
        bar(bw_span, fftshift(abs( squeeze(H_hat(ibs, iue, : ) ) ) ),1,'LineWidth', 1);
        axis([min(bw_span) max(bw_span) 0 1.1*max(abs( squeeze(H_hat(ibs, iue, :) ) ) )])
        grid on;
        title(sprintf('UE %d -> BS ant. %d Channel Estimates (Magnitude)', iue, ibs))
        xlabel('Baseband Frequency (MHz)')

bh = bar(bw_span, fftshift(channel_condition_db) ,1, 'LineWidth', 1);
axis([min(bw_span) max(bw_span) 0 max(channel_condition_db)+1])
grid on;
title('Channel Condition (dB)')
xlabel('Baseband Frequency (MHz)')
  1. Error Vector Magnitude (EVM) as a function of the data symbol index and the corresponding effective SNR. We also look at EVM as a function of both subcarrier index and OFDM symbol index.
%% EVM & SNR
sym_errs = sum(tx_data(:) ~= rx_data(:));
bit_errs = length(find(dec2bin(bitxor(tx_data(:), rx_data(:)),8) == '1'));
evm_mat = double.empty();
aevms = zeros(N_UE,1);
snr_mat = zeros(N_UE,1);

cf = cf + 1;
figure(cf); clf;
for sp = 1:N_UE
    tx_vec = tx_syms_mat(:,:,sp);
    evm_mat(:,sp)  = abs(tx_vec(:) - syms_eq_pc(:,sp) ).^2;
    aevms(sp) = mean(evm_mat(:,sp));
    snr_mat(sp) = -10*log10(aevms (sp));
    axis tight
    hold on
    plot([1 length(evm_mat(:,sp) )], 100*[aevms(sp), aevms(sp)],'color', sec_clr,'LineWidth',4)
    hold off
    xlabel('Data Symbol Index')
    ylabel('EVM (%)');
    legend('Per-Symbol EVM','Average EVM','Location','NorthWest');
    h = text(round(.05*length(evm_mat(:,sp))), 100*aevms(sp), sprintf('Effective SINR: %.1f dB', snr_mat(sp)));
    set(h,'Color',[1 0 0])
    set(h,'EdgeColor',[1 0 0])
    set(h,'BackgroundColor',[1 1 1])
    title(sprintf('Stream from UE %d', sp));
    grid on
for sp=1:N_UE
   imagesc(1:N_OFDM_SYM, (SC_IND_DATA - N_SC/2), 100*fftshift( reshape(evm_mat(:,sp), [], N_OFDM_SYM), 1));
    hold on;
    h = line([1,N_OFDM_SYM],[0,0]);
    set(h,'Color',[1 0 0]);
    hold off;
    grid on
    xlabel('OFDM Symbol Index');
    ylabel('Subcarrier Index');
    title(sprintf('Stream from UE %d', sp));
    h = colorbar;
    set(get(h,'title'),'string','EVM (%)'); 

And we can finally print out some of the computed numbers as well:

fprintf('\n MIMO Results:\n');
fprintf('Num Bytes:   %d\n', N_UE*N_DATA_SYMS * log2(MOD_ORDER) / 8);
fprintf('Sym Errors:  %d (of %d total symbols)\n', sym_errs, N_UE * N_DATA_SYMS);
fprintf('Bit Errors:  %d (of %d total bits)\n', bit_errs, N_UE*N_DATA_SYMS * log2(MOD_ORDER));

fprintf("SNRs (dB): \n");
fprintf("%.2f\t", snr_mat);

Last updated on 12 Feb 2020 / Published on 18 Mar 2019