Сгенерируйте и передайте формы волны 5G Используя тестовое оборудование и измерительное оборудование

Этот пример использует:

В этом примере показано, как сгенерировать и передать беспроводную форму волны 5G с помощью 5G Toolbox™, Instrument Control Toolbox™ и инструментов Keysight Technologies® RF.

Введение

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

Требования

Чтобы запустить этот пример, вам нужно:

Генератор сигнала Keysight Technologies® E4438C
Сигнал Keysight Technologies® N9030A анализатор
5G Toolbox™
Instrument Control Toolbox™
Signal Processing Toolbox™
Пакет поддержки Instrument Control Toolbox™ для библиотек Keysight™ IO и интерфейса VISA

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

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

Для получения дополнительной информации о Приложении Wireless Waveform Generator обратитесь к 5G Waveform Generator.

На панели инструментов MATLAB Apps выберите Wireless Waveform Generator App. Откройте приложение Wireless Waveform Generator и сконфигурируйте, чтобы сгенерировать 5G NR-TM (Тестовая модель) форма волны. Использование этой функции приложения требует 5G Toolbox™.

В разделе Waveform Type выберите формы волны Test Models (NR TM). На левой панели приложения, на вкладке Waveform, можно установить параметры выбранной формы волны. Для этого примера:

Установите Частотный диапазон как FR1 (GHz на 450 МГц - 6).
Установите Тестовую модель как полосу Full, универсальная форма 64QAM представленный как NR-FR1-TM3.1.
Установите пропускную способность Канала как 10 МГц, Поднесущая, располагающая с интервалами как 30 кГц и Дуплексный режим как FDD.
Нажмите на Generate.

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

Передайте беспроводной сигнал

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

Читайте синфазный и квадратура (IQ) данные из Signal Analyzer по TCP/IP

Для анализа беспроводной передачи в MATLAB Instrument Control Toolbox используется, чтобы сконфигурировать 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");

Инициируйте инструмент, чтобы сделать измерение. Ожидайте операции измерения, которая будет завершена и считана в форме волны. Прежде, чем обработать данные, разделите меня & компоненты 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

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

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;

Представление спектра может иметь больше информации, чем представление области времени данных. Например, можно использовать представление спектра, чтобы идентифицировать основные диапазоны частот, пропускную способность сигнала, и т.д. У вас должен быть Signal Processing Toolbox, чтобы построить представление спектра.

% 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 передача формы волны NR-TM

Измерьте EVM 5G форма волны NR-TM

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

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

Каждая форма волны NR-TM задана комбинацией:

Имя ТМ NR
Пропускная способность канала
Интервал поднесущей
Режим Duplexing

% 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 выбирает целевой PDSCHs, чтобы анализировать. Это основано на RNTI. По умолчанию следующий RNTI рассматривается для вычисления EVM:

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

Эти значения могут быть заменены. Функциональный hListTargetPDSCHs также принимает, что третий параметр с набором RNTIs рассматривает для измерения 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, ищущий смещения с шагом 10 кГц до 200 кГц

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, с гранулярностью на 5 Гц.

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 (RMS и максимум) значения:

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

Документация