%hDecodeSubframe Subframe decoding and EVM calculation
%   [EVM,PDSCHSYMEVM] = hDecodeSubframe(SUBFRAMENO,WAVEFORM,GNB,RESAMPLEDWAVEINFO,SIMPARAMS,COUNT)
%   decodes the WAVEFORM with number SUBFRAMENO given parameters in GNB,
%   RESAMPLEDWAVEINFO and SIMPARAMS. It returns the
%   decoded symbols PDSCHSYM and the average EVM. Considers the COUNT of
%   received subframes to avoid decoding the first one due to delay.

%   Copyright 2019-2020 The MathWorks, Inc.

function [EVM,pdschSym] = hDecodeSubframe(subframeNo,waveform,gnb,resampledWaveInfo,simParams,count)
 
    % Get parameters from simParams
    nrtm = simParams.ConfiguredModel.nrtm;
    bw = simParams.ConfiguredModel.bw;
    scs = simParams.ConfiguredModel.scs;
    dm = simParams.ConfiguredModel.duplexing;
    ncellid = simParams.ConfiguredModel.ncellid;
    sv = simParams.ConfiguredModel.sv;
    stopTime = simParams.StopTime;

    % Create generator object for the above NR-TM
    tmwavegen = hNRReferenceWaveformGenerator(nrtm,bw,scs,dm,ncellid,sv);
    [~,tmWaveInfo,resourcesInfo] = generateWaveform(tmwavegen);

    % Generate list of PDSCH config indices where to measure EVM
    pdschArray = hListTargetPDSCHs(tmwavegen.Config,resourcesInfo.WaveformResources);

    % Set general parameters
    symbolsPerSlot = tmWaveInfo.Info.SymbolsPerSlot;
    slotNumber = 1:floor(stopTime/tmWaveInfo.Info.SubframePeriod)*tmWaveInfo.Info.SlotsPerSubframe;
    slotsPerSubframe = tmWaveInfo.Info.SlotsPerSubframe;

    % Subframe number
    sN = mod(double(subframeNo),10); % 10 subframes per frame

    % Calculate the first symbol number for the slot within the subframe
    gnb.NSymbol = 0;

    % Downsample the oversampled waveform 
    waveform = resample(waveform,1,resampledWaveInfo.Info.OSR);

    % OFDM demodulation
    rxGrid = hOFDMDemodulate(gnb,waveform);

    % Initialize the PDSCH symbols array
    pdschSym = [];

    %% Compute EVM
    % The EVM is calculated per:
    % - Slot (time)
    % - OFDM symbol (time)
    % - Subcarrier (frequency)
    % The overall EVM is computed as well.

    % Initialize local and persistent parameters
    noSubcarriers = size(rxGrid,1);
    EVMrmsSlot = zeros(1,tmWaveInfo.Info.SlotsPerSubframe);
    EVMpeakSlot = zeros(1,tmWaveInfo.Info.SlotsPerSubframe);
    EVMrmsSymbol = [];
    EVMpeakSymbol = [];
    evSqSubframe = [];
    refSymSqSubframe = []; 
    persistent evSqSubcarriers;
    persistent refSymSqSubcarriers;
    persistent evSqOverall; 
    persistent refSymSqOverall;

    for n = 0:(tmWaveInfo.Info.SlotsPerSubframe - 1)

        % Slot number
        NSlot = sN*tmWaveInfo.Info.SlotsPerSubframe+n;

        % Retrieve resources
        [pdschIndices,refSymbols,dmrsIndices,dmrsSymbols] = hSlotResources(pdschArray,NSlot);

        % Grid for current slot in subframe
        currentGrid = rxGrid(:,n*symbolsPerSlot+(1:symbolsPerSlot),:);

        % Do not include first two slot symbols for PDSCH EVM (used for control
        % 38.141-1 table 4.9.2.2-2)
        idx = pdschIndices <= 2*noSubcarriers; 
        pdschIndices(idx) = [];
        refSymbols(idx) = [];

        % Channel estimation
        [estChannelGrid,noiseEst] = nrChannelEstimate(currentGrid,dmrsIndices,dmrsSymbols);

        % Get PDSCH resource elements from the received grid
        [pdschRx,pdschHest] = nrExtractResources(pdschIndices,currentGrid,estChannelGrid);

        % Equalization
        pdschEq = nrEqualizeMMSE(pdschRx,pdschHest,noiseEst);
        pdschSym = [pdschSym; pdschEq(:)]; %#ok<AGROW>

        % RMS and peak EVM calculation per slot
        if ~any(isnan(pdschEq)) && ~isempty(pdschEq)
            evSq = abs(refSymbols-pdschEq).^2; % Norm square of elements in error vector
            refSymSq = abs(refSymbols).^2; % Norm square of reference symbols
            EVMrmsSlot(n+1) = 100*sqrt(sum(evSq)/sum(refSymSq)); % RMS EVM
            EVMinst = 100*sqrt(evSq/mean(refSymSq)); % Instantaneous EVM
            EVMpeakSlot(n+1) = max(EVMinst); % Peak EVM 
        else
            evSq = []; % Norm square of elements in error vector
            refSymSq = []; % Norm square of reference symbols
            EVMrmsSlot(n+1) = NaN; % RMS EVM
            EVMpeakSlot(n+1) = NaN; % Peak EVM 
        end

        % Store results for overall EVM calculation
        evSqSubframe = [evSqSubframe evSq']; %#ok<AGROW>
        refSymSqSubframe = [refSymSqSubframe refSymSq'];  %#ok<AGROW>

        % RMS and peak EVM calculation per OFDM symbol
        gridEvSq = zeros(size(currentGrid));
        gridEvSq(pdschIndices) = evSq;
        gridRefSymSq = zeros(size(currentGrid));
        gridRefSymSq(pdschIndices) = refSymSq;
        EVMrmsSym = zeros(1,symbolsPerSlot);
        EVMpeakSym = zeros(1,symbolsPerSlot);
        for symbolNo = 1:symbolsPerSlot
            evSqSymbol = nonzeros(gridEvSq(:,symbolNo))';
            refSymSqSymbol = nonzeros(gridRefSymSq(:,symbolNo))'; 
            if ~any(isnan(evSqSymbol)) && ~isempty(evSqSymbol)
                EVMrmsSym(symbolNo) = 100*sqrt(sum(evSqSymbol)/sum(refSymSqSymbol));
                EVMpeakSym(symbolNo) = max(100*sqrt(evSqSymbol/mean(refSymSqSymbol)));
            else
                EVMrmsSym(symbolNo) = NaN;
                EVMpeakSym(symbolNo) = NaN;
            end
        end
        EVMrmsSymbol = [EVMrmsSymbol EVMrmsSym]; %#ok<AGROW>
        EVMpeakSymbol = [EVMpeakSymbol EVMpeakSym]; %#ok<AGROW>

        % Store results for EVM calculation per subcarrier
        for subcarrierNo = 1:noSubcarriers
            evSqSubcarrier{subcarrierNo,:} = nonzeros(gridEvSq(subcarrierNo,:))'; %#ok<AGROW>
            refSymSqSubcarrier{subcarrierNo,:} = nonzeros(gridRefSymSq(subcarrierNo,:))';  %#ok<AGROW>
        end
        evSqSubcarriers = [evSqSubcarriers evSqSubcarrier]; %#ok<AGROW>
        refSymSqSubcarriers = [refSymSqSubcarriers refSymSqSubcarrier]; %#ok<AGROW>
    end

    % Overall RMS EVM
    evSqOverall = [evSqOverall evSqSubframe];
    refSymSqOverall = [refSymSqOverall refSymSqSubframe];
    EVM = 100*sqrt(sum(evSqOverall)/sum(refSymSqOverall));

    % RMS and peak EVM calculation per subcarrier
    for subcarrierNo = 1:noSubcarriers
        evSqSub = [];
        refSymSqSub = [];
        for slotNo = 1:count*tmWaveInfo.Info.SlotsPerSubframe
            evSqSub = [evSqSub evSqSubcarriers{subcarrierNo,slotNo}]; %#ok<AGROW>
            refSymSqSub = [refSymSqSub refSymSqSubcarriers{subcarrierNo,slotNo}]; %#ok<AGROW>
        end 
        if ~any(isnan(evSqSub)) && ~isempty(evSqSub)
            EVMrmsSubcarrier(subcarrierNo) = 100*sqrt(sum(evSqSub)/sum(refSymSqSub)); %#ok<AGROW>
            EVMpeakSubcarrier(subcarrierNo) = max(100*sqrt(evSqSub/mean(refSymSqSub))); %#ok<AGROW>
        else
            EVMrmsSubcarrier(subcarrierNo) = NaN; %#ok<AGROW>
            EVMpeakSubcarrier(subcarrierNo) = NaN; %#ok<AGROW>
        end
    end

    %% Plot EVM graphs for all the received subframes

    % Create a new figure when starting simulation
    if count == 1
        figure
        figEVM = gcf; %#ok<NASGU>
    end

    % Plot EVM per slot
    subplot(3,1,1);
    % Set 10 different colors, one per subframe
    colors = [0, 0.4470, 0.7410; 0.8500, 0.3250, 0.0980; 0.9290, 0.6940, 0.1250; 0.4940, ...
        0.1840, 0.5560; 0.4660, 0.6740, 0.1880; 0.3010, 0.7450, 0.9330; 0.6350, 0.0780, ...
        0.1840; 0.9290, 0.9, 0.3; 0.9290, 0.5, 0.9; 0.6660, 0.3740, 0.2880;]; 
    bar(slotNumber(((slotsPerSubframe*(count-1))+1):((slotsPerSubframe*(count-1))+1) ...
        +slotsPerSubframe-1)-1,EVMpeakSlot,'FaceColor',colors(subframeNo+1,:),'FaceAlpha',0.3,'EdgeColor','none');
    hold on;
    bar(slotNumber(((slotsPerSubframe*(count-1))+1):((slotsPerSubframe*(count-1))+1) ...
        +slotsPerSubframe-1)-1,EVMrmsSlot,'FaceColor',colors(subframeNo+1,:),'EdgeColor','none');
    grid on;
    if count > 10 % Maximum number of subframes
        xlim([slotNumber(slotsPerSubframe*(count-10)) slotNumber((slotsPerSubframe*count)-1)]);
    end
    xticks('auto')
    title('EVM per Slot');
    xlabel('Slot Number');
    ylabel('EVM (%)');
    if count <= 10 % Maximum number of subframes
        legendCount = count;
    else
        legendCount = 10;
    end
    subframeLegend = "";
    for i=1:legendCount
        subframeLegend(i*2-1) = "Peak - Subframe " + (i-1); 
        subframeLegend(i*2) = "RMS - Subframe " + (i-1);
    end
    legend(subframeLegend,'Location','EastOutside','NumColumns',2);

    % Plot EVM per OFDM symbol
    subplot(3,1,2);
    firstSymbol = (((slotsPerSubframe*(count-1))+1)-1)*symbolsPerSlot;
    lastSymbol = ((((slotsPerSubframe*(count-1))+1)+slotsPerSubframe)-1)*symbolsPerSlot-1;
    bar(firstSymbol:lastSymbol,EVMpeakSymbol,'FaceColor',colors(subframeNo+1,:),'FaceAlpha',0.3,'EdgeColor','none');
    hold on;
    bar(firstSymbol:lastSymbol,EVMrmsSymbol,'FaceColor',colors(subframeNo+1,:),'EdgeColor','none');
    grid on;
    if count > 10 % Maximum number of subframes
        xlim([slotNumber(slotsPerSubframe*(count-10))*symbolsPerSlot slotNumber((slotsPerSubframe*count)-1)*symbolsPerSlot]);
    end
    xlabel('OFDM Symbol Number');
    ylabel('EVM (%)'); 
    title('EVM per OFDM Symbol');

    % Plot EVM per subcarrier
    subplot(3,1,3);
    bar(1:noSubcarriers,EVMpeakSubcarrier,'FaceColor',[0.1660, 0.3740, 0.0080],'FaceAlpha',0.3,'EdgeColor','none')
    hold on;
    bar(1:noSubcarriers,EVMrmsSubcarrier,'FaceColor',[0.1660, 0.3740, 0.0080],'EdgeColor','none')
    grid on;
    xlabel('Subcarrier Number');
    ylabel('EVM (%)');
    legend('Peak','RMS','Location','EastOutside');
    title('EVM per Subcarrier');
    xlim([0 noSubcarriers+1]);
    hold off;
    set(gcf,'units','inches','position',[6,2,11,7.5]);

    %% Save overall EVM to workspace

    assignin('base','EVMrms',EVM);

end