Measure Impact of Sub-THz Hardware Impairments on 6G Waveforms

Since R2024a

This example shows how to explore the impact of hardware impairments at sub-THz frequencies on a candidate 6G waveform. The hardware impairments are: phase noise, power amplifier (PA) nonlinearity, and a filter to limit spectral emissions outside of the channel bandwidth. The examples measures the adjacent channel power ratio (ACPR) and error vector magnitude (EVM) of the impaired waveform to demonstrate the impact of impairments.

Set System Parameters

This example models a 6G-like OFDM waveform with system parameters defined in Hexa-X Deliverable D2.3 Table 5-7 [1] to investigate the impact of hardware impairments on candidate 6G waveforms.

carrierFrequency = 140; % GHz
bandwidth = 2160; % MHz - transmission bandwidth and guardbands
subcarrierSpacing = 3840; % kHz

Set the bandwidth occupancy. This is the ratio between the transmission bandwidth and channel bandwidth. A lower bandwidth occupancy increases the size of guardbands, reducing spectral emissions at the cost of reduced spectral efficiency.

bandwidthOccupancy = 0.8; % Ratio of transmission bandwidth to channel bandwidth

Create a carrier configuration. Calculate the transmission bandwidth and required number of resource blocks (NSizeGrid) given the bandwidth occupancy.

carrier = pre6GCarrierConfig;
carrier.SubcarrierSpacing = subcarrierSpacing;
channelBandwidth = bandwidth*1e6; % Bandwidth of RBs and guard carriers in Hz
carrier.NSizeGrid = floor((channelBandwidth/(carrier.SubcarrierSpacing*1e3)*bandwidthOccupancy)/12); % 12 subcarriers per RB
transmissionBandwidth = carrier.SubcarrierSpacing*1e3*carrier.NSizeGrid*12; % Bandwidth of RBs in Hz

This example creates a time division duplex (TDD) waveform consisting of a pattern of downlink slots followed by uplink slots with a specified period. Set the length of the waveform and TDD parameters.

numSubframes = 0.5; % Number of subframes to simulate
tddConfig.TDDPeriod = 10; % TDD period in slots
tddConfig.NumDownlinkSlots = 9; % Number of slots in TDD period containing PDSCH

numSlots = numSubframes*carrier.SlotsPerSubframe;
tddConfig.SlotAllocation = 0:tddConfig.NumDownlinkSlots-1; % Assume downlink slots at start of TDD period
disp(join(["Simulating" num2str(numSlots) "slots"]))

Simulating 128 slots

visualizeTDDAllocation(tddConfig);

Configure a physical downlink shared channel (PDSCH) for full resource block (RB) allocation and enable the phase tracking reference signal (PT-RS).

pdsch = pre6GPDSCHConfig;
pdsch.PRBSet = 0:(carrier.NSizeGrid-1);
pdsch.Modulation = "16QAM";
pdsch.EnablePTRS = true;

Generate Waveform

Create a resource grid containing a PDSCH transmission, demodulation reference signal (DM-RS), and PT-RS in appropriate TDD slots.

waveGrid = []; % Store resource grid to transmit
for nSlot = 0:(numSlots-1)
    carrier.NSlot = nSlot;
    slotGrid = hpre6GResourceGrid(carrier,pdsch.NumLayers);

    % Generate content if slot allocated to PDSCH
    if isSlotActive(carrier.NSlot,tddConfig)
        [ind,indInfo] = hpre6GPDSCHIndices(carrier,pdsch);
        % Generate and map PDSCH for a random codeword
        cw = randi([0 1],indInfo.G,1);
        sym = hpre6GPDSCH(carrier,pdsch,cw);
        slotGrid(ind) = sym;

        % Generate and map DMRS
        dmrsInd = hpre6GPDSCHDMRSIndices(carrier,pdsch);
        dmrsSym = hpre6GPDSCHDMRS(carrier,pdsch);
        slotGrid(dmrsInd) = dmrsSym;

        % Generate and map PT-RS
        ptrsInd = hpre6GPDSCHPTRSIndices(carrier,pdsch);
        ptrsSym = hpre6GPDSCHPTRS(carrier,pdsch);
        slotGrid(ptrsInd) = ptrsSym;
    end
    waveGrid = [waveGrid slotGrid];
end

OFDM-modulate the resource grid with an oversampling factor. Oversample the waveform to measure spectral emissions outside of the channel bandwidth. The example sets the FFT size such that it spans the specified oversampled channel bandwidth and is an integer multiple of 128 for a normal cyclic prefix length.

% OFDM-modulate with specified oversampling factor
oversamplingFactor = 3; % Multiple of the channel bandwidth
carrier.NSlot = 0; % OFDM-modulate from slot 0
nFFTWaveform = ceil((channelBandwidth/(carrier.SubcarrierSpacing*1e3))*oversamplingFactor/128)*128;
[txWaveform,ofdmInfo] = hpre6GOFDMModulate(carrier,waveGrid,Nfft=nFFTWaveform);

% Normalize waveform to maximum amplitude for PA modeling
txWaveform = txWaveform/max(abs(txWaveform),[],"all");

% Prepend slot to the waveform to allow to allow for filter delay
waveformLength = height(txWaveform);
slotLength = waveformLength/numSlots;
txWaveform = [txWaveform(end-slotLength+1:end,:); txWaveform];

Add Impairments

In this section, you model phase noise and power amplifier (PA) nonlinearity impairments, and apply a filter to limit spectral emissions outside of the channel bandwidth.

Phase Noise

Introduce phase noise distortion. This example uses a phase noise model defined in Hexa-X Deliverable D2.3 Section 4.2.1.1, scaled to 140 GHz as specified in Hexa-X Deliverable D2.3 Table 5-7. This Hexa-X model is based on measurements of 20 GHz wideband RF synthesizer TI LMX2596. Set the minimum frequency offset the phase noise model uses to compute the phase noise spectrum mask. By default, this example uses a 1 MHz minimum frequency offset, therefore the phase noise level below 1 MHz does not match the Hexa-X model. Simulating with a lower minimum frequency offset will more accurately model the low frequency phase noise of the Hexa-X model but result in a longer time to design the phase noise filter.

enablePN = true; % Enable phase noise model
if enablePN
    minimumFrequencyOffset = 1e6; % Hz
    mdl = "HexaX Model 2";
    phaseNoise = hpre6GPhaseNoise(carrierFrequency*1e9,ofdmInfo.SampleRate, ...
        MinFrequencyOffset=minimumFrequencyOffset,Model=mdl, ...
        RandomStream="mt19937ar with seed");
    visualize(phaseNoise); % Show phase noise PSD
    rxWaveform = phaseNoise(txWaveform); % Apply phase noise model
else
    rxWaveform = txWaveform; %#ok<UNRCH>
end

Low-Pass Filter

Filter the baseband waveform to limit emissions outside the channel bandwidth. If the current passband and stopband frequencies, PassbandFrequency and StopbandFrequency, result in high EVM values, use a wider filter by increasing PassbandFrequency and StopbandFrequency. To use a narrower filter, reduce PassbandFrequency and StopbandFrequency. You can also modify the PassbandRipple and the StopbandAttenuation.

enableLPF = true; % Enable low-pass filter
if enableLPF
    % Create low-pass filter object
    LPF = dsp.LowpassFilter;
    LPF.SampleRate = ofdmInfo.SampleRate;
    LPF.FilterType = "IIR";
    LPF.PassbandFrequency = (transmissionBandwidth + 24*carrier.SubcarrierSpacing*1e3)/2;
    LPF.StopbandFrequency = channelBandwidth/2; % TransmissionBandwidth and guards
    LPF.PassbandRipple = 0.2;
    LPF.StopbandAttenuation = 40;
    figure;
    freqz(LPF); % Plot the response of the low-pass filter

    % Filter the waveform
    rxWaveform = LPF(rxWaveform);
    release(LPF);
end

Power Amplifier

Apply a PA impairment model. This example applies a memoryless GaN PA model as described in Hexa-X Deliverable D2.3 Table 5-7 and Hexa-X Deliverable D2.2 Section 3.6.6.1 [2]. This model assumes the ACPR characteristic with respect to PA output power shift according to the expected downscaling of maximum PA output power when increasing frequency.

Set the model backoff to reduce the amplitude of the input signal and reduce distortion.

enablePA = true; % Enable power amplifier model
if enablePA
    backoff = 6; % In dB
    rxWaveform = db2mag(-backoff)*rxWaveform; % Apply PA backoff
    visualizeAMAMCharacteristic(@paMemorylessGaN,"GaN");
    
    rxWaveform = paMemorylessGaN(rxWaveform);    
end

Measure ACPR

Measure the ACPR to study spectral regrowth caused by the nonlinear PA model.

measureACPR = true;
if measureACPR
    % Calculate the number of adjacent channels which can be measured given
    % the sample rate
    numAdjacentChannels = floor((ofdmInfo.SampleRate/channelBandwidth-1)/2);
    if numAdjacentChannels>0        
        sa = spectrumAnalyzer;
        sa.SampleRate = ofdmInfo.SampleRate;
        sa.ChannelMeasurements.Type = "acpr";
        sa.ChannelMeasurements.Enabled = true;
        sa.ChannelMeasurements.Span = transmissionBandwidth;
        sa.ChannelMeasurements.NumOffsets = numAdjacentChannels;
        sa.ChannelMeasurements.AdjacentBW = transmissionBandwidth;
        sa.ChannelMeasurements.ACPROffsets = (1:numAdjacentChannels)*channelBandwidth;
        sa(rxWaveform);
    else
        warning('Sample rate too low to measure ACPR, increase oversamplingFactor')
    end
end

Measure EVM

Perform the following steps to measure the EVM:

1. Synchronize the received waveform

2. OFDM-demodulate the received waveform

3. Perform channel estimation

4. Equalize the PDSCH symbols

5. Estimate and compensate the common phase error (CPE)

6. Compute PDSCH EVM

Show the equalized PDSCH symbols for each slot.

% Ignore transients in first part of waveform
rxWaveform(1:slotLength,:) = [];

% Timing synchronization
refSlotGrid = waveGrid(:,1:carrier.SymbolsPerSlot,:);
offset = hpre6GTimingEstimate(carrier,rxWaveform,refSlotGrid,Nfft=ofdmInfo.Nfft);
rxWaveform = rxWaveform(1+offset:end,:);

% Demodulate and measure EVM of each slot
rxWaveGrid = hpre6GOFDMDemodulate(carrier,rxWaveform,Nfft=ofdmInfo.Nfft);

% Setup constellation diagram
constDiagram = comm.ConstellationDiagram;
constDiagram.ShowReferenceConstellation = true;

L = carrier.SymbolsPerSlot;
numRxSlots = floor(size(rxWaveGrid,2)/L);
if numRxSlots<1
    error('Not enough data to measure EVM, increase numSubframes')
end
evGrid = []; % Store error vector for each resource element
evm = []; % Store EVM measurement for each slot containing PDSCH
for nSlot = 0:(numRxSlots-1)
    carrier.NSlot = nSlot;
    evSlotGrid = NaN*hpre6GResourceGrid(carrier,pdsch.NumLayers); % Store error vector
    if isSlotActive(nSlot,tddConfig)
        % Extract slot from grid
        rxSlotGrid = rxWaveGrid(:,nSlot*L+(1:L),:);

        % Estimate channel
        pdschDMRSInd = hpre6GPDSCHDMRSIndices(carrier,pdsch);
        pdschDMRSSymbols = hpre6GPDSCHDMRS(carrier,pdsch);
        [hest,nVar] = hpre6GChannelEstimate(carrier,rxSlotGrid,pdschDMRSInd,pdschDMRSSymbols,CDMLengths=pdsch.DMRS.CDMLengths);

        % Equalize PDSCH
        [pdschInd,pdschIndInfo] = hpre6GPDSCHIndices(carrier,pdsch);
        [pdschRxSym,pdschHest] = nrExtractResources(pdschInd,rxSlotGrid,hest);
        pdschEqSym = nrEqualizeMMSE(pdschRxSym,pdschHest,nVar);

        % Correct common phase error
        pdschEqSym = correctCPE(carrier,pdsch,pdschInd,pdschIndInfo,rxSlotGrid,hest,nVar,pdschEqSym);

        % Demap, decode, and re-modulate to obtain layer demapped equalized
        % symbols (rxSym) and reference symbols (refSym) for EVM
        % measurement.
        [cw,rxSym] = hpre6GPDSCHDecode(carrier,pdsch,pdschEqSym,nVar);
        expectedCW = cw{1}<0; % LLR<0 is bit 1, otherwise bit 0
        expectedPDSCHSym = hpre6GPDSCH(carrier,pdsch,expectedCW);
        [~,refSym] = hpre6GPDSCHDecode(carrier,pdsch,expectedPDSCHSym,0);

        % Show equalized PDSCH constellation
        constDiagram.ReferenceConstellation = unique(refSym{1});
        constDiagram.Title = join(["Equalized PDSCH slot" num2str(nSlot)]);
        constDiagram(rxSym{1});

        % Calculate error vector and EVM
        slotEVM = measureEVM(refSym{1},rxSym{1});
        fprintf("Slot %d PDSCH EVM, RMS: %0.3f%% Peak: %0.3f%%\n",nSlot,slotEVM.RMS*100,slotEVM.Peak*100);
        evm = [evm slotEVM];

        % Store error vector in grid
        evSlotGrid = NaN(size(rxSlotGrid));
        evSlotGrid(pdschInd) = slotEVM.EV;
    end
    evGrid = [evGrid evSlotGrid]; %#ok<*AGROW>
end

Slot 0 PDSCH EVM, RMS: 3.755% Peak: 16.844%
Slot 1 PDSCH EVM, RMS: 3.100% Peak: 15.873%
Slot 2 PDSCH EVM, RMS: 3.036% Peak: 12.505%
Slot 3 PDSCH EVM, RMS: 2.989% Peak: 11.587%
Slot 4 PDSCH EVM, RMS: 2.901% Peak: 11.346%
Slot 5 PDSCH EVM, RMS: 2.573% Peak: 8.262%
Slot 6 PDSCH EVM, RMS: 2.675% Peak: 8.874%
Slot 7 PDSCH EVM, RMS: 2.788% Peak: 10.717%
Slot 8 PDSCH EVM, RMS: 2.911% Peak: 13.210%
Slot 10 PDSCH EVM, RMS: 3.032% Peak: 8.778%
Slot 11 PDSCH EVM, RMS: 2.729% Peak: 8.867%
Slot 12 PDSCH EVM, RMS: 3.036% Peak: 11.233%
Slot 13 PDSCH EVM, RMS: 3.365% Peak: 11.645%
Slot 14 PDSCH EVM, RMS: 2.763% Peak: 10.647%
Slot 15 PDSCH EVM, RMS: 2.611% Peak: 8.835%
Slot 16 PDSCH EVM, RMS: 3.086% Peak: 12.029%
Slot 17 PDSCH EVM, RMS: 3.428% Peak: 12.731%
Slot 18 PDSCH EVM, RMS: 2.794% Peak: 8.649%
Slot 20 PDSCH EVM, RMS: 3.044% Peak: 13.066%
Slot 21 PDSCH EVM, RMS: 2.572% Peak: 8.695%
Slot 22 PDSCH EVM, RMS: 2.999% Peak: 11.147%
Slot 23 PDSCH EVM, RMS: 3.131% Peak: 13.562%
Slot 24 PDSCH EVM, RMS: 3.420% Peak: 12.666%
Slot 25 PDSCH EVM, RMS: 4.019% Peak: 16.212%
Slot 26 PDSCH EVM, RMS: 2.689% Peak: 10.144%
Slot 27 PDSCH EVM, RMS: 2.467% Peak: 7.728%
Slot 28 PDSCH EVM, RMS: 2.804% Peak: 9.291%
Slot 30 PDSCH EVM, RMS: 2.796% Peak: 13.387%
Slot 31 PDSCH EVM, RMS: 3.154% Peak: 11.011%
Slot 32 PDSCH EVM, RMS: 2.747% Peak: 9.354%
Slot 33 PDSCH EVM, RMS: 2.636% Peak: 9.154%
Slot 34 PDSCH EVM, RMS: 2.934% Peak: 11.586%
Slot 35 PDSCH EVM, RMS: 2.540% Peak: 8.860%
Slot 36 PDSCH EVM, RMS: 3.053% Peak: 9.888%
Slot 37 PDSCH EVM, RMS: 3.156% Peak: 12.925%
Slot 38 PDSCH EVM, RMS: 3.190% Peak: 15.825%
Slot 40 PDSCH EVM, RMS: 3.005% Peak: 10.599%
Slot 41 PDSCH EVM, RMS: 2.528% Peak: 8.828%
Slot 42 PDSCH EVM, RMS: 2.824% Peak: 10.614%
Slot 43 PDSCH EVM, RMS: 2.726% Peak: 10.576%
Slot 44 PDSCH EVM, RMS: 2.735% Peak: 10.777%
Slot 45 PDSCH EVM, RMS: 2.748% Peak: 9.195%
Slot 46 PDSCH EVM, RMS: 2.985% Peak: 9.066%
Slot 47 PDSCH EVM, RMS: 2.539% Peak: 8.803%
Slot 48 PDSCH EVM, RMS: 3.371% Peak: 13.485%
Slot 50 PDSCH EVM, RMS: 2.805% Peak: 8.677%
Slot 51 PDSCH EVM, RMS: 2.428% Peak: 9.136%
Slot 52 PDSCH EVM, RMS: 3.684% Peak: 17.042%
Slot 53 PDSCH EVM, RMS: 2.802% Peak: 10.523%
Slot 54 PDSCH EVM, RMS: 2.804% Peak: 9.368%
Slot 55 PDSCH EVM, RMS: 2.788% Peak: 10.153%
Slot 56 PDSCH EVM, RMS: 3.136% Peak: 12.831%
Slot 57 PDSCH EVM, RMS: 2.739% Peak: 8.204%
Slot 58 PDSCH EVM, RMS: 2.984% Peak: 11.302%
Slot 60 PDSCH EVM, RMS: 3.007% Peak: 10.761%
Slot 61 PDSCH EVM, RMS: 2.622% Peak: 8.684%
Slot 62 PDSCH EVM, RMS: 3.276% Peak: 11.094%
Slot 63 PDSCH EVM, RMS: 3.290% Peak: 13.024%
Slot 64 PDSCH EVM, RMS: 2.962% Peak: 10.227%
Slot 65 PDSCH EVM, RMS: 2.817% Peak: 9.630%
Slot 66 PDSCH EVM, RMS: 3.134% Peak: 11.550%
Slot 67 PDSCH EVM, RMS: 3.220% Peak: 12.505%
Slot 68 PDSCH EVM, RMS: 3.263% Peak: 13.111%
Slot 70 PDSCH EVM, RMS: 2.905% Peak: 10.428%
Slot 71 PDSCH EVM, RMS: 2.768% Peak: 10.422%
Slot 72 PDSCH EVM, RMS: 2.759% Peak: 9.903%
Slot 73 PDSCH EVM, RMS: 3.405% Peak: 15.982%
Slot 74 PDSCH EVM, RMS: 2.293% Peak: 6.821%
Slot 75 PDSCH EVM, RMS: 2.835% Peak: 8.695%
Slot 76 PDSCH EVM, RMS: 3.074% Peak: 10.207%
Slot 77 PDSCH EVM, RMS: 2.463% Peak: 8.090%
Slot 78 PDSCH EVM, RMS: 2.713% Peak: 8.993%
Slot 80 PDSCH EVM, RMS: 3.273% Peak: 12.462%
Slot 81 PDSCH EVM, RMS: 2.916% Peak: 10.495%
Slot 82 PDSCH EVM, RMS: 3.408% Peak: 14.847%
Slot 83 PDSCH EVM, RMS: 3.049% Peak: 10.221%
Slot 84 PDSCH EVM, RMS: 3.153% Peak: 11.632%
Slot 85 PDSCH EVM, RMS: 3.357% Peak: 11.146%
Slot 86 PDSCH EVM, RMS: 3.099% Peak: 11.029%
Slot 87 PDSCH EVM, RMS: 3.073% Peak: 12.753%
Slot 88 PDSCH EVM, RMS: 3.166% Peak: 13.436%
Slot 90 PDSCH EVM, RMS: 2.896% Peak: 11.020%
Slot 91 PDSCH EVM, RMS: 2.797% Peak: 9.538%
Slot 92 PDSCH EVM, RMS: 2.919% Peak: 13.055%
Slot 93 PDSCH EVM, RMS: 2.505% Peak: 9.229%
Slot 94 PDSCH EVM, RMS: 2.592% Peak: 8.633%
Slot 95 PDSCH EVM, RMS: 2.651% Peak: 9.568%
Slot 96 PDSCH EVM, RMS: 3.260% Peak: 10.910%
Slot 97 PDSCH EVM, RMS: 2.999% Peak: 9.820%
Slot 98 PDSCH EVM, RMS: 2.925% Peak: 8.608%
Slot 100 PDSCH EVM, RMS: 3.065% Peak: 11.772%
Slot 101 PDSCH EVM, RMS: 2.958% Peak: 9.612%
Slot 102 PDSCH EVM, RMS: 2.985% Peak: 12.338%
Slot 103 PDSCH EVM, RMS: 2.939% Peak: 14.105%
Slot 104 PDSCH EVM, RMS: 2.889% Peak: 10.277%
Slot 105 PDSCH EVM, RMS: 2.746% Peak: 11.697%
Slot 106 PDSCH EVM, RMS: 2.597% Peak: 8.738%
Slot 107 PDSCH EVM, RMS: 3.008% Peak: 13.573%
Slot 108 PDSCH EVM, RMS: 3.126% Peak: 10.318%
Slot 110 PDSCH EVM, RMS: 3.155% Peak: 12.611%
Slot 111 PDSCH EVM, RMS: 3.385% Peak: 13.498%
Slot 112 PDSCH EVM, RMS: 2.775% Peak: 9.324%
Slot 113 PDSCH EVM, RMS: 3.258% Peak: 12.520%
Slot 114 PDSCH EVM, RMS: 2.816% Peak: 13.212%
Slot 115 PDSCH EVM, RMS: 3.293% Peak: 12.403%
Slot 116 PDSCH EVM, RMS: 2.683% Peak: 9.193%
Slot 117 PDSCH EVM, RMS: 3.166% Peak: 13.103%
Slot 118 PDSCH EVM, RMS: 3.128% Peak: 11.523%
Slot 120 PDSCH EVM, RMS: 3.040% Peak: 12.156%
Slot 121 PDSCH EVM, RMS: 3.374% Peak: 13.415%
Slot 122 PDSCH EVM, RMS: 3.069% Peak: 13.670%
Slot 123 PDSCH EVM, RMS: 2.866% Peak: 10.491%
Slot 124 PDSCH EVM, RMS: 2.876% Peak: 11.730%
Slot 125 PDSCH EVM, RMS: 2.654% Peak: 9.544%
Slot 126 PDSCH EVM, RMS: 2.971% Peak: 10.974%

Measure RMS and peak EVM for whole waveform. Elements of evGrid contain NaN when no PDSCH is present, therefore, ignore these for the measurement (omitmissing).

rmsEVM = sqrt(mean(abs(evGrid).^2,"all","omitmissing"))*100;
peakEVM = max(abs(evGrid),[],"all","omitmissing")*100;

fprintf("Overall PDSCH EVM, RMS: %0.3f%% Peak: %0.3f%%\n",rmsEVM,peakEVM);

Overall PDSCH EVM, RMS: 2.972% Peak: 17.042%

Plot EVM per subcarrier, averaged over OFDM symbols.

evmSubcarrierRMS = plotEVMPerSubcarrier(evGrid);

Plot EVM per symbol, averaged over subcarriers.

evmSymbolRMS = plotEVMPerSymbol(evGrid);

References

[1] Hexa-X Deliverable D2.3 - Radio models and enabling techniques towards ultra-high data rate links and capacity in 6G

[2] Hexa-X Deliverable D2.2 - Initial radio models and analysis towards ultra-high data rate links in 6G

Local Functions

function y = paMemorylessGaN(x)
    % Model 28 GHz GaN power amplifier as defined in R4-165901 Section 2.1.4

    ak = [-0.334697-0.942326i; 0.89015-0.72633i; -2.58056+4.81215i; 4.81548-9.54837i; -4.41452+8.63164i; 1.54271-2.94034i];
    y = x.*abs(x).^(2*(0:5))*ak;
end

function visualizeAMAMCharacteristic(paModel,modelName)
    % Plot the nonlinear characteristic of the power amplifier impairment

    x = logspace(-1.5,0).'; % Input samples with normalized input power
    y = paModel(x); % Nonlinearity

    figure; 
    plot(pow2db(abs(x).^2),pow2db(abs(y).^2),'-');
    hold on;
    plot(pow2db(abs(x).^2),pow2db(abs(x).^2),'-.'); 
    grid on
    xlabel("Instantaneous normalized Input Power (dBW)"); 
    ylabel("Instanteneous Output Power (dBW)"); 
    title(join(["AM/AM," modelName]))
    legend(join([modelName "characteristic"]),"Linear characteristic",Location="Northwest");
end

function visualizeTDDAllocation(tddConfig)
    % Show the TDD allocation pattern

    x = [0 1 1 0];
    y = [0 0 1 1];
    xTxtPos = (x(2)-x(1))/2+x(1);
    yTxtPos = (y(3)-y(1))/2+y(1)-1; % Position text below patch blocks
    
    hp = [];
    t = [];
    f = figure;
    for i = 0:tddConfig.TDDPeriod-1
        if isSlotActive(i,tddConfig)
            patchColorIdx = 1; % DL slot
            txt = "D";
        else
            patchColorIdx = 2; % UL slot
            txt = "U";
        end
        hp = [hp patch(x+i,y,0,CDataMode="auto",SeriesIndex=patchColorIdx)]; % Use patch color from axis color series
        t = [t txt];
        hold on;
        text(xTxtPos+i,yTxtPos,txt,HorizontalAlignment="center")
    end
    title(join(["TDD allocation pattern," num2str(tddConfig.TDDPeriod) "slot period"]))
    axis off
    set(gca,"DataAspectRatio",[1 1 1],"YLim",[yTxtPos*2 y(3)+0.25]); % Set limits to create space between data and title
    
    pos = get(gcf,"Position");
    yScaleFactor = 3; % Scale figure to avoid excesive empty space above and below
    set(f,"Position",[pos(1) pos(2) pos(3) pos(4)/yScaleFactor]);

    hpdl = hp(find(t=="D",1)); % Handle of first downlink slot
    hpul = hp(find(t=="U",1)); % Handle of first uplink slot
    legendTxt = ["Downlink slot" "Uplink slot"];
    legendTxt = legendTxt([any(t=="D") any(t=="U")]); % Remove entry if no downlink or uplink slots
    legend([hpdl hpul],legendTxt,location="eastoutside");
end

function isActive = isSlotActive(nSlot,tddConfig)
    % Returns true if the carrier slot number contains the PDSCH

    isActive = any(mod(nSlot,tddConfig.TDDPeriod) == mod(tddConfig.SlotAllocation,tddConfig.TDDPeriod));
end

function pdschEqSym = correctCPE(carrier,pdsch,pdschInd,pdschIndInfo,rxSlotGrid,hest,nVar,pdschEqSym)
    % Correct PDSCH common phase error

    ptrsInd = hpre6GPDSCHPTRSIndices(carrier,pdsch);
    if isempty(ptrsInd)
        % No CPE correction
        return
    end

    % Map equalized PT-RS symbols to tempGrid
    tempGrid = hpre6GResourceGrid(carrier,pdsch.NumLayers);
    [ptrsRxSym,ptrsHest,~,~,~,ptrsLayerIndices] = nrExtractResources(ptrsInd,rxSlotGrid,hest,tempGrid);
    ptrsEqSym = nrEqualizeMMSE(ptrsRxSym,ptrsHest,nVar);
    tempGrid(ptrsLayerIndices) = ptrsEqSym;

    % Estimate the residual channel at the PT-RS locations in tempGrid
    ptrsSymbols = hpre6GPDSCHPTRS(carrier,pdsch);
    cpe = hpre6GChannelEstimate(carrier,tempGrid,ptrsInd,ptrsSymbols);

    % Sum estimates across subcarriers, receive antennas, and layers
    cpe = angle(sum(cpe,[1 3 4]));

    % Map equalized PDSCH symbols to tempGrid and correct CPE in each OFDM symbol
    tempGrid(pdschInd) = pdschEqSym;
    symLoc = pdschIndInfo.PTRSSymbolSet(1)+1:pdschIndInfo.PTRSSymbolSet(end)+1;
    tempGrid(:,symLoc,:) = tempGrid(:,symLoc,:).*exp(-1i*cpe(symLoc));
    pdschEqSym = tempGrid(pdschInd);
end

function m = measureEVM(refSym,rxSym)
    %   Returns a structure EVM containing:
    %     EV   - The error vector
    %     RMS  - Root Mean Square EVM (%)
    %     Peak - Peak EVM (%)

    evm = comm.EVM;
    evm.AveragingDimensions = [1 2]; % Average across subcarriers and symbols
    evm.Normalization = "Average constellation power";
    evm.AverageConstellationPower = 1; % All constellations have average unit power
    evm.MaximumEVMOutputPort = true;
    [rms,peak] = evm(refSym,rxSym);

    % Convert EVM percentages to decimal values
    m.EV = rxSym-refSym; % Error vector
    m.RMS = single(rms)/100;
    m.Peak = single(peak)/100;
end

function [evmSubcarrierRMS,evmSubcarrierPeak] = plotEVMPerSubcarrier(evGrid)
    % Plot EVM per subcarrier

    % Resource elements are NaN when they contain no PDSCH, so ignore for
    % measurement (omitmissing)
    evmSubcarrierRMS = sqrt(mean(abs(evGrid).^2,2,"omitmissing"))*100;
    evmSubcarrierPeak = sqrt(max(abs(evGrid).^2,[],2,"omitmissing"))*100; 
    numSC = height(evmSubcarrierRMS);
    figure;
    plot(0:numSC-1,evmSubcarrierRMS,'.-');
    hold on;
    plot(0:numSC-1,evmSubcarrierPeak,'.-');
    legend(["RMS" "Peak"]);
    xlabel("Subcarrier Number");
    ylabel("EVM (%)");
    grid on
    title("PDSCH EVM per subcarrier")
end

function [evmSymbolRMS,evmSymbolPeak] = plotEVMPerSymbol(evGrid)
    % Plot EVM per OFDM symbol

    % Resource elements are NaN when they contain no PDSCH, so ignore for
    % measurement (omitmissing)
    evmSymbolRMS = sqrt(mean(abs(evGrid).^2,1,"omitmissing"))*100;
    evmSymbolPeak = sqrt(max(abs(evGrid).^2,[],1,"omitmissing"))*100; 
    figure;
    plot(0:width(evmSymbolRMS)-1,evmSymbolRMS,'.-');
    hold on;
    plot(0:width(evmSymbolRMS)-1,evmSymbolPeak,'.-');
    legend(["RMS" "Peak"]);
    xlabel("Symbol Number");
    ylabel("EVM (%)");
    grid on
    title("PDSCH EVM per symbol")
end