Integrated Sensing and Communication Using 6G Waveform on NI USRP Radio

Since R2026a

This example demonstrates integrated sensing and communication (ISAC) capabilities of a candidate 6G waveform using a target emulation design deployed on the FPGA of an NI™ USRP™ radio.

Associated Radar Target Emulation Examples

This example is one in a group of related examples that demonstrate how to use a radar target emulation design running on the FPGA of an NI USRP radio in different deployment scenarios. For more information, see Radar Target Emulation Applications Overview.

Introduction

The dual use of 6G candidate signals for communication and sensing applications is of interest for use cases ranging from intruder detection in smart home or industrial applications to tracking moving objects such as unmanned aerial vehicles (UAVs). This example illustrates how you can use a candidate 6G waveform to track a moving target by modeling the signal paths using a target emulator deployed on the FPGA of an NI USRP radio.

Channel matrix estimates obtained from physical downlink shared channel (PDSCH) frames inherently capture information about time delays and Doppler shifts experienced by the transmitted waveform as it travels to the receiver. By processing this information, the receiver can form a range-Doppler plot that tracks the target movement and simultaneously calculates the error-vector magnitude (EVM).

In this example, you explore these concepts by following these steps:

1. Set up an ISAC scenario with a single transmitter and receiver that represent the base station (gNodeB) and user equipment (UE), respectively. Configure multiple targets to move in the space around the receiver. For more information on general ISAC modeling and scenario configuration, see Integrated Sensing and Communication Using 5G Waveform (Phased Array System Toolbox).

2. Configure a 6G waveform. For more information on transmitting and receiving a 6G waveform with an NI USRP radio, see Explore 400 MHz 6G Waveform Using USRP X410 (5G Toolbox).

3. Send the waveform to the target emulator logic running on the FPGA of your NI USRP radio.

4. Receive the retransmitted waveform from the target emulator and calculate the EVM and sensing information.

Prerequisites

This example uses the target emulator bitstream from the Radar Target Emulation on NI USRP Radio example. Before you run this example, follow these steps:

1. Run the Radar Target Emulation on NI USRP Radio example to generate a bitstream for the target emulation design, deploy the bitstream on the FPGA, and verify the design running on your radio.

2. Open this example in MATLAB® by using the openExample function: openExample("wt/ISACUsing6GWaveformOnNIUSRPRadioExample")

3. Copy these files into the new example working directory:

   • The hand-off information file, wtRadarTargetEmulatorSL_wthandoffinfo.mat

   • The bitstream (.bit) and device tree (.dts) file, which are generated in the project folder

   • The generated host interface setup function file, gs_wtRADARTargetEmulatorSL_setup.m

Set Up Radio

Call the radioConfigurations function. The function returns all available radio setup configurations that you saved using the Radio Setup wizard.

savedRadioConfigurations = radioConfigurations;

To update the dropdown menu with your saved radio setup configuration names, click Update. Then select the radio to use with this example.

savedRadioConfigurationNames = [string({savedRadioConfigurations.Name})];
radioName = savedRadioConfigurationNames(1) ; 

Evaluate the transmit and receive antennas available on your radio device. You use these values to configure transmit and capture antennas later in the script.

availableReceiveAntennas = hCaptureAntennas(radioName);
availableTransmitAntennas = hTransmitAntennas(radioName);

Create a radio object. This object enables you to synchronize transmit and receive operations.

radio = radioConfigurations(radioName);

Select the radio sampling rate to be a 6G sampling rate compatible with your hardware. Because the target emulator bitstream is generated with the Reference Design Optimization parameter set to High and the Sample Rate parameter set to a master clock rate (MCR) by default, you must select an MCR. For example 245.76 MHz for a USRP X410 or N320.

RadioSamplingRate = 245.76e6;

Set Up Sensing Scenario

Create a simple 2-D scenario with an update rate of 500 ms and run the scenario for 10 steps to demonstrate the target motion.

nSteps = 10;
tUpdate = 0.5;
tRun = (nSteps*tUpdate)-tUpdate;
Fc = 0.706e9;
wavelength = freq2wavelen(Fc);
scenario = radarScenario(UpdateRate=1/tUpdate,StopTime=tRun);

Set the transmitter position to the origin and the receiver position at 500 m away.

Tx = platform(scenario,Position=[0 0 0]);
Rx = platform(scenario,Position=[500 0 0]);

Set up five targets by defining starting and ending positions.

% First target
targetPosStart = [300 400 0];
targetPosEnd = [400 500 0];
targets{1} = platform(scenario);
targets{1}.Trajectory = waypointTrajectory([targetPosStart;targetPosEnd],[0 tRun]);
% Second target
targetPosStart = [600 200 0];
targetPosEnd = [300 300 0];
targets{2} = platform(scenario);
targets{2}.Trajectory = waypointTrajectory([targetPosStart;targetPosEnd],[0 tRun]);
% Third target
targetPosStart = [800 200 0];
targetPosEnd = [400 100 0];
targets{3} = platform(scenario);
targets{3}.Trajectory = waypointTrajectory([targetPosStart;targetPosEnd],[0 tRun]);
% Fourth target
targetPosStart = [-200 -200 0];
targetPosEnd = [500 -500 0];
targets{4} = platform(scenario);
targets{4}.Trajectory = waypointTrajectory([targetPosStart;targetPosEnd],[0 tRun]);
% Fifth target
targetPosStart = [800 800 0];
targetPosEnd = [700 700 0];
targets{5} = platform(scenario);
targets{5}.Trajectory = waypointTrajectory([targetPosStart;targetPosEnd],[0 tRun]);

Plot the transmitter and receiver positions, the starting positions, and the paths of the moving targets.

hISACPlotScenario(scenario);

Configure 6G Waveform Parameters

Set up 6G waveform parameters, with control over the bandwidth and subcarrier spacing.

bandwidth = 100; % MHz
subcarrierSpacing = 15; % kHz

Create a carrier configuration by using the pre6GCarrierConfig (5G Toolbox) object. Calculate the required number of resource blocks (NSizeGrid) given the bandwidth.

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

Specify the number of subframes to generate.

numSubframes = 10;
numSlots2Generate = numSubframes*carrier.SlotsPerSubframe;
disp("Generating "+numSlots2Generate+" slots")

Generating 10 slots

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 = "QPSK";

6G Waveform Configuration for Sensing

PDSCH frames include demodulation reference signal (DM-RS) symbols that help the receiver estimate the channel characteristics. The communication function uses this information to correctly demodulate the transmitted data. The sensing function can use the same information to estimate positions and velocities of the targets. The configuration of DM-RS pilot signals directly influences the sensing performance. For more information about DM-RS configurations, see the NR PDSCH Resource Allocation and DM-RS and PT-RS Reference Signals (5G Toolbox) example.

Set the DM-RS parameters.

pdsch.DMRS.DMRSTypeAPosition = 2;           % Mapping type A only. First DM-RS symbol position is 2
pdsch.DMRS.DMRSLength = 1;                  % Single-symbols DM-RS
pdsch.DMRS.DMRSAdditionalPosition = 2;      % Additional DM-RS symbol positions (max range 0...3)
pdsch.DMRS.NumCDMGroupsWithoutData = 1;     % Number of CDM groups without data
pdsch.DMRS.DMRSConfigurationType = 1;       % DM-RS configuration type (1,2)

Visualize one resource block of the configured resource grid.

resourceGrid = hpre6GResourceGrid(carrier, pdsch.NumLayers);
ind = nrPDSCHDMRSIndices(carrier, pdsch);
resourceGrid(ind) = nrPDSCHDMRS(carrier, pdsch);

visualizeResourceGrid(resourceGrid(1:12,1:14,1));
title({'PDSCH DM-RS Resource Elements in the Carrier Resource Grid', '(single resource block)'});

The density of pilot signals in the time domain determines the maximum Doppler shift that can be measured unambiguously, without causing aliasing. In the configured PDSCH frame, the maximum separation between pilot signals in the time domain is $M_t$ = 5 symbols. The corresponding maximum unambiguous Doppler shift can be computed using the formula

$$f_{\text{Dmax}}\le \frac{1}{2M_t{T}_{\text{OFDM}}}$$

where $T_{\text{OFDM}}$is the symbol duration.

% Compute maximum pilot signal separation in time
[~,pdschInfo] = hpre6GPDSCHIndices(carrier,pdsch);
Mt = max(diff(pdschInfo.DMRSSymbolSet));

% Duration of an OFDM symbol in seconds
ofdmInfo = hpre6GOFDMInfo(carrier,Nfft=RadioSamplingRate/(subcarrierSpacing*1e3));
ofdmSymbolDuration = mode(ofdmInfo.SymbolLengths)/ofdmInfo.SampleRate;

maxDoppler = 1/(2*Mt*ofdmSymbolDuration);
fprintf('Maximum unambiguous Doppler shift: %.2f kHz\n', maxDoppler*1e-3);

Maximum unambiguous Doppler shift: 1.40 kHz

Compute the corresponding maximum velocity and compare it to the maximum ground speeds of the targets in the considered ISAC scenario.

maxVelocity = dop2speed(maxDoppler, wavelength)/2;
fprintf('Maximum unambiguous velocity: %.2f m/s\n', maxVelocity);

Maximum unambiguous velocity: 297.55 m/s

fprintf("Target velocities are: %s m/s",join(string(cell2mat(cellfun(@(x) norm(x.Trajectory.Velocities(1,:)),targets,'UniformOutput',false))),', '))

Target velocities are: 31.42697, 70.27284, 91.62457, 169.2394, 31.42697 m/s

The velocities of the targets in the scenario are within the maximum unambiguous velocity.

Similarly, the density of pilot signals in the frequency domain determines the maximum time delay that can be measured without causing aliasing. In the frequency domain the pilot signals repeat every $M_f$= 2 subcarriers. The corresponding maximum unambiguous time delay can be computed using the formula

$${\tau }_{\max }\le \frac{1}{2M_f\Delta f}$$

where $\Delta f$is the subcarrier spacing.

% Compute pilot signal separation in frequency
Mf = min(diff(pdsch.DMRS.DMRSSubcarrierLocations));
maxDelay = 1/(2*Mf*subcarrierSpacing*1e3);
fprintf('Maximum time delay: %.2fe-6 s\n', maxDelay*1e6);

Maximum time delay: 16.67e-6 s

Compute the corresponding maximum range.

maxRange = time2range(maxDelay)*2; % time2range assumes 2-way propagation
fprintf('Maximum unambiguous range: %.2f m\n', maxRange);

Maximum unambiguous range: 4996.54 m

All targets in the scenario are within the maximum unambiguous range.

Finally, calculate the range resolution, which is determined by the transmission bandwidth.

% Range resolution
transmissionBandwidth = subcarrierSpacing*1e3*carrier.NSizeGrid*12; % Bandwidth of resource blocks in Hz
rangeResolution = bw2rangeres(transmissionBandwidth);
fprintf('Range resolution: %.2f m\n', rangeResolution);

Range resolution: 1.50 m

Generate Waveform

Create a resource grid containing a PDSCH transmission, DM-RS, and PT-RS.

waveGrid = []; % Store resource grid to transmit
for nSlot = 0:(numSlots2Generate-1)
    carrier.NSlot = nSlot;
    txSlotGrid = hpre6GResourceGrid(carrier,pdsch.NumLayers);
    [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;

    % Generate and map DM-RS
    dmrsInd = hpre6GPDSCHDMRSIndices(carrier,pdsch);
    dmrsSym = hpre6GPDSCHDMRS(carrier,pdsch);
    txSlotGrid(dmrsInd) = dmrsSym;
    waveGrid = [waveGrid txSlotGrid]; %#ok<AGROW>
end

OFDM-modulate the resource grid and plot the spectrum. This is the 6G candidate waveform that you transmit.

carrier.NSlot = 0; % OFDM-Modulate from slot 0
[txWaveform,ofdmInfo] = hpre6GOFDMModulate(carrier,waveGrid,Nfft=RadioSamplingRate/(subcarrierSpacing*1e3));
% Scale the transmit waveform for the hardware.
txWaveform = txWaveform / max(abs(txWaveform));
txWaveform = 0.7*txWaveform*double(intmax('int16'));
txSA = spectrumAnalyzer;
txSA.SampleRate = ofdmInfo.SampleRate;
txSA.Title = "Generated 6G Candidate Waveform";
txSA(txWaveform)

Set Up Radio

Create a usrp System object™. This System object enables you to configure the radio front end properties, stream samples from the radio front end to your DUT, and transmit and receive IQ data to and from the host.

device = usrp(radio);

If you have not yet programmed your device with the bitstream, select the program bitstream option. Update the code with your bitstream and hand-off information files that you generated during the radar target emulation bitstream build. If your radio device is a USRP X310, you need only a bitstream file to program the FPGA.

programBitstream = false;
if(programBitstream)
programFPGA(device, ...
    "x4xx.bit", ...  % replace with your .bit file
    "usrp_x410_fpga_X1_200.dts");  % replace with your .dts file
end

Use the describeFPGA function to configure the DUT interfaces according to the hand-off information file.

% replace with your .mat handoff file
describeFPGA(device,"wtRADARTargetEmulatorSL_wthandoffinfo.mat"); 

Configure Device

Configure your usrp System object to loop back the waveform through the target emulator. You can loop back over the air or directly on the FPGA.

Set the sample rate.

device.SampleRate = ofdmInfo.SampleRate;

Set the LoopbackMode property to FPGA to loop back each transmit antenna to an associated receive antenna on the FPGA. Alternatively, you can select Disabled to transmit and receive data over the air.

device.LoopbackMode = "FPGA";

Because the bitstream is generated with the Reference Design Optimization parameter set to High, the radio antennas connected to the DUT are fixed to the options you configured in the Interface Mapping table in Simulink®. To use FPGA loopback, your DUT must be configured with an input and output antenna connection that is a valid loopback antenna pair. For details, see the table in the LoopbackMode property description.

Specify a transmit antenna and a capture antenna to exercise the target emulation DUT from MATLAB. Since the bitstream is built with the Reference Design Optimization set to High, only one transmit and capture antenna is available, respectively.

% MATLAB IO Antennas
device.TransmitAntennas = availableTransmitAntennas(1);
device.CaptureAntennas = availableReceiveAntennas(2);

Specify the center frequency and gain for the selected transmit and receive antennas. To loop back over the air, ensure the transmit and receive center frequencies are equal. If FPGA loopback is enabled, these properties have no effect. To isolate the two pairs of antennas, either use loopback cables or separate the antennas in frequency by specifying multiple center frequency values. If you have a USRP N310 radio, specifying multiple center frequencies is only possible using independent channels. For details, see the ReceiveCenterFrequency and TransmitCenterFrequency property descriptions.

device.ReceiveCenterFrequency = [Fc,Fc+200e6];
device.TransmitCenterFrequency = [Fc+200e6,Fc];

Specify the desired RF gain for the target emulator input and output. The output gain should be specified as a reference corresponding to the maximum output power of the target emulator in the initial state.

targetEmulatorInputRadioGain = 30; % Change according to your RF setup and hardware
targetEmulatorOutputRadioGainReference = 30;

Specify the radio gains for the transmit and capture antennas.

mlCaptureRadioGain = 30; % Change according to your RF setup and hardware
mlTransmitRadioGain = 30;

Specify the radio gains which are applied in the radio front end for each antenna. The DUT input and output antennas are specified first, followed by the capture and transmit antennas.

device.ReceiveRadioGain = [targetEmulatorInputRadioGain, mlCaptureRadioGain];
device.TransmitRadioGain = [targetEmulatorOutputRadioGainReference, mlTransmitRadioGain];

Create and Set Up fpga Object

To interface with the target emulator, create an fpga object with the device System object.

device.reset;
dut = fpga(device);

Set up the fpga object using the setup function you generated in the Simulink toolstrip workflow.

gs_wtRADARTargetEmulatorSL_setup(dut);

Set Up Device

To establish a connection with the radio hardware, call setup on the device System object.

setup(device);

Prepare Scenario and Configure Initial Target Emulator State

Specify the front end delay so that it can be accounted for in the simulation. When using FPGA loopback mode, set the front end delay to zero. If the front end delay is not accounted for when you transmit and receive the data over the air, there is a small shift in the estimated ranges. Calibrate the additional delay for your radio by doing a timed loopback.

calibratedFrontEndDelay = 0;
reTxOffset = 512;

Create a phased.ScatteringMIMOChannel (Phased Array System Toolbox) object to model the channel between the transmitter and receiver, including the targets. The channel models the direct path between the transmitter and the receiver as well as the paths from the transmitter to the receiver through the targets.

channel = hISACChannel(scenario,Fc,ofdmInfo.SampleRate);

Start the scenario and map the channel characteristics for the current scenario state to the target emulator.

restart(scenario);
advance(scenario);
[TargetRange, TargetDoppler, TargetEnabled] = hISACUpdateTargetEmulator(dut, scenario, channel, ofdmInfo, calibratedFrontEndDelay, reTxOffset);

Prepare for Target Emulation

Configure the receive waveform processing. Set do6GEVM to calculate an EVM metric. Set displayConst to display a constellation plot of the received signal.

do6GEVM = false;
displayConst = false;
numFramesToCapture = 2;
captureLength = 10e-3*ofdmInfo.SampleRate*numFramesToCapture;

Prepare 6G timing synchronization signal.

refSlotGrid = waveGrid(:,1:carrier.SymbolsPerSlot,:);
[refWaveform,~] = hpre6GOFDMModulate(carrier,refSlotGrid,Nfft=RadioSamplingRate/(subcarrierSpacing*1e3));

Set up a constellation diagram.

constDiagram = [];
if displayConst
    constDiagram = comm.ConstellationDiagram;
    constDiagram.Title = "Equalized PDSCH Constellation";
    constDiagram.ShowReferenceConstellation = true;
    modList = {'QPSK','16QAM','64QAM','256QAM','1024QAM','4096QAM'};numSymbolsList = [4 16 64 256 1024 4096];
    numSym = numSymbolsList(strcmpi(modList,pdsch.Modulation));
    constDiagram.ReferenceConstellation = qammod((0:numSym-1)',numSym,UnitAveragePower=1);   
end

Calculate the times at which the individual symbols within a frame are transmitted. This information is used for sensing.

tau = cumsum([0 repmat(ofdmInfo.SymbolLengths, 1, ofdmInfo.SlotsPerFrame/ofdmInfo.SlotsPerSubframe)])/ofdmInfo.SampleRate;
tau = tau(1:end-1);

Run the Scenario

Prepare a figure to display the range-Doppler plot.

sensingFigure = figure;
sensingAxes = axes(sensingFigure);
title('Estimated Range-Doppler Map');
xlabel('Frequency (kHz)')
ylabel('Bistatic Range (m)')
set(sensingAxes, 'Ydir', 'normal')
hcb = colorbar;
hcb.Label.String = '(dB)';
ylim([0 2000]);
xlim([-1 1]);

Transmit the 6G candidate waveform continuously from the host.

transmit(device,txWaveform,"continuous");

Run the scenario. For each loop, capture the output of the target emulator and process the data with sensing and communication functions. Display the range-Doppler plot of the targets.

% Stream captureLength samples through the target emulator
device(captureLength);

% Run the scenario
while advance(scenario)

    % Capture the data to the host for processing
    [rxWaveform,nSamps] = capture(device,captureLength);
    if nSamps==0
        error('No data received');
    end
    if ~nnz(rxWaveform)
        error('Waveform all zeros');
    end
    rxWaveform = single(rxWaveform);

    % Timing synchronization
    [offset, corro] = timingEstimate(rxWaveform,refWaveform);
    rxWaveform = rxWaveform(1+offset:end,:);

    % Demodulate
    rxWaveGrid = hpre6GOFDMDemodulate(carrier,rxWaveform,Nfft=RadioSamplingRate/(subcarrierSpacing*1e3));

    numRxSlots = floor(size(rxWaveGrid,2)/carrier.SymbolsPerSlot);
    if numRxSlots<1
        error("Not enough data to measure EVM, increase numSubframes")
    end
    slots2decode = min(numSlots2Generate,numRxSlots);

    numSubcarriers = carrier.NSizeGrid*12;
    numSymbolsPerFrame = carrier.SlotsPerFrame*carrier.SymbolsPerSlot;

    timeSynchronizedChannelMatrix = zeros(numSubcarriers, numSymbolsPerFrame);
    % Decode each slot
    for nSlot = 0:(slots2decode-1)
        carrier.NSlot = nSlot;

        % Extract slot from grid
        rxSlotGrid = rxWaveGrid(:,nSlot*carrier.SymbolsPerSlot+(1:carrier.SymbolsPerSlot),:);

        % Estimate channel
        pdschDMRSInd = hpre6GPDSCHDMRSIndices(carrier,pdsch);
        pdschDMRSSymbols = hpre6GPDSCHDMRS(carrier,pdsch);
        [hest,nVar] = hpre6GChannelEstimate(carrier,rxSlotGrid,pdschDMRSInd,pdschDMRSSymbols,CDMLengths=pdsch.DMRS.CDMLengths);
        
        % Decode data and calculate and print EVM
        if do6GEVM
            [slotEVM] = get6GMetrics(constDiagram,carrier,pdsch,rxSlotGrid,hest,nVar);
            fprintf("\nSlot %d PDSCH EVM, RMS: %0.3f%%, Peak: %0.3f%%",nSlot,slotEVM.RMS*100,slotEVM.Peak*100);
        end

        % Sensing
        timeSynchronizedChannelMatrix(:, nSlot*carrier.SymbolsPerSlot+1:(nSlot+1)*carrier.SymbolsPerSlot) = hest;
    end
    % To filter the static scatterers out, put a zero in the first Doppler
    X = fft(timeSynchronizedChannelMatrix, [], 2);
    X(:, 1) = 0;
    x = ifft(X, [], 2);
    title(sensingAxes,"Estimated Range-Doppler Map")
    hISACRangeDopplerMap(x, subcarrierSpacing * 1e3, tau, sensingAxes, TargetRange, TargetDoppler, TargetEnabled);

    [TargetRange, TargetDoppler, TargetEnabled] = hISACUpdateTargetEmulator(dut, scenario, channel, ofdmInfo, calibratedFrontEndDelay, reTxOffset);
    device(captureLength);
end

release(device);
release(dut);

Local Functions

function evm = get6GMetrics(constDiagram,carrier,pdsch,rxSlotGrid,hest,nVar)
    % 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;
    
    % 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
    if ~isempty(constDiagram)
        constDiagram(rxSym{1});
    end
    
    % % Calculate error vector and EVM
    evm = measureEVM(refSym{1},rxSym{1});
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 visualizeResourceGrid(G)

x = 0:(size(G, 2)-1);
y = 0:(size(G, 1)-1);

fig = figure;
ax = axes(fig);
imagesc(ax, x-0.5, y-0.5, abs(G));
grid(ax, 'on'); 
ax.GridColor=[1 1 1];
ax.YDir = 'normal';

ax.XTick = x-1;
ax.YTick = y-1;

ax.XTickLabel = compose('%d', x);
ax.YTickLabel = compose('%d', y);

xlabel(ax, 'OFDM Symbols');
ylabel(ax, 'Subcarriers');
end