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

carrierFrequency = 140; % GHz
bandwidth = 2160; % MHz - transmission bandwidth and guardbands
subcarrierSpacing = 3840; % kHz
waveform = "OFDM"; % Waveform modulation

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

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

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("Simulating "+numSlots+" slots")

visualizeTDDAllocation(tddConfig);

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

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

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

        % Generate and map PT-RS
        ptrsInd = hpre6GPDSCHPTRSIndices(carrier,pdsch);
        ptrsSym = hpre6GPDSCHPTRS(carrier,pdsch);
        txSlotGrid(ptrsInd) = ptrsSym;

        if waveform=="DFT-s-OFDM"
            txSlotGrid = transformPrecode(pdsch,txSlotGrid,pdschInd,ptrsInd);
        end
        
        % Generate and map DM-RS
        dmrsInd = hpre6GPDSCHDMRSIndices(carrier,pdsch);
        dmrsSym = hpre6GPDSCHDMRS(carrier,pdsch);
        txSlotGrid(dmrsInd) = dmrsSym;
    end
    waveGrid = [waveGrid txSlotGrid];
end

% 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 for filter delay
waveformLength = height(txWaveform);
slotLength = waveformLength/numSlots;
txWaveform = [txWaveform(end-slotLength+1:end,:); txWaveform];

impariredWaveform = txWaveform;
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
    impariredWaveform = phaseNoise(impariredWaveform); % Apply phase noise model
end

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
    impariredWaveform = LPF(impariredWaveform);
    release(LPF);
end

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

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(impariredWaveform);
    else
        warning("Sample rate too low to measure ACPR, increase oversamplingFactor")
    end
end

measureCCDF = true;
if measureCCDF
    pm = powermeter(ComputeCCDF=true);
    averagePower = pm(impariredWaveform);
    disp("Average power: "+averagePower+" dBm")
    figure;
    plotCCDF(pm,GaussianReference=true)
end

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

% Timing synchronization
refSlotGrid = waveGrid(:,1:carrier.SymbolsPerSlot,:);
refWaveform = hpre6GOFDMModulate(carrier,refSlotGrid,Nfft=ofdmInfo.Nfft);
offset = timingEstimate(rxWaveform,refWaveform);
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
cpe = [];
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 and PT-RS
        [pdschInd,pdschIndInfo] = hpre6GPDSCHIndices(carrier,pdsch);
        ptrsInd = hpre6GPDSCHPTRSIndices(carrier,pdsch);
        ind = [pdschInd; ptrsInd];
        [rxSym,rxSymHest] = nrExtractResources(ind,rxSlotGrid,hest);
        eqSym = nrEqualizeMMSE(rxSym,rxSymHest,nVar);

        eqSlotGrid = hpre6GResourceGrid(carrier,pdsch.NumLayers);
        eqSlotGrid(ind) = eqSym;

        if waveform=="DFT-s-OFDM"
            eqSlotGrid = transformDeprecode(pdsch,eqSlotGrid,pdschInd,ptrsInd);
        end

        % Correct common phase error
        eqSlotGrid = correctCPE(carrier,pdsch,pdschIndInfo,eqSlotGrid);

        % Demap, decode, and re-modulate to obtain layer demapped equalized
        % symbols (rxSym) and reference symbols (refSym) for EVM
        % measurement.
        pdschEqSym = eqSlotGrid(pdschInd);
        [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 = "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

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

evmSubcarrierRMS = plotEVMPerSubcarrier(evGrid);

evmSymbolRMS = plotEVMPerSymbol(evGrid);

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("AM/AM, "+modelName)
    legend(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 excessive 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 slotGrid = transformPrecode(pdsch,slotGrid,pdschInd,ptrsInd)
    % Perfrom transform precoding

    % Apply scaling to PTRS symbols
    ptrsSym = slotGrid(ptrsInd)*hPre6GPTRSPowerFactorDFTsOFDM(pdsch.Modulation);
    
    % Extract only PDSCH and PTRS and sort into absolute frequency/symbol
    % number
    ind = [pdschInd; ptrsInd];
    [sortedCombInd,sortIdx] = sort(ind);
    sym = [slotGrid(pdschInd); ptrsSym];
    sortedSym = sym(sortIdx);
    
    MRB = numel(pdsch.PRBSet);
    tp = nrTransformPrecode(sortedSym,MRB);
    slotGrid(sortedCombInd) = tp;
end

function slotGrid = transformDeprecode(pdsch,slotGrid,pdschInd,ptrsInd)
    % Perform transform de-precoding

    % De-precode
    ind = [pdschInd; ptrsInd];
    sortedInd = sort(ind);
    eqSym = slotGrid(sortedInd);
    MRB = numel(pdsch.PRBSet);
    detransformed = nrTransformDeprecode(eqSym,MRB);
    slotGrid(sortedInd) = detransformed;

    % Undo scaling on PT-RS
    slotGrid(ptrsInd) = slotGrid(ptrsInd)/hPre6GPTRSPowerFactorDFTsOFDM(pdsch.Modulation);
end

function [eqSymGrid,cpe] = correctCPE(carrier,pdsch,pdschIndInfo,eqSymGrid)
    % Correct PDSCH common phase error

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

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

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

    % Correct CPE in each OFDM symbol
    symLoc = pdschIndInfo.PTRSSymbolSet(1)+1:pdschIndInfo.PTRSSymbolSet(end)+1;
    eqSymGrid(:,symLoc,:) = eqSymGrid(:,symLoc,:).*exp(-1i*cpe(symLoc));
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