nFrames = 5; % Number of 10 ms frames
SNRdB = 10;  % SNR in dB

carrier = nrCarrierConfig;
carrier.SubcarrierSpacing = 15;
carrier.NSizeGrid = 52;

NStartBWP = 0;
NSizeBWP = 52;

csirs = nrCSIRSConfig;
csirs.CSIRSType = {'nzp','nzp','nzp'};
csirs.RowNumber = [4 4 4];
csirs.NumRB = 52;
csirs.RBOffset = 0;
csirs.CSIRSPeriod = [4 0];
csirs.SymbolLocations = {0, 0, 0};
csirs.SubcarrierLocations = {0, 4, 8};
csirs.Density = {'one','one','one'};

nTxAnts = 4;
nRxAnts = 4;

validateCSIRSConfig(carrier,csirs,nTxAnts);

% Configure CSI reporting configuration parameters
reportConfig.NStartBWP = NStartBWP;
reportConfig.NSizeBWP = NSizeBWP;
reportConfig.CodebookType = 'Type1SinglePanel';
reportConfig.PanelDimensions = [2 1];
reportConfig.CQIMode = 'Subband';
reportConfig.PMIMode = 'Subband';
reportConfig.SubbandSize = 4;
reportConfig.PRGSize = [];
reportConfig.CodebookMode = 1;
reportConfig.CodebookSubsetRestriction = [];
reportConfig.i2Restriction = [];
reportConfig.RIRestriction = [];

% Create channel object
channel = nrTDLChannel;
channel.NumTransmitAntennas = nTxAnts;
channel.NumReceiveAntennas = nRxAnts;
channel.DelayProfile = 'TDL-C';
channel.MaximumDopplerShift = 50;
channel.DelaySpread = 300e-9;

waveformInfo = nrOFDMInfo(carrier);
channel.SampleRate = waveformInfo.SampleRate;

chInfo = info(channel);
maxChDelay = ceil(max(chInfo.PathDelays*channel.SampleRate)) + chInfo.ChannelFilterDelay;

% Calculate total number of slots
totSlots = nFrames*carrier.SlotsPerFrame;

% Initialize variables to store CQI, PMI, RI, and subband SINR values for
% practical and perfect channel estimation scenarios
cqiPracticalPerSlot = [];
subbandCQIPractical = [];
pmiPracticalPerSlot = struct('i1',[],'i2',[]);
SINRPerSubbandPerCWPractical = [];
cqiPerfectPerSlot = [];
subbandCQIPerfect = [];
pmiPerfectPerSlot = struct('i1',[],'i2',[]);
SINRPerSubbandPerCWPerfect = [];
riPracticalPerSlot = [];
riPerfectPerSlot = [];

% Get number of CSI-RS ports
csirsPorts = csirs.NumCSIRSPorts(1);

% Get CDM lengths corresponding to configured CSI-RS resources
cdmLengths = getCDMLengths(csirs);

% Initialize the practical timing offset as zero. It is updated in the
% slots when the correlation is strong for practical synchronization.
offsetPractical = 0;

% Initialize a variable to store the slot numbers (0-based) where CSI-RS
% resources are present
activeSlots = [];

% Loop over all slots
for slotIdx = 0:totSlots - 1
    % Create carrier resource grid for one slot
    csirsSlotGrid = nrResourceGrid(carrier,csirsPorts);

    % Update slot number in carrier configuration object
    carrier.NSlot = slotIdx;

    % Generate CSI-RS indices and symbols
    csirsInd = nrCSIRSIndices(carrier,csirs);
    csirsSym = nrCSIRS(carrier,csirs);

    % Map CSI-RS to slot grid
    csirsSlotGrid(csirsInd) = csirsSym;

    % Map CSI-RS ports to transmit antennas
    wtx = eye(csirsPorts,nTxAnts);
    txGrid = reshape(reshape(csirsSlotGrid,[],csirsPorts)*wtx,size(csirsSlotGrid,1),size(csirsSlotGrid,2),nTxAnts);

    % Perform OFDM modulation to generate time-domain waveform
    txWaveform = nrOFDMModulate(carrier,txGrid);

    % Append zeros at the end of the transmitted waveform to flush channel
    % content. These zeros take into account any delay introduced in the
    % channel. The channel delay is a mix of multipath delay and
    % implementation delay. This value may change depending on the sampling
    % rate, delay profile, and delay spread.
    txWaveform = [txWaveform; zeros(maxChDelay,size(txWaveform,2))]; %#ok<AGROW>

    % Transmit waveform through channel
    [rxWaveform,pathGains,sampleTimes] = channel(txWaveform);

    % Generate and add AWGN to received waveform
    SNR = 10^(SNRdB/20);                                                % Linear SNR value
    N0 = 1/(sqrt(2.0*nRxAnts*double(waveformInfo.Nfft))*SNR);           % Noise standard deviation
    rng('default');                                                     % Set RNG state for repeatability
    noise = N0*complex(randn(size(rxWaveform)),randn(size(rxWaveform)));
    rxWaveform = rxWaveform + noise;

    % Perform practical timing estimation. Correlate the received waveform
    % with the CSI-RS to obtain the timing offset estimate and the
    % correlation magnitude. Use the function hSkipWeakTimingOffset to
    % update the receiver timing offset. If the correlation peak is weak,
    % the current timing estimate is ignored and the previous offset is
    % used.
    [t,mag] = nrTimingEstimate(carrier,rxWaveform,csirsInd,csirsSym);
    offsetPractical = hSkipWeakTimingOffset(offsetPractical,t,mag);

    % Get path filters
    pathFilters = getPathFilters(channel);
    % Perform perfect timing estimation
    offsetPerfect = nrPerfectTimingEstimate(pathGains,pathFilters);

    % Perform time-domain offset correction for practical and
    % perfect timing estimation scenarios
    rxWaveformPractical = rxWaveform(1+offsetPractical:end,:);
    rxWaveformPerfect = rxWaveform(1+offsetPerfect:end,:);

    % Perform OFDM demodulation on previously synchronized waveforms
    rxGridPractical = nrOFDMDemodulate(carrier,rxWaveformPractical);
    rxGridPerfect = nrOFDMDemodulate(carrier,rxWaveformPerfect);

    % Append zeros in the event where the timing synchronization results in
    % an incomplete slot
    symbPerSlot = carrier.SymbolsPerSlot;
    K = size(rxGridPractical,1);
    LPractical = size(rxGridPractical,2);
    LPerfect = size(rxGridPerfect,2);
    if LPractical < symbPerSlot
        rxGridPractical = cat(2,rxGridPractical,zeros(K,symbPerSlot-LPractical,nRxAnts));
    end
    if LPerfect < symbPerSlot
        rxGridPerfect = cat(2,rxGridPerfect,zeros(K,symbPerSlot-LPerfect,nRxAnts));
    end
    rxGridPractical = rxGridPractical(:,1:symbPerSlot,:);
    rxGridPerfect = rxGridPerfect(:,1:symbPerSlot,:);

    % Consider only the NZP-CSI-RS symbols and indices for channel estimation
    nzpCSIRSSym = csirsSym(csirsSym ~= 0);
    nzpCSIRSInd = csirsInd(csirsSym ~= 0);

    % Calculate practical channel estimate. Use a time-averaging window
    % that covers all of the transmitted CSI-RS symbols.
    [PracticalHest,nVarPractical] = nrChannelEstimate(carrier,rxGridPractical, ...
        nzpCSIRSInd,nzpCSIRSSym,'CDMLengths',cdmLengths,'AveragingWindow',[0 5]);

    % Perform perfect channel estimation
    PerfectHest = nrPerfectChannelEstimate(carrier,pathGains,pathFilters,offsetPerfect,sampleTimes);

    % Get perfect noise estimate value from noise realization
    noiseGrid = nrOFDMDemodulate(carrier,noise(1+offsetPerfect:end,:));
    nVarPerfect = var(noiseGrid(:));
    if ~isempty(nzpCSIRSInd)
        activeSlots = [activeSlots slotIdx]; %#ok<AGROW>

        % Calculate the RI value using practical channel estimate
        numLayersPractical = hRISelect(carrier,csirs,reportConfig,PracticalHest,nVarPractical);

        % Calculate CQI and PMI values using practical channel estimate
        [cqiPractical,pmiPractical,cqiInfoPractical] = hCQISelect(carrier,csirs, ...
            reportConfig,numLayersPractical,PracticalHest,nVarPractical);
        numCodeWordsPr = size(cqiPractical,2);
        numSBs = size(cqiPractical,1);

        % Store CQI, PMI, RI, and subband SINR values of each slot for the
        % practical channel estimation scenario. Because the number of
        % codewords can vary based on the rank, append NaNs to the
        % CQI-related variables to account for the second codeword
        % information in the slots where only one codeword is present.
        riPracticalPerSlot(1,slotIdx+1) = numLayersPractical; %#ok<SAGROW> 
        cqiPracticalPerSlot(:,:,slotIdx+1) = [cqiPractical NaN(numSBs,2-numCodeWordsPr)]; %#ok<SAGROW> 
        pmiPracticalPerSlot(slotIdx+1) = pmiPractical;
        subbandCQIPractical(:,:,slotIdx+1) = [cqiInfoPractical.SubbandCQI NaN(numSBs,2-numCodeWordsPr)]; %#ok<SAGROW> 
        SINRPerSubbandPerCWPractical(:,:,slotIdx+1) = [cqiInfoPractical.SINRPerSubbandPerCW NaN(numSBs,2-numCodeWordsPr)]; %#ok<SAGROW> 

        % Calculate the RI value using perfect channel estimate
        numLayersPerfect = hRISelect(carrier,csirs,reportConfig,PerfectHest,nVarPerfect);

        % Calculate CQI and PMI values using perfect channel estimate
        [cqiPerfect,pmiPerfect,cqiInfoPerfect] = hCQISelect(carrier,csirs, ...
            reportConfig,numLayersPerfect,PerfectHest,nVarPerfect);
        numCodeWordsPe = size(cqiPerfect,2);

        % Store CQI, PMI, RI, and subband SINR values of each slot for the
        % perfect channel estimation scenario. Because the number of
        % codewords can vary based on the rank, append NaNs to the
        % CQI-related variables to account for the second codeword
        % information in the slots where only one codeword is present.
        riPerfectPerSlot(1,slotIdx+1) = numLayersPerfect; %#ok<SAGROW> 
        cqiPerfectPerSlot(:,:,slotIdx+1) = [cqiPerfect NaN(numSBs,2-numCodeWordsPe)]; %#ok<SAGROW> 
        subbandCQIPerfect(:,:,slotIdx+1) = [cqiInfoPerfect.SubbandCQI NaN(numSBs,2-numCodeWordsPe)]; %#ok<SAGROW> 
        pmiPerfectPerSlot(slotIdx+1) = pmiPerfect;
        SINRPerSubbandPerCWPerfect(:,:,slotIdx+1) = [cqiInfoPerfect.SINRPerSubbandPerCW NaN(numSBs,2-numCodeWordsPe)]; %#ok<SAGROW> 
    end
end

plotWidebandCQIAndSINR(cqiPracticalPerSlot,cqiPerfectPerSlot,SINRPerSubbandPerCWPractical,SINRPerSubbandPerCWPerfect,activeSlots);

% Provide slot number for which subband CQI values and PMI i2 indices must be plotted
slotNum = 0; % Consider slot number as 0 (0-based) here, because CSI is
             % reported in that slot for the configured CSI-RS resource(s)
plotSubbandCQIAndSINR(subbandCQIPractical,subbandCQIPerfect,SINRPerSubbandPerCWPractical,SINRPerSubbandPerCWPerfect,activeSlots,slotNum);

% Set the showPMIandRIPlots flag to true, if you want to see the plots for
% the PMI and RI 
showPMIandRIPlots = true;
if showPMIandRIPlots
    plotPMIAndRI(pmiPracticalPerSlot,pmiPerfectPerSlot,riPracticalPerSlot,riPerfectPerSlot,activeSlots,slotNum);
end

function validateCSIRSConfig(carrier,csirs,nTxAnts)
%   validateCSIRSConfig(CARRIER,CSIRS,NTXANTS) validates the CSI-RS
%   configuration, given the carrier specific configuration object CARRIER,
%   CSI-RS configuration object CSIRS, and the number of transmit antennas
%   NTXANTS.

    % Validate the number of CSI-RS ports
    if ~isscalar(unique(csirs.NumCSIRSPorts))
        error('nr5g:InvalidCSIRSPorts',...
            'All the CSI-RS resources must be configured to have the same number of CSI-RS ports.');
    end

    % Validate the CDM lengths
    if ~iscell(csirs.CDMType)
        cdmType = {csirs.CDMType};
    else
        cdmType = csirs.CDMType;
    end
    if (~all(strcmpi(cdmType,cdmType{1})))
        error('nr5g:InvalidCSIRSCDMTypes',...
            'All the CSI-RS resources must be configured to have the same CDM lengths.');
    end
    if nTxAnts < csirs.NumCSIRSPorts(1)
        error('nr5g:InvalidNumTxAnts',['Number of transmit antennas (' num2str(nTxAnts)...
            ') must be greater than or equal to the number of CSI-RS ports (' num2str(csirs.NumCSIRSPorts(1)) ').']);
    end

    % Check for the overlap between the CSI-RS indices
    csirsInd = nrCSIRSIndices(carrier,csirs,"OutputResourceFormat",'cell');
    numRes = numel(csirsInd);
    csirsIndAll = cell(1,numRes);
    ratioVal = csirs.NumCSIRSPorts(1)/prod(getCDMLengths(csirs));
    for resIdx = 1:numRes
        if ~isempty(csirsInd{resIdx})
            grid = nrResourceGrid(carrier,csirs.NumCSIRSPorts(1));
            [~,tempInd] = nrExtractResources(csirsInd{resIdx},grid);
            if numel(tempInd)/numel(csirsInd{resIdx}) ~= ratioVal
                error('nr5g:OverlappedCSIRSREsSingleResource',['CSI-RS indices of resource '...
                    num2str(resIdx) ' must be unique. Try changing the symbol or subcarrier locations.']);
            end
            csirsIndAll{resIdx} = tempInd(:);
            for idx = 1:resIdx-1
                overlappedInd = ismember(csirsIndAll{idx},csirsIndAll{resIdx});
                if any(overlappedInd)
                    error('nr5g:OverlappedCSIRSREsMultipleResources',['The resource elements of the '...
                        'configured CSI-RS resources must not overlap. Try changing the symbol or '...
                        'subcarrier locations of CSI-RS resource ' num2str(idx) ' and resource ' num2str(resIdx) '.']);
                end
            end
        end
    end
end

function cdmLengths = getCDMLengths(csirs)
%   CDMLENGTHS = getCDMLengths(CSIRS) returns the CDM lengths given the
%   CSI-RS configuration object CSIRS.

    CDMType = csirs.CDMType;
    if ~iscell(csirs.CDMType)
        CDMType = {csirs.CDMType};
    end
    CDMTypeOpts = {'noCDM','fd-CDM2','CDM4','CDM8'};
    CDMLengthOpts = {[1 1],[2 1],[2 2],[2 4]};
    cdmLengths = CDMLengthOpts{strcmpi(CDMTypeOpts,CDMType{1})};
end

function plotWidebandCQIAndSINR(cqiPracticalPerSlot,cqiPerfectPerSlot,SINRPerSubbandPerCWPractical,SINRPerSubbandPerCWPerfect,activeSlots)
%   plotWidebandCQIAndSINR(CQIPRACTICALPERSLOT,CQIPERFECTPERSLOT,SINRPERSUBBANDPERCWPRACTICAL,SINRPERSUBBANDPERCWPERFECT,ACTIVESLOTS)
%   plots the wideband SINR and wideband CQI values for each codeword
%   across all specified active slots ACTIVESLOTS (0-based) (in which the
%   CQI is reported as other than NaN) for practical and perfect channel
%   estimation cases.

    % Check if there are no slots in which CSI-RS is present
    if isempty(activeSlots)
        disp('No CQI data to plot, because there are no slots in which CSI-RS is present.');
        return;
    end

    % Get slot numbers (1-based) in which CQI is reported
    activeSlotNum = activeSlots + 1;

    if all(isnan(cqiPracticalPerSlot(1,2,activeSlotNum)))
        numCW = 1;
    else
        numCW = 2;
    end
    cqiPracticalPerCW = permute(cqiPracticalPerSlot(1,:,:),[1 3 2]);
    cqiPerfectPerCW = permute(cqiPerfectPerSlot(1,:,:),[1 3 2]);
    SINRPerCWPractical = permute(SINRPerSubbandPerCWPractical(1,:,:),[1 3 2]);
    SINRPerCWPerfect = permute(SINRPerSubbandPerCWPerfect(1,:,:),[1 3 2]);

    % Extract wideband CQI indices for slots where CSI-RS is present
    cqiPracticalPerCWActiveSlots = cqiPracticalPerCW(1,activeSlotNum,:);
    cqiPerfectPerCWActiveSlots = cqiPerfectPerCW(1,activeSlotNum,:);
    widebandSINRPractical = 10*log10(SINRPerCWPractical(1,activeSlotNum,:));
    widebandSINRPerfect = 10*log10(SINRPerCWPerfect(1,activeSlotNum,:));
    xTickLables = num2cell(activeSlots(:));

    if isempty(reshape(cqiPracticalPerCWActiveSlots(:,:,1),1,[]))
        disp('No CQI data to plot, because all CQI values are NaNs.');
        return;
    end

    figure();
    subplot(211)
    plot(widebandSINRPerfect(:,:,1),'ro-')
    hold on;
    plot(widebandSINRPractical(:,:,1),'b*-')
    if numCW == 2
        hold on;
        plot(widebandSINRPerfect(:,:,2),'rs:')
        hold on;
        plot(widebandSINRPractical(:,:,2),'bd:')
        title('Estimated Wideband SINR Values for Codeword 1&2');
        xlabel('Slots');
        ylabel('Wideband SINR Values in dB');
        xticks(1:size(cqiPracticalPerCWActiveSlots,2));
        xticklabels(xTickLables);
        [lowerBound,upperBound] = bounds([widebandSINRPractical(:);widebandSINRPerfect(:)],'omitnan');
        ylim([lowerBound-1 upperBound+3.5]);
        legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.','Codeword 2:Perfect channel est.','Codeword 2:Practical channel est.'});
    else
        title('Estimated Wideband SINR Values for Codeword 1');
        xlabel('Slots');
        ylabel('Wideband SINR Values in dB');
        xticks(1:size(cqiPracticalPerCWActiveSlots,2));
        xticklabels(xTickLables);
        [lowerBound,upperBound] = bounds([widebandSINRPractical(:);widebandSINRPerfect(:)]);
        ylim([lowerBound-1 upperBound+3.5]);
        legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.'});
    end

    subplot(212)
    plot(cqiPerfectPerCWActiveSlots(:,:,1),'r-o');
    hold on;
    plot(cqiPracticalPerCWActiveSlots(:,:,1),'b-*');
    if numCW == 2
        hold on;
        plot(cqiPerfectPerCWActiveSlots(:,:,2),'r:s');
        hold on;
        plot(cqiPracticalPerCWActiveSlots(:,:,2),'b:d');
        xlabel('Slots');
        xticks(1:size(cqiPracticalPerCWActiveSlots,2));
        xticklabels(xTickLables);
        ylabel('Wideband CQI');
        [lowerBound,upperBound] = bounds([cqiPracticalPerCWActiveSlots(:);cqiPerfectPerCWActiveSlots(:)]);
        ylim([lowerBound-1 upperBound+3.5]);
        title('Wideband CQI Values for Codeword 1&2');
        legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.','Codeword 2:Perfect channel est.','Codeword 2:Practical channel est.'});
    else
        xlabel('Slots');
        xticks(1:size(cqiPracticalPerCWActiveSlots,2));
        xticklabels(xTickLables);
        ylabel('Wideband CQI');
        [lowerBound,upperBound] = bounds([cqiPracticalPerCWActiveSlots(:,:,1)';cqiPerfectPerCWActiveSlots(:,:,1)']);
        ylim([lowerBound-1 upperBound+3.5]);
        title('Wideband CQI Values for Codeword 1');
        legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.'});
    end
end

function plotSubbandCQIAndSINR(subbandCQIPractical,subbandCQIPerfect,SINRPerCWPractical,SINRPerCWPerfect,activeSlots,nslot)
%   plotSUBBANDCQIAndSINR(SUBBANDCQIPRACTICAL,SUBBANDCQIPERFECT,SINRPERCWPRACTICAL,SINRPERCWPERFECT,ACTIVESLOTS,NSLOT)
%   plots the SINR and CQI values for each codeword across all of the
%   subbands for practical and perfect channel estimation cases for the
%   given slot NSLOT (0-based) among all specified active slots ACTIVESLOTS
%   (0-based). The function does not plot the values if CQIMode is
%   'Wideband' or if the CQI and SINR values are all NaNs in the given
%   NSLOT.

    % Check if there are no slots in which CSI-RS is present
    if isempty(activeSlots)
        disp('No CQI data to plot, because there are no slots in which CSI-RS is present.');
        return;
    end

    if all(isnan(subbandCQIPractical(:,2,activeSlots+1)))
        numCW = 1;
    else
        numCW = 2;
    end
    numSubbands = size(subbandCQIPractical,1);
    if numSubbands > 1
        Mode = 'Subband';
    else
        Mode = 'Wideband';
    end

    % Plot subband CQI values
    if strcmpi(Mode,'subband')
        subbandCQIPerCWPractical = subbandCQIPractical(2:end,:,nslot+1);
        subbandCQIPerCWPerfect = subbandCQIPerfect(2:end,:,nslot+1);
        subbandSINRPerCWPractical = 10*log10(SINRPerCWPractical(2:end,:,nslot+1));
        subbandSINRPerCWPerfect = 10*log10(SINRPerCWPerfect(2:end,:,nslot+1));

        xTickLables = num2cell(1:size(subbandCQIPerfect,1)-1); % Excluding the value of wideband
        if ~any(nslot == activeSlots)
            disp(['For the specified slot (' num2str(nslot) '), CQI values are not reported. Please choose another slot number.'])
        else
            figure();
            subplot(211)
            plot(subbandSINRPerCWPerfect(:,1),'ro-')
            hold on;
            plot(subbandSINRPerCWPractical(:,1),'b*-')
            if numCW == 2
                hold on;
                plot(subbandSINRPerCWPerfect(:,2),'rs:')
                hold on;
                plot(subbandSINRPerCWPractical(:,2),'bd:')
                title(['Estimated Subband SINR Values for Codeword 1&2 in Slot ' num2str(nslot)]);
                xlabel('Subbands');
                ylabel('Subband SINR Values in dB');
                xticks(1:numSubbands);
                xticklabels(xTickLables);
                xlim([0 numSubbands+1]);
                [lowerBound,upperBound] = bounds([subbandSINRPerCWPerfect(:);subbandSINRPerCWPractical(:)],'omitnan');
                ylim([lowerBound-1 upperBound+3.5]);
                legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.','Codeword 2:Perfect channel est.','Codeword 2:Practical channel est.'});
            else
                title(['Estimated Subband SINR Values for Codeword 1 in Slot ' num2str(nslot)]);
                xlabel('Subbands');
                ylabel('Subband SINR Values in dB');
                xticks(1:numSubbands);
                xticklabels(xTickLables);
                xlim([0 numSubbands+1]);
                [lowerBound,upperBound] = bounds([subbandSINRPerCWPerfect(:,1);subbandSINRPerCWPractical(:,1)]);
                ylim([lowerBound-1 upperBound+3.5]);
                legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.'});
            end

            subplot(212)
            plot(subbandCQIPerCWPerfect(:,1),'r-o');
            hold on;
            plot(subbandCQIPerCWPractical(:,1),'b-*');
            if numCW == 2
                hold on;
                plot(subbandCQIPerCWPerfect(:,2),'r:s');
                hold on;
                plot(subbandCQIPerCWPractical(:,2),'b:d');
                title(['Subband CQI Values for Codeword 1&2 in Slot ' num2str(nslot)]);
                xlabel('Subbands');
                ylabel('Subbands CQI');
                xticks(1:numSubbands);
                xticklabels(xTickLables);
                xlim([0 numSubbands+1]);
                [lowerBound,upperBound] = bounds([subbandCQIPerCWPerfect(:);subbandCQIPerCWPractical(:)],'omitnan');
                ylim([lowerBound-1 upperBound+3.5]);
                legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.','Codeword 2:Perfect channel est.','Codeword 2:Practical channel est.'});
            else
                title(['Subband CQI Values for Codeword 1 in Slot ' num2str(nslot)]);
                xlabel('Subbands');
                ylabel('Subbands CQI');
                xticks(1:numSubbands);
                xticklabels(xTickLables);
                xlim([0 numSubbands+1]);
                [lowerBound,upperBound] = bounds([subbandCQIPerCWPerfect(:,1);subbandCQIPerCWPractical(:,1)]);
                ylim([lowerBound-1 upperBound+3.5]);
                legend({'Codeword 1:Perfect channel est.','Codeword 1:Practical channel est.'});
            end
        end
    end
end

function plotPMIAndRI(pmiPracticalPerSlot,pmiPerfectPerSlot,riPracticalPerSlot,riPerfectPerSlot,activeSlots,nslot)
%   plotPMIAndRI(PMIPRACTICALPERSLOT,PMIPERFECTPERSLOT,RIPRACTICALPERSLOT,RIPERFECTPERSLOT,ACTIVESLOTS,NSLOT)
%   plots the RI and PMI i1 indices across all specified active slots
%   ACTIVESLOTS (0-based), for practical and perfect channel estimation
%   scenarios. The function also plots the i2 indices of practical and
%   perfect channel estimation scenarios across all specified active slots
%   when the PMI mode is 'Wideband' or plots i2 indices across all of the
%   subbands for the specified slot number NSLOT (0-based) when the PMI
%   mode is 'Subband'.

    % Check if there are no slots in which CSI-RS is present
    if isempty(activeSlots)
        disp('No PMI and RI data to plot, because there are no slots in which CSI-RS is present.');
        return;
    end

    % Get slot numbers (1-based) in which PMI is reported
    activeSlotNum = activeSlots+1;
    
    % Get number of active slots
    numActiveSlots = numel(activeSlotNum);
    
    numi1Indices = numel(pmiPracticalPerSlot(activeSlotNum(1)).i1);
    
    if numi1Indices == 6
        codebookType = 'Type1MultiPanel';
    else
        codebookType = 'Type1SinglePanel';
    end
    
    % Extract wideband PMI indices (i1 values) for slots where CSI-RS is
    % present
    i1PerfectValsActiveSlots = reshape([pmiPerfectPerSlot(activeSlotNum).i1],numi1Indices,[])';
    i1PracticalValsActiveSlots = reshape([pmiPracticalPerSlot(activeSlotNum).i1],numi1Indices,[])';
    
    % Extract RI values for slots where CSI-RS is present
    RIPerfectValsActiveSlots = riPerfectPerSlot(activeSlotNum)';
    RIPracticalValsActiveSlots = riPracticalPerSlot(activeSlotNum)';
    
    if isempty(i1PerfectValsActiveSlots)
        disp('No PMI and RI data to plot, because all PMI and RI values are NaNs.');
        return;
    end
    
    figure;
    subplot(4,1,1)
    plot(RIPerfectValsActiveSlots,'r-o');
    hold on;
    plot(RIPracticalValsActiveSlots,'b-*');
    xlabel('Slots')
    ylabel('RI Values');
    xticks(1:numActiveSlots);
    xTickLables = num2cell(activeSlots(:));
    xticklabels(xTickLables);
    [~,upperBound] = bounds([RIPerfectValsActiveSlots; RIPracticalValsActiveSlots]);
    xlim([0 numActiveSlots+8]);
    ylim([0 upperBound+2]);
    yticks(0:upperBound+2);
    title('RI Values')
    legend({'Perfect channel est.','Practical channel est.'});
    
    % Get i11 indices
    i11PerfectVals = i1PerfectValsActiveSlots(:,1);
    i11PracticalVals = i1PracticalValsActiveSlots(:,1);
    
    % Plot i11 values
    subplot(4,1,2)
    plot(i11PerfectVals,'r-o');
    hold on;
    plot(i11PracticalVals,'b-*');
    xlabel('Slots')
    ylabel('i11 Indices');
    xticks(1:numActiveSlots);
    xticklabels(xTickLables);
    [lowerBound,upperBound] = bounds([i11PerfectVals; i11PracticalVals]);
    xlim([0 numActiveSlots+8]);
    ylim([lowerBound-2 upperBound+2]);
    title('PMI: i11 Indices')
    legend({'Perfect channel est.','Practical channel est.'});
    
    % Get i12 indices
    i12PerfectVals = i1PerfectValsActiveSlots(:,2);
    i12PracticalVals = i1PracticalValsActiveSlots(:,2);
    
    % Plot i12 values
    subplot(4,1,3)
    plot(i12PerfectVals,'r-o');
    hold on;
    plot(i12PracticalVals,'b-*');
    xlabel('Slots')
    ylabel('i12 Indices');
    xticks(1:numActiveSlots);
    xticklabels(xTickLables);
    [lowerBound,upperBound] = bounds([i12PerfectVals; i12PracticalVals]);
    xlim([0 numActiveSlots+8]);
    ylim([lowerBound-2 upperBound+2]);
    title('PMI: i12 Indices')
    legend({'Perfect channel est.','Practical channel est.'});
    
    % Get i13 indices
    i13PerfectVals = i1PerfectValsActiveSlots(:,3);
    i13PracticalVals = i1PracticalValsActiveSlots(:,3);
    
    % Plot i13 values
    subplot(4,1,4)
    plot(i13PerfectVals,'r-o');
    hold on;
    plot(i13PracticalVals,'b-*');
    xlabel('Slots')
    ylabel('i13 Indices');
    xticks(1:numActiveSlots);
    xticklabels(xTickLables);
    [lowerBound,upperBound] = bounds([i13PerfectVals; i13PracticalVals]);
    xlim([0 numActiveSlots+8]);
    ylim([lowerBound-2 upperBound+2]);
    title('PMI: i13 Indices')
    legend({'Perfect channel est.','Practical channel est.'});
    
    % Plot the i141, i142 and i143 indices in type 1 multipanel case
    if strcmpi(codebookType,'Type1MultiPanel')
        figure()
        % Get i141 indices
        i141PerfectVals = i1PerfectValsActiveSlots(:,4);
        i141PracticalVals = i1PracticalValsActiveSlots(:,4);
    
        % Plot i141 values
        subplot(3,1,1)
        plot(i141PerfectVals,'r-o');
        hold on;
        plot(i141PracticalVals,'b-*');
        xlabel('Slots')
        ylabel('i141 Indices');
        xticks(1:numActiveSlots);
        xticklabels(xTickLables);
        [lowerBound,upperBound] = bounds([i141PerfectVals; i141PracticalVals]);
        xlim([0 numActiveSlots+8]);
        ylim([lowerBound-2 upperBound+2]);
        title('PMI: i141 Indices')
        legend({'Perfect channel est.','Practical channel est.'});
    
        % Get i142 indices
        i142PerfectVals = i1PerfectValsActiveSlots(:,5);
        i142PracticalVals = i1PracticalValsActiveSlots(:,5);
    
        % Plot i142 values
        subplot(3,1,2)
        plot(i142PerfectVals,'r-o');
        hold on;
        plot(i142PracticalVals,'b-*');
        xlabel('Slots')
        ylabel('i142 Indices');
        xticks(1:numActiveSlots);
        xticklabels(xTickLables);
        [lowerBound,upperBound] = bounds([i142PerfectVals; i142PracticalVals]);
        xlim([0 numActiveSlots+8]);
        ylim([lowerBound-2 upperBound+2]);
        title('PMI: i142 Indices')
        legend({'Perfect channel est.','Practical channel est.'});
    
        % Get i143 indices
        i143PerfectVals = i1PerfectValsActiveSlots(:,6);
        i143PracticalVals = i1PracticalValsActiveSlots(:,6);
    
        % Plot i143 values
        subplot(3,1,3)
        plot(i143PerfectVals,'r-o');
        hold on;
        plot(i143PracticalVals,'b-*');
        xlabel('Slots')
        ylabel('i143 Indices');
        xticks(1:numActiveSlots);
        xticklabels(xTickLables);
        [lowerBound,upperBound] = bounds([i143PerfectVals; i143PracticalVals]);
        xlim([0 numActiveSlots+8]);
        ylim([lowerBound-2 upperBound+2]);
        title('PMI: i143 Indices')
        legend({'Perfect channel est.','Practical channel est.'});
    end
    
    numSubbands = size(pmiPracticalPerSlot(activeSlotNum(1)).i2,2);
    
    % Get the number of i2 indices according to codebook type
    numi2Indices = 1;
    if strcmpi(codebookType,'Type1MultiPanel')
        numi2Indices = 3;
    end
    
    % Extract i2 values
    i2PerfectVals = reshape([pmiPerfectPerSlot(activeSlotNum).i2],[numSubbands,numi2Indices,numActiveSlots]);     % Of size numActiveSlots-by-numi2Indices-numSubbands
    i2PracticalVals = reshape([pmiPracticalPerSlot(activeSlotNum).i2],[numSubbands,numi2Indices,numActiveSlots]); % Of size numActiveSlots-by-numi2Indices-numSubbands
    
    
    if numSubbands > 1
        mode = 'Subband';
    else
        mode = 'Wideband';
    end
    
    % Plot i2 values
    if strcmpi(mode,'Wideband')
        figure;
    
        % In type 1 single panel case, there is only one i2 index. The
        % first column of i2PerfectVals and i2PracticalVals corresponds to
        % i2 index. In type 1 multipanel case, there i2 values are a set
        % of three indices i20, i21 and i22. The first column of in this
        % case i2PerfectVals and i2PracticalVals corresponds to i20 index.
        % Extract the respective index values
        i20PerfectVals = reshape(i2PerfectVals(:,1,:),[],numActiveSlots).';
        i20PracticalVals = reshape(i2PracticalVals(:,1,:),[],numActiveSlots).';
    
        % Plot i2 values or i20 values respectively
        subplot(numi2Indices,1,1)
        plot(i20PerfectVals,'r-o');
        hold on;
        plot(i20PracticalVals,'b-*');
        if strcmpi(codebookType,'Type1SinglePanel')
            title('PMI: i2 Indices')
            xlabel('Slots')
            ylabel('i2 Indices');
        else
            title('PMI: i20 Indices')
            xlabel('Slots')
            ylabel('i20 Indices');
        end
        xticks(1:numActiveSlots);
        xticklabels(xTickLables);
        [lowerBound,upperBound] = bounds([i20PerfectVals; i20PracticalVals]);
        xlim([0 numActiveSlots+8]);
        ylim([lowerBound-1 upperBound+1]);
        legend({'Perfect channel est.','Practical channel est.'});
    
        if strcmpi(codebookType,'Type1MultiPanel')
            % Get i21 values in each slot
            i21PerfectVals = reshape(i2PerfectVals(:,2,:),[],numActiveSlots).';
            i21PracticalVals = reshape(i2PracticalVals(:,2,:),[],numActiveSlots).';
    
            % Plot i21 values
            subplot(numi2Indices,1,2)
            plot(i21PerfectVals,'r-o');
            hold on;
            plot(i21PracticalVals,'b-*');
            title('PMI: i21 Indices')
            xlabel('Slots')
            ylabel('i21 Indices');
            xticks(1:numActiveSlots);
            xticklabels(xTickLables);
            [lowerBound,upperBound] = bounds([i21PerfectVals; i21PracticalVals]);
            xlim([0 numActiveSlots+8]);
            ylim([lowerBound-1 upperBound+1]);
            legend({'Perfect channel est.','Practical channel est.'});
    
            % Get i22 values in each slot
            i22PerfectVals = reshape(i2PerfectVals(:,3,:),[],numActiveSlots).';
            i22PracticalVals = reshape(i2PracticalVals(:,3,:),[],numActiveSlots).';
    
            % Plot i22 values
            subplot(numi2Indices,1,3)
            plot(i22PerfectVals,'r-o');
            hold on;
            plot(i22PracticalVals,'b-*');
            title('PMI: i22 Indices')
            xlabel('Slots')
            ylabel('i22 Indices');
            xticks(1:numActiveSlots);
            xticklabels(xTickLables);
            [lowerBound,upperBound] = bounds([i22PerfectVals; i22PracticalVals]);
            xlim([0 numActiveSlots+8]);
            ylim([lowerBound-1 upperBound+1]);
            legend({'Perfect channel est.','Practical channel est.'});
        end
    else
        if any(nslot == activeSlots)
    
            % In subband mode, plot the PMI i2 indices corresponding to the
            % specified slot number
            figure;
            % Plot i20 values
            subplot(numi2Indices,1,1)
            plot(pmiPerfectPerSlot(nslot+1).i2(1,:),'r-o');
            hold on;
            plot(pmiPracticalPerSlot(nslot+1).i2(1,:),'b-*');
            if strcmpi(codebookType,'Type1SinglePanel')
                title(['PMI: i2 Indices for All Subbands in Slot ' num2str(nslot)]);
                xlabel('Subbands')
                ylabel('i2 Indices');
            else
                title(['PMI: i20 Indices for All Subbands in Slot ' num2str(nslot)]);
                xlabel('Subbands')
                ylabel('i20 Indices');
            end
            xticks(1:numSubbands);
            xticklabels(num2cell(1:numSubbands));
            [lowerBound,upperBound] = bounds([pmiPerfectPerSlot(nslot+1).i2(1,:) pmiPracticalPerSlot(nslot+1).i2(1,:)]);
            yticks(lowerBound:upperBound);
            yticklabels(num2cell(lowerBound:upperBound));
            xlim([0 numSubbands+1])
            ylim([lowerBound-1 upperBound+1]);
            legend({'Perfect channel est.','Practical channel est.'});
    
            if strcmpi(codebookType,'Type1MultiPanel')
                % Plot i21 values
                subplot(numi2Indices,1,2)
                plot(pmiPerfectPerSlot(nslot+1).i2(2,:),'r-o');
                hold on;
                plot(pmiPracticalPerSlot(nslot+1).i2(2,:),'b-*');
                title(['PMI: i21 Indices for All Subbands in Slot ' num2str(nslot)]);
                xlabel('Subbands')
                ylabel('i21 Indices');
                xticks(1:numSubbands);
                xticklabels(num2cell(1:numSubbands));
                [lowerBound,upperBound] = bounds([pmiPerfectPerSlot(nslot+1).i2(2,:) pmiPracticalPerSlot(nslot+1).i2(2,:)]);
                yticks(lowerBound:upperBound);
                yticklabels(num2cell(lowerBound:upperBound));
                xlim([0 numSubbands+1])
                ylim([lowerBound-1 upperBound+1]);
                legend({'Perfect channel est.','Practical channel est.'});
    
                % Plot i22 values
                subplot(numi2Indices,1,3)
                plot(pmiPerfectPerSlot(nslot+1).i2(3,:),'r-o');
                hold on;
                plot(pmiPracticalPerSlot(nslot+1).i2(3,:),'b-*');
                title(['PMI: i22 Indices for All Subbands in Slot ' num2str(nslot)]);
                xlabel('Subbands')
                ylabel('i22 Indices');
                xticks(1:numSubbands);
                xticklabels(num2cell(1:numSubbands));
                [lowerBound,upperBound] = bounds([pmiPerfectPerSlot(nslot+1).i2(3,:) pmiPracticalPerSlot(nslot+1).i2(3,:)]);
                yticks(lowerBound:upperBound);
                yticklabels(num2cell(lowerBound:upperBound));
                xlim([0 numSubbands+1])
                ylim([lowerBound-1 upperBound+1]);
                legend({'Perfect channel est.','Practical channel est.'});
            end
        else
            disp(['For the specified slot (' num2str(nslot) '), PMI i2 indices are not reported. Please choose another slot number.'])
        end
    end
end