exponenta event banner

Генерация и передача сигналов 5G с использованием испытательного и измерительного оборудования

Этот пример показывает, как произвести и передать беспроводное 5G использование формы волны 5G Toolbox™, Контроль за Инструментом инструменты Keysight Technologies® RF и Toolbox™.

Введение

В этом примере взаимодействие с генератором радиочастотных сигналов и анализатором с помощью инструментария управления приборами. Затем создайте форму сигнала 5G NR-TM с помощью приложения Wireless Waveform Generator из панели инструментов 5G. Загрузите сгенерированный сигнал в генератор сигналов Keysight Technologies E4438C для передачи по воздуху. Наконец, зафиксируйте сигнал в воздухе с помощью анализатора сигналов Keysight Technologies N9030A и проанализируйте его в MATLAB.

Требования

Для выполнения этого примера необходимо:

  • Keysight Technologies ® E4438C генератор сигналов

  • Анализатор сигналов Keysight Technologies ® N9030A

  • 5G Toolbox™

  • Toolbox™ управления приборами

  • Toolbox™ обработки сигналов

  • Контроль за инструментом пакет поддержки Toolbox™ для библиотек Keysight™ IO и интерфейса ВИЗЫ

Создание формы сигнала основной полосы частот с помощью приложения Wireless Waveform Generator

Создание формы сигнала 5G NR-TM с помощью приложения Wireless Waveform Generator App. Поскольку генерация формы сигнала основана на приложении, вам не нужно писать код MATLAB.

Дополнительные сведения о приложении Wireless Waveform Generator см. в разделе 5G Waveform Generator.

На панели инструментов MATLAB Apps выберите Беспроводное Приложение для Генератора Формы волны. Откройте Беспроводное приложение для Генератора Формы волны и настройте, чтобы произвести 5G НОМЕР ТМ (Экспериментальная модель) форма волны. Для использования этой функции приложения требуется 5G Toolbox™.

В разделе «Тип сигнала» выберите «Тестовые модели» (NR-TM). На левой панели приложения на вкладке Форма волны можно задать параметры выбранной формы волны. Для этого примера:

  • Установите диапазон частот FR1 (450 МГц - 6 ГГц).

  • Установите тестовую модель как полную полосу, однородные 64QAM, представленные как NR-FR1-TM3.1.

  • Установите полосу пропускания канала 10 МГц, интервал между поднесущими 30 кГц и дуплексный режим FDD.

  • Щелкните Создать (Generate).

Как и ожидалось, полоса пропускания сигнала 10 МГц видна на основной полосе частот.

Передача сигнала по воздуху

После генерации сигнала основной полосы подключитесь к генератору радиочастотного сигнала через один из поддерживаемых коммуникационных интерфейсов с помощью инструментария управления приборами. В приложении Wireless Waveform Generator перейдите на вкладку «Передатчик» для передачи сгенерированного сигнала по воздуху. Приложение автоматически находит генератор сигналов, подключенный через интерфейс TCP/IP через кабель LAN. Выберите драйвер SCPI генератора сигналов Agilent/Keysight в раскрывающемся меню Driver. Установите центральную частоту 3,4 ГГц и выходную мощность -15 дБм. Частота дискретизации основной полосы уже вычислена на основе генерируемого сигнала в MATLAB. Для передачи нажмите кнопку Transmit Waveform.

Считывание синфазных и квадратурных (IQ) данных из анализатора сигналов по TCP/IP

Для анализа передачи по воздуху в MATLAB используется инструментарий управления приборами для настройки анализатора сигналов Keysight Technologies N9030A и считывания синфазных и квадратурных фазовых данных.

Определите параметры, используемые для конфигурирования прибора перед выполнением измерений. На основе измеряемого сигнала может потребоваться изменить некоторые из следующих параметров.

% Center frequency of the modulated waveform (Hz)
centerFrequency = 3.4e9; 

% Bandwidth of the signal (Hz)
bandwidth = 12.287999e6;

% Measurement time (s)
measureTime = 20e-3;

% Mechanical attenuation in the signal analyzer (dB)
mechAttenuation = 0;

% Start frequency for Spectrum Analyzer 
startFrequency = 3.39e9; 

% Stop frequency for Spectrum Analyzer  
stopFrequency = 3.41e9;

% Resolution Bandwidth for Spectrum Analyzer 
resolutionBandwidth = 220e3; 

% Video Bandwidth for Spectrum Analyzer view
videoBandwidth = 220000;

Перед подключением к анализатору спектра выполните следующие действия:

  • Настройка подключения КИП с использованием TCP/IP-соединения.

  • Отрегулируйте размер входного буфера так, чтобы он мог удерживать данные, возвращаемые прибором.

  • Установите время ожидания, достаточное для измерения и передачи данных.

  • Подсоедините к прибору.

sigAnalyzerObj = visadev("TCPIP0::A-N9030A-71181.dhcp.mathworks.com::inst0::INSTR");
sigAnalyzerObj.ByteOrder = "big-endian";
configureCallback(sigAnalyzerObj,"byte", 85e6, @callbackFcn) 
sigAnalyzerObj.Timeout = 20;

Сбросьте прибор в известное состояние с помощью соответствующей команды SCPI. Запросите идентификационные данные прибора, чтобы убедиться, что нужный прибор подключен.

writeline(sigAnalyzerObj, "*RST");
instrumentInfo = writeread(sigAnalyzerObj, "*IDN?");
fprintf("Instrument identification information: %s", instrumentInfo);
Instrument identification information: Agilent Technologies,N9030A,US00071181,A.14.16

Анализаторы сигналов и спектров серии X выполняют IQ-измерения, а также спектральные измерения. В этом примере выполняется сбор данных IQ временной области, визуализация данных с помощью MATLAB и анализ полученных данных. Команды SCPI используются для конфигурирования прибора и определения формата передачи данных после выполнения измерений.

% Set up signal analyzer mode to Basic IQ mode
writeline(sigAnalyzerObj,":INSTrument:SELect BASIC");

% Set the center frequency
writeline(sigAnalyzerObj,[':SENSe:FREQuency:CENTer %s' num2str(centerFrequency)]);

% Set the resolution bandwidth
writeline(sigAnalyzerObj,[':SENSe:WAVEform:BANDwidth:RESolution ' num2str(bandwidth)]);

% Turn off averaging
writeline(sigAnalyzerObj,":SENSe:WAVeform:AVER OFF");

% Set to take one single measurement once the trigger line goes high
writeline(sigAnalyzerObj,":INIT:CONT OFF"); 

% Set the trigger to external source 1 with positive slope triggering
writeline(sigAnalyzerObj,":TRIGger:WAVeform:SOURce IMMediate");
writeline(sigAnalyzerObj,":TRIGger:LINE:SLOPe POSitive");

% Set the time for which measurement needs to be made
writeline(sigAnalyzerObj,[':WAVeform:SWE:TIME '  num2str(measureTime)]);

% Turn off electrical attenuation 
writeline(sigAnalyzerObj,":SENSe:POWer:RF:EATTenuation:STATe OFF");

% Set mechanical attenuation level
writeline(sigAnalyzerObj,[':SENSe:POWer:RF:ATTenuation ' num2str(mechAttenuation)]);

% Turn IQ signal ranging to auto
writeline(sigAnalyzerObj,":SENSe:VOLTage:IQ:RANGe:AUTO ON");

% Set the endianness of returned data
writeline(sigAnalyzerObj,":FORMat:BORDer NORMal");

% Set the format of the returned data
writeline(sigAnalyzerObj,":FORMat:DATA REAL,64");

Запустите прибор для выполнения измерения. Дождитесь завершения операции измерения и прочтите сигнал. Перед обработкой данных отделите компоненты I&Q от перемеженных данных, возвращаемых прибором, и создайте сложный вектор в MATLAB.

% Trigger the instrument and initiate measurement
writeline(sigAnalyzerObj,"*TRG");
writeline(sigAnalyzerObj,":INITiate:WAVeform");

% Wait till measure operation is complete
measureComplete = writeread(sigAnalyzerObj,"*OPC?");

% Read the IQ data
writeline(sigAnalyzerObj,":READ:WAV0?");
data = readbinblock(sigAnalyzerObj,"double");

% Separate the data and build the complex IQ vector
inphase = data(1:2:end);
quadrature = data(2:2:end);
IQData = inphase+1i*quadrature;

Инструмент предоставляет информацию о последних полученных данных. Зафиксируйте эту информацию и отобразите ее.

writeline(sigAnalyzerObj,":FETCH:WAV1?");
signalSpec = readbinblock(sigAnalyzerObj,"double");

% Display the measurement information
sampleRate = 1/signalSpec(1);
fprintf("Sample Rate (Hz) = %s", num2str(sampleRate));
Sample Rate (Hz) = 15359998.75
fprintf("Number of points read = %s", num2str(signalSpec(4)));
Number of points read = 307200
fprintf("Max value of signal (dBm) = %s", num2str(signalSpec(6)));
Max value of signal (dBm) = -39.4151
fprintf("Min value of signal (dBm) = %s", num2str(signalSpec(7)));
Min value of signal (dBm) = -104.1205

Постройте первые 1000 точек полученных данных временной области и аннотируйте рисунок.

plot(real(IQData(1:1000)),'b'); 
hold on;
plot(imag(IQData(1:1000)),'g');
legend('Inphase signal', 'Quadrature signal');
title('IQ Data for the first 1000 points of acquired signal');
xlabel('Sample number');
ylabel('Voltage');
hold off;

Представление спектра может содержать больше информации, чем представление данных во временной области. Например, можно использовать представление спектра для идентификации основных полос частот, полосы пропускания сигнала и т.д. Для печати вида спектра необходимо иметь панель инструментов обработки сигналов.

% Plot the PSD of the acquired signal
periodogram(IQData, hamming(length(IQData)),[], sampleRate, 'centered');

Переведите прибор в режим анализатора спектра и сравните вид спектра, сформированный в MATLAB, с видом на анализаторе сигнала. Дополнительные команды SCPI используются для настройки параметров измерения и отображения приборов.

% Switch back to the spectrum analyzer view
writeline(sigAnalyzerObj,':INSTrument:SELect SA');

% Set mechanical attenuation level
writeline(sigAnalyzerObj,[':SENSe:POWer:RF:ATTenuation ' num2str(mechAttenuation)]);

% Set the center frequency, RBW and VBW and trigger
writeline(sigAnalyzerObj,[':SENSe:FREQuency:CENTer ' num2str(centerFrequency)]);
writeline(sigAnalyzerObj,[':SENSe:FREQuency:STARt ' num2str(startFrequency)]);
writeline(sigAnalyzerObj,[':SENSe:FREQuency:STOP ' num2str(stopFrequency)]);
writeline(sigAnalyzerObj,[':SENSe:BANDwidth:RESolution ' num2str(resolutionBandwidth)]);
writeline(sigAnalyzerObj,[':SENSe:BANDwidth:VIDeo ' num2str(videoBandwidth)]);

% Continuous measurement
writeline(sigAnalyzerObj,':INIT:CONT ON'); 

% Trigger
writeline(sigAnalyzerObj,'*TRG');

Для очистки КИП очистите соединение КИП.

% Clear instrument connection
clear sigAnalyzerObj;

Нажмите на кнопку Stop Transmission в Приложении для Генератора Формы волны, чтобы остановиться 5G НОМЕР передачи формы волны ТМ

Измерьте EVM 5G НОМЕР формы волны ТМ

После получения полученных данных IQ основной полосы из анализатора сигналов можно визуализировать, декодировать и анализировать данные с помощью 5G Toolbox.

% Ensure the IQData is a column vector
if ~iscolumn(IQData)
    IQData = IQData(:);
end

Каждый сигнал NR-TM определяется комбинацией:

  • Имя NR-TM

  • Полоса пропускания канала

  • Интервал между поднесущими

  • Дуплексный режим

% Select the NR-TM simulation parameters
nrtm = "NR-FR1-TM3.1"; % reference channel
bw   = "10MHz"; % channel bandwidth
scs  = "30kHz"; % subcarrier spacing
dm   = "FDD";   % duplexing mode

Генерация форм сигнала TM

Комбинация полосы пропускания канала и интервала между поднесущими должна быть действительной парой из соответствующей таблицы конфигурации полосы пропускания FR.

% Create generator object for the selected NR-TM
tmwavegen = hNRReferenceWaveformGenerator(nrtm,bw,scs,dm);

Формирование формы сигнала txWaveform. По умолчанию один кадр создается для FDD, а два - для TDD. tmwaveinfo вывод содержит информацию о сформированном сигнале. resourcesinfo вывод содержит информацию о ресурсах, выделенных PDSCH и PDSCH DM-RS в сгенерированном сигнале.

[txWaveform,tmwaveinfo,resourcesinfo] = generateWaveform(tmwavegen);
txWaveform = IQData(1:2*numel(txWaveform));

Сгенерированный сигнал TM может содержать более одного PDSCH. Функция hListTargetPDSCHs выбирает целевые PDSCH для анализа. Это основано на RNTI. По умолчанию для расчета EVM используется следующий RNTI:

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

Эти значения могут быть переопределены. Функция hListTargetPDSCHs также принимает третий параметр с набором RNTI для рассмотрения для измерения EVM.

pdschArray - массив структур конфигурации PDSCH для анализа. Каждая структура содержит набор ресурсов, используемых для PDSCH и связанных DM-RS и PT-RS. Он также содержит набор параметров для генерации PDSCH.

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

Ухудшение: фазовый шум и нелинейность

Фазовый шум и нелинейность без памяти могут быть включены или отключены путем переключения флагов phaseNoiseOn и nonLinearityModelOn.

phaseNoiseOn = false;
nonLinearityModelOn = false;

Нормализуйте форму сигнала в соответствии с динамическим диапазоном нелинейности.

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

Приемник

Приемник в этом примере выполняет следующие шаги:

  • Синхронизация с использованием DM-RS через один кадр для FDD (два кадра для TDD)

  • OFDM-демодуляция принятого сигнала

  • Оценка канала

  • Уравнивание

Получение параметров формы сигнала для синхронизации и демодуляции OFDM

gnb.NRB = tmwavegen.Config.SCSCarriers{1}.NSizeGrid;
gnb.CyclicPrefix = tmwavegen.Config.BandwidthParts{1}.CyclicPrefix;
gnb.SubcarrierSpacing = tmwavegen.Config.SCSCarriers{1}.SubcarrierSpacing;

Компенсация грубого смещения несущей

DMRS, поиск смещений с приращением 10kHz до 200kHz

peak = [];
t = 0:numel(txWaveform)-1;
foffsetRange = -20:20;
for foffsetnr = foffsetRange
    % Correction
    phase = 2*pi*foffsetnr*1e4*t'/tmwaveinfo.Info.SamplingRate;
    rxWaveform = txWaveform .* exp(1j*phase);

    % 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);
    peak = [peak; max(mag)];
end
[~,ind] = max(peak);
coarseOffset = foffsetRange(ind);
fprintf('Coarse frequency offset = %.0f *10kHz\n', coarseOffset)
Coarse frequency offset = 0 *10kHz
% Correction
phase = 2*pi*coarseOffset*1e4*t'/tmwaveinfo.Info.SamplingRate;
rxWaveform = txWaveform .* exp(1j*phase);

Более точная коррекция смещения частоты

Оценка сдвига частоты на основе корреляции циклического префикса. Это может оценивать сдвиг частоты (abs (foffset)), который меньше половины интервала между поднесущими.

foffset = simpleFreqEstimation(rxWaveform,tmwaveinfo.Info);
fprintf('Cyclic-prefix-based frequency offset = %f Hz\n', foffset)
Cyclic-prefix-based frequency offset = -5.913695 Hz
% Correction
t = 0:numel(txWaveform)-1;
phase = 2*pi*foffset*t'/tmwaveinfo.Info.SamplingRate;
rxWaveform = rxWaveform .* exp(1j*phase);

Более тонкая коррекция

Более точная коррекция на основе DMRS, с 5Hz гранулярностью.

peak = [];
t = 0:numel(rxWaveform)-1;
foffsetRange = -100:5:100;
for foffsetnr = foffsetRange
    % Correction
    phase = 2*pi*foffsetnr*t'/tmwaveinfo.Info.SamplingRate;
    rxWaveform0 = rxWaveform .* exp(1j*phase);
    
    % 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(rxWaveform0,gnb.NRB,gnb.SubcarrierSpacing,NSlot,refGrid);
    peak = [peak; max(mag)];
end
[~,ind] = max(peak);
fineOffset = foffsetRange(ind);
fprintf('DMRS-based frequency offset = %.1f Hz\n', fineOffset)
DMRS-based frequency offset = -25.0 Hz
% Correction
phase = 2*pi*fineOffset*t'/tmwaveinfo.Info.SamplingRate;
rxWaveform = rxWaveform .* exp(1j*phase);

Создайте опорную сетку, охватывающую 10 мсек (один кадр). Эта сетка содержит только DM-RS и используется для синхронизации.

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

Демодуляция OFDM

rxGrid = nrOFDMDemodulate(tmWaveformSync,gnb.NRB,gnb.SubcarrierSpacing,0);

Для всех слотов получите ресурсы PDSCH и DM-RS (местоположения элементов ресурсов). Эта информация присутствует в pdschArray. Затем проводят оценку и выравнивание канала. Полученные данные сохраняются для вычисления 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,...
        'AveragingWindow',[3 3]);
    
    % 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

Отображение диаграммы созвездий.

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

Расчет EVM

Рассчитайте следующие значения EVM (среднеквадратичные и максимальные):

  • EVM на символ OFDM

  • EVM на слот

  • EVM на поднесущую

  • Итого EVM для всего сигнала

Расчет составного EVM по всем каналам. Показанные значения EVM нормализуются силой опорных символов. Эта мощность рассчитывается для рассматриваемого интервала измерения. Например, при вычислении EVM на интервал все опорные символы в этом интервале используются для вычисления мощности, используемой при нормализации. В случае EVM на поднесущую интервал измерения составляет один кадр для FDD и два кадра для TDD. Общее значение EVM измеряется в течение одного кадра для FDD и двух кадров для TDD.

Отображение EVM на символ OFDM, слот и поднесущую.

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

Расчет общего EVM.

[evmRMS, evmMax] = evm(evSqGrid(:),refSqGrid(:));
fprintf("PDSCH RNTIs: %s", mat2str(targetRNTIs));
PDSCH RNTIs: [0 2]
fprintf("RMS EVM = %s", num2str(evmRMS));
RMS EVM = 2.7521
fprintf("Max EVM = %s", num2str(evmMax));
Max EVM = 13.9067

Локальные функции:

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 all(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 = NaN;
else
    evmMax = max(100*sqrt(evSq./mean(refSq)));
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;

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