EVM Measurement of 5G NR-TM Waveforms

This example measures the error vector magnitude (EVM) of an NR test model (NR-TM) waveform with phase noise and memoryless nonlinearity impairments.


The 3GPP 5G NR standard defines a set of downlink waveforms, known as test models (NR-TM), for the purpose of base station RF testing. The NR-TM for frequency range 1 (FR1) are defined in TS 38.141-1, and the NR-TM for frequency range 2 (FR2) are defined in TS 38.141-2.

This example shows how to generate an NR-TM waveform and add impairments to it. Here, we consider phase noise and memoryless nonlinearities. We then calculate the EVM of the resulting signal. The RMS and peak EVM are plotted per OFDM symbol, slot and subcarrier. The overall EVM is also calculated. The figure below shows the processing chain implemented in this example.

Simulation parameters

Each NR-TM waveform is defined by a combination of:

  • NR-TM name

  • Channel bandwidth

  • Subcarrier spacing

  • Duplexing mode

% Release 15 NR-TM test models for FR1 and FR2
% Select one of the NR-TMs among:
% "NR-FR1-TM1.1","NR-FR1-TM1.2","NR-FR1-TM2",...
% "NR-FR1-TM2a","NR-FR1-TM3.1","NR-FR1-TM3.1a",...
% "NR-FR1-TM3.2","NR-FR1-TM3.3","NR-FR2-TM1.1",...
% "NR-FR2-TM2","NR-FR2-TM3.1"
nrtm = "NR-FR1-TM3.2";

% Select the NR-TM waveform parameters
bw   = "10MHz"; % Channel bandwidth
scs  = "30kHz"; % Subcarrier spacing
dm   = "FDD";   % Duplexing mode

fprintf('TM = %s\n', nrtm);
TM = NR-FR1-TM3.2

Generate TM waveform

The channel bandwidth and subcarrier spacing combination must be a valid pair from the associated FR bandwidth configuration table. The specifications (TS 38.141-2) define FR2 NR-TM for TDD only but this example also allows FDD waveforms to be created.

Create generator object for the selected NR-TM.

tmwavegen = hNRReferenceWaveformGenerator(nrtm,bw,scs,dm);

Generate the waveform txWaveform. By default, one frame is created for FDD and two for TDD. The tmwaveinfo output contains information about the generated signal. The resourcesInfo output contains information about the resources allocated to the PDSCH and PDSCH DM-RS in the generated signal. We generate three plots:

  • Channel location in the resource grids

  • Bandwidth occupancy and guard bands plot

  • Power levels resource grids

[txWaveform,tmwaveinfo,resourcesinfo] = generateWaveform(tmwavegen);

% Visualize the associated PRB and subcarrier resource grids

The generated TM waveform may contain more than one PDSCH. The function hListTargetPDSCHs() selects the target PDSCHs to analyse. This is based on the RNTI. By default, the following RNTIs are considered for EVM calculation:

  • NR-FR1-TM2: RNTI = 2 (64QAM EVM)

  • NR-FR1-TM2a: RNTI = 2 (256QAM EVM)

  • NR-FR1-TM3.1: RNTI = 0 and 2 (64QAM EVM)

  • NR-FR1-TM3.1a: RNTI = 0 and 2 (256QAM EVM)

  • NR-FR1-TM3.2: RNTI = 2 (16QAM EVM)

  • NR-FR1-TM3.3: RNTI = 1 (QPSK EVM)

  • NR-FR2-TM2: RNTI = 2 (64QAM EVM)

  • NR-FR2-TM3.1: RNTI = 0 and 2 (64QAM EVM)

Note that, as per the specifications (TS 38.141-1, TS 38.141-2), the following TMs are not designed to perform EVM measurements: NR-FR1-TM1.1, NR-FR1-TM1.2, NR-FR2-TM1.1. However, in this example, if these TMs are generated, the EVM is measured for the following RNTIs:

  • NR-FR1-TM1.1: RNTI = 0 (QPSK EVM)

  • NR-FR1-TM1.2: RNTI = 2 (QPSK EVM)

  • NR-FR2-TM1.1: RNTI = 0 (QPSK EVM)

These values can be overridden. The function hListTargetPDSCHs() also accepts a third parameter with the set of RNTIs to consider for EVM measurement.

pdschArray is an array of the PDSCH configuration structures to analyse. Each structure contains the set of resources used for the PDSCH and the associated DM-RS and PT-RS. It also contains the set of parameters to generate the PDSCH.

[pdschArray, targetRNTIs] = hListTargetPDSCHs(tmwavegen.Config,resourcesinfo.WaveformResources);

Impairment: Phase Noise and Nonlinearity

This example considers two impairments: phase noise and memoryless nonlinearity. The impairments can be enabled or disabled by toggling the flags phaseNoiseOn and nonLinearityModelOn.

phaseNoiseOn = true;
nonLinearityModelOn = true;

Normalize the waveform to fit the dynamic range of the nonlinearity.

txWaveform = txWaveform/max(abs(txWaveform));

The waveform consists of one frame for FDD and two for TDD. Repeat the signal twice. We will remove the first half of the resulting waveform to avoid the transient introduced by the phase noise model.

txWaveform = repmat(txWaveform,2,1);

Introduce phase noise distortion. The generated figure shows the phase noise characteristic. The carrier frequency considered depends on the frequency range. We use values of 4 GHz and 28 GHz for FR1 and FR2 respectively. The phase noise characteristic is generated with the pole/zero model described in R1-163984, "Discussion on phase noise modeling".

if phaseNoiseOn
    sr = tmwaveinfo.Info.SamplingRate;

    % Carrier frequency
    if tmwavegen.Config.FrequencyRange == "FR1" % carrier frequency for FR1
        fc = 4e9;
    else % carrier frequency for FR2
        fc = 28e9;

    % Phase noise level
    foffsetLog = (3:0.1:log10(sr/2)); % model offset from 1e3Hz to sr/2
    foffset = 10.^foffsetLog;    % linear freq offset
    PN_dBc_Hz = PNmodelPoleZero(foffset,fc);
    figure; semilogx(foffset,PN_dBc_Hz); xlabel('Frequency offset (Hz)'); ylabel('dBc/Hz'); title('Phase noise model'); grid on

    % Apply phase noise to waveform
    pnoise = comm.PhaseNoise('FrequencyOffset',foffset,'Level',PN_dBc_Hz,'SampleRate',sr);
    rxWaveform = pnoise(txWaveform);
    rxWaveform = txWaveform;

Introduce non-linear distortion. We use the Rapp model in this example. The figure shows the nonlinearity introduced. The parameters for the Rapp model are chosen to match the characteristic of the memoryless model from TR 38.803 (Memoryless polynomial model - Annex A.1).

if nonLinearityModelOn
    rapp = comm.MemorylessNonlinearity('Method','Rapp model');
    rapp.Smoothness = 1.55;
    rapp.OutputSaturationLevel = 1;

    % Plot non-linear characteristic

    % Apply nonlinearity
    rxWaveform = rapp(rxWaveform);

The signal was previously repeated twice. We now remove the first half of this signal. This avoids any transient introduced by the impairment models.

if dm == "FDD"
    nFrames = 1;
else % TDD
    nFrames = 2;
rxWaveform(1:nFrames*tmwaveinfo.Info.SamplesPerSubframe*10) = [];


The receiver in this example performs the following steps:

  • Synchronization using the DM-RS over one frame for FDD (two frames for TDD)

  • OFDM demodulation of the received waveform

  • Channel estimation

  • Equalisation

Get waveform parameters for synchronization and OFDM demodulation

gnb.NRB = tmwavegen.Config.Carriers.NRB;
gnb.CyclicPrefix = tmwavegen.Config.BWP.CyclicPrefix;
gnb.SubcarrierSpacing = tmwavegen.Config.Carriers.SubcarrierSpacing;

Generate a reference grid spanning 10 msec (one frame). This grid contains only the DM-RS and is used for synchronization.

refGrid = referenceGrid(tmwaveinfo,pdschArray);
NSlot = 0;
[offset,mag] = nrTimingEstimate(rxWaveform,gnb.NRB,gnb.SubcarrierSpacing,NSlot,refGrid);
tmWaveformSync = rxWaveform(1+offset:end,:);

OFDM demodulation.

rxGrid = hOFDMDemodulate(gnb, tmWaveformSync);

For all slots, get the PDSCH and DM-RS resources (resource element locations). This information is present in pdschArray. Then perform channel estimation and equalization. The resulting data is stored to calculate the EVM.

% Get total number of slots and the number of subcarriers and symbols per
% slot
symbolsPerSlot = tmwaveinfo.Info.SymbolsPerSlot;
totalNrSlots = floor(size(rxGrid,2)/symbolsPerSlot);
noSubcarriers = size(rxGrid,1);

% Declare storage variables
constellationSymbols = [];  % equalized symbols for constellation plot
constellationRef = [];      % reference symbols for constellation plot
refSqGrid = []; % Grid with magnitude square of reference symbols
evSqGrid = [];  % Grid with magnitude square of error vector (symbols - reference)

for NSlot=0:totalNrSlots-1
    % Extract grid for current slot
    currentGrid = rxGrid(:,NSlot*symbolsPerSlot+(1:symbolsPerSlot),:);

    % Retrieve resources
    [pdschIndices,refSymbols,dmrsIndices,dmrsSymbols] = hSlotResources(pdschArray,NSlot);

    % Do not include first two slot symbols for PDSCH EVM (used for control
    % TS 38.141-1 table
    idx = pdschIndices <= 2*noSubcarriers;
    pdschIndices(idx) = [];
    refSymbols(idx) = [];

    % Channel estimation
    [estChannelGrid,~] = nrChannelEstimate(currentGrid,dmrsIndices,dmrsSymbols);

    % Get PDSCH resource elements from the received grid
    [pdschRx,pdschHest] = nrExtractResources(pdschIndices,currentGrid,estChannelGrid);

    % Equalization: set noiseEst to 0 for zero-forcing equalization
    noiseEst = 0;
    [pdschEq,csi] = nrEqualizeMMSE(pdschRx,pdschHest,noiseEst);

    refSqGridSlot = zeros(size(currentGrid)); % slot grid magnitude square of reference symbols
    evSqGridSlot = zeros(size(currentGrid));  % slot grid magnitude square of error vector
    if ~isempty(refSymbols) % for slots with data

        % Error vector magnitude square
        evSq = abs(refSymbols-pdschEq).^2;
        % Reference symbols magnitude square
        refSymSq = abs(refSymbols).^2;

        % Store constellation symbols, reference symbols and square error vector
        constellationSymbols = [constellationSymbols; pdschEq];
        constellationRef = [constellationRef; refSymbols];
        evSqGridSlot(pdschIndices) = evSq;
        refSqGridSlot(pdschIndices) = refSymSq;

    refSqGrid = [refSqGrid refSqGridSlot];
    evSqGrid = [evSqGrid evSqGridSlot];


Display the constellation diagram.

plot(constellationSymbols,'.'); hold on
title('Equalized symbols constellation')
grid on; xlabel('In-Phase'); ylabel('Quadrature');

EVM Calculation

Calculate the following EVM (RMS and maximum) values:

  • EVM per OFDM symbol

  • EVM per slot

  • EVM per subcarrier

  • Total EVM for the whole signal

The EVM values shown are normalized by the power of the reference symbols. This power is calculated for the measurement interval considered. For example, when calculating the EVM per slot, all the reference symbols in that slot are used to calculate the power used in the normalization. In the case of the EVM per subcarrier, the measurement interval is one frame for FDD and two frames for TDD. The total EVM is measured over one frame for FDD and two frames for TDD.

Display the EVM per OFDM symbol, slot and subcarrier.

% EVM per symbol
[evmPerSymbol, evmMaxSymbol] = evm(evSqGrid,refSqGrid);

figure; subplot(3,1,1)
title('EVM vs OFDM symbol')
grid on; xlabel('Symbol number'); ylabel('EVM (%)'); legend('rms EVM','peak EVM','Location','bestoutside')

% EVM per slot
% reshape grids to one column per slot
evSqGridSlot = reshape(evSqGrid,size(evSqGrid).*symbolsPerSlot.^[1 -1]);
refSqGridSlot = reshape(refSqGrid,size(refSqGrid).*symbolsPerSlot.^[1 -1]);

[evmPerSlot, evmMaxSlot] = evm(evSqGridSlot,refSqGridSlot);

slotNo = 0:length(evmPerSlot)-1;
title('EVM vs slot')
grid on; xlabel('Slot number'); ylabel('EVM (%)'); legend('rms EVM','peak EVM','Location','bestoutside')

% EVM per subcarrier
[evmPerSubcarrier, evmMaxSubcarrier] = evm(evSqGrid.',refSqGrid.');

subcarrierNo = 0:length(evmPerSubcarrier)-1;
title('EVM vs subcarrier')
grid on; xlabel('Subcarrier number'); ylabel('EVM (%)'); legend('rms EVM','max EVM','Location','bestoutside')

Calculate overall EVM.

[evmRMS, evmMax] = evm(evSqGrid(:),refSqGrid(:));
disp("PDSCH RNTIs: " + mat2str(targetRNTIs))
disp("RMS EVM = " + string(evmRMS) + "%");
disp("Max EVM = " + string(evmMax) + "%");

Local functions

function [evmRMS, evmMax] = evm(evSqGrid,refSqGrid)
% Calculates the evm per column of the input matrix
% inputs (both of the same size):
%   - evSqGrid = matrix of error vector magnitude square |refSymb-rxSymb|.^2
%   - refSqGrid = matrix of ref symbols magnitude square |refSymb|.^2

    if size(evSqGrid)~=size(refSqGrid)
        error('Input matrices must have the same size');

    % RMS EVM (%)
    evmRMS = rmsEVM(sum(evSqGrid,1),sum(refSqGrid,1));

    % max EVM (%)
    evmMax = zeros(1,size(refSqGrid,2));
    % find non zero REs in reference grid
    nonZeroREs = (refSqGrid~=0);
    for ncol = 1:length(evmMax)
        % remove reference symbols with value zero
        nzEvSqREsSlot = evSqGrid(nonZeroREs(:,ncol),ncol);
        nzRefSqREsSlot = refSqGrid(nonZeroREs(:,ncol),ncol);
        evmMax(ncol) = maxEVM(nzEvSqREsSlot,nzRefSqREsSlot);

function evmRMS = rmsEVM(evSq,refSq)
% Calcuates RMS EVM. Inputs
%   - evSq = matrix of error vector magnitude square
%   - refSq = matrix of ref symbols magnitude square
% They have the same size or refSq can be a scalar.
    evmRMS = 100*sqrt(evSq./refSq);

function evmMax = maxEVM(evSq,refSq)
% inputs (both of the same size):
%   - evSq = array of error vector magnitude square
%   - refSq = array of ref symbols magnitude square
    if ~isempty(evSq) && ~isempty(refSq)
        evmMax = max(100*sqrt(evSq./mean(refSq)));
        evmMax = NaN;

function refGrid = referenceGrid(tmwaveinfo,pdschArray)
% Create a reference grid for synchronization. This should span a whole
% frame, starting in slot 0. It contains the DM-RSs specified in
% pdschArray.

    L = tmwaveinfo.Info.SymbolsPerSubframe*10; % symbols in a frame
    refGrid = zeros(tmwaveinfo.Info.NSubcarriers,L); % empty grid
    rbsPerSlot = tmwaveinfo.Info.NSubcarriers*tmwaveinfo.Info.SymbolsPerSlot;
    % Polpulate the DM-RSs in the reference grid for all slots
    for NSlot=0:(10*tmwaveinfo.Info.SlotsPerSubframe)-1 % 1 frame worth of subframes
        [~,~,dmrsIndices,dmrsSymbols] = hSlotResources(pdschArray,NSlot);
        refGrid(dmrsIndices+NSlot*rbsPerSlot) = dmrsSymbols;


function PN_dBC_Hz = PNmodelPoleZero(f,fc)
% Generate the phase noise characteristic in dBc/Hz for the frequency
% offset values specified by vector f for the carrier frequency fc.
% The model used here is the multi-pole/zero model as proposed in:
% Yinan Qi et al. "On the Phase Tracking Reference Signal (PT-RS) Design
% for 5G New Radio (NR)". Vehicular Technology Conference (VTC-Fall), Aug
% 2018.

% Pole/zeros and PSD0 as specified in the paper mentioned above and in 3GPP
% R1-163984, "Discussion on phase noise modeling". This corresponds to
% "parameter set A" for a base frequency of 30GHz.
fcBase = 30e9;
fz= [1.8 2.2 40]*1e6;
fp = [0.1 0.2 8]*1e6;
PSD0 = -79.4;

% Compute numerator
num = ones(size(f));
for ii=1:numel(fz)
    num = num.* ( 1 +(f./fz(ii)).^2);

% Compute denominator
den = ones(size(f));
for ii=1:numel(fz)
    den = den.* ( 1 +(f./fp(ii)).^2);

% Compute phase noise and apply a shift for carrier frequencies different
% from the base frequency.
PN_dBC_Hz = 10*log10(num./den) + PSD0 + 20*log10(fc/fcBase);


function plotNonLinearCharacteristic(memoryLessNonlinearity)
% Plot the non linear characteristic of the PA impairment represented by
% memoryLessNonlinearity system object. This input parameter is a
% comm.MemorylessNonlinearity system object.

% input samples
x = complex((1/sqrt(2))*(-1+2*rand(1000,1)),(1/sqrt(2))*(-1+2*rand(1000,1)));

% Non linearity
yRapp = memoryLessNonlinearity(x);
% release object so we can feed it different number of samples

% Plot characteristic
figure; plot(10*log10(abs(x).^2),10*log10(abs(x).^2)); hold on; grid on
xlabel('Input power (dBW)'); ylabel('Output power (dBW)'); title('nonlinearity impairment')
legend('Linear characteristic', 'Rapp nonlinearity','Location','Northwest');


Related Topics