%% 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.

% Copyright 2019 The MathWorks, Inc.

%% Introduction
% 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.
%
% <<../5G_NRTM_EVM_ProcessingChain.png>>

%% 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);

%% 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
displayResourceGrid(tmwavegen);

%%
% 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;
    end    

    % 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);
else
    rxWaveform = txWaveform;
end

%%
% 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
     plotNonLinearCharacteristic(rapp);
    
    % Apply nonlinearity
    rxWaveform = rapp(rxWaveform);
end

%% 
% 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;
end
rxWaveform(1:nFrames*tmwaveinfo.Info.SamplesPerSubframe*10) = [];

%% Receiver
% 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 4.9.2.2-2)
    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;
        
    end
    refSqGrid = [refSqGrid refSqGridSlot];
    evSqGrid = [evSqGrid evSqGridSlot];
    
end

%% 
% Display the constellation diagram.
figure
plot(constellationSymbols,'.'); hold on
plot(constellationRef,'+')
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)
plot(0:length(evmPerSymbol)-1,evmPerSymbol,'.-',0:length(evmPerSymbol)-1,evmMaxSymbol,'.:');
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);

subplot(3,1,2)
slotNo = 0:length(evmPerSlot)-1;
plot(slotNo,evmPerSlot,'.-',slotNo,evmMaxSlot,'.:');
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.');

subplot(3,1,3)
subcarrierNo = 0:length(evmPerSubcarrier)-1;
plot(subcarrierNo,evmPerSubcarrier,'.-',subcarrierNo,evmMaxSubcarrier,'.:')
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');
    end

    % 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);
    end
end

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);
end

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)));
    else
        evmMax = NaN;
    end
end

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;    
    end

end

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);
end

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

% 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);

end

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
release(memoryLessNonlinearity); 

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

end
