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

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

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

Введение

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

Требования

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

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

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

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

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

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

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

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

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

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

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

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

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

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

% 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.
Отрегулируйте вход buffer size так, чтобы он мог хранить данные, возвраты прибором.
Установите тайм-аут, чтобы обеспечить достаточное время для измерения и передачи данных.
Подключите к прибору.

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;

Представление спектра может иметь больше информации, чем представление данных во временном интервале. Для примера можно использовать представление спектра, чтобы идентифицировать основные полосы частот, полосу пропускания сигнала и т. Д. Вы должны иметь 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 в приложении Waveform Generator, чтобы остановить 5G передачу сигнала NR-TM

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

После того, как вы получите захваченные данные baseband 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 - два. The tmwaveinfo выход содержит информацию о сгенерированном сигнале. The 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 (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

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