% Configure three carriers for the aggregation. You can select a different
% number of carriers by modifying the number of elements in the
% |channelBandwidths|, |SCSs|, |modulations|, and |nCellIDs| vectors.
frequencyRange = 'FR1'; % (FR1 or FR2)
channelRaster = 100; % Channel raster in kHz (15, 60 or 100)
channelBandwidths = [60 40 40]; % Channel bandwidths in MHz 
                                % (5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 
                                % 90 or 100 for FR1)
                                % (50, 100, 200 or 400 for FR2)
SCSs = [30 30 30]; % Subcarrier spacings in kHz 
                   % (15, 30 or 60 for FR1) 
                   % (60 or 120 for FR2)
modulations = {'QPSK' '64QAM' '256QAM'}; % ('QPSK', '64QAM' or '256QAM' for FR1)
                                         % ('QPSK', '16QAM' or '64QAM' for FR2)
nCellIDs = [1 2 3]; % Cell IDs          
duplexMode = 'TDD'; % Duplex mode ('TDD' or 'FDD')
sv = '15.2.0'; % Standard version ('15.1.0', '15.2.0' or '15.7.0')
numSubframes = 10; % Number of subframes to generate per carrier

% Choose the spacings between consecutive carriers, |ccSpacings|, or leave
% the vector empty for calculating the spacings, as described in TS 38.104
% Section 5.3A
ccSpacings = []; % MHz

% Select the CC to demodulate
CCofInterest = 2;

% Verify the number of CCs, channel bandwidths, subcarrier spacings,
% modulations, and spacings. Additionally, check that the spacings are
% greater than 0.
hNRVerifyCarrierParameters(length(channelBandwidths),length(SCSs),...
    length(modulations),length(nCellIDs),ccSpacings);

% Establish the number of component carriers 
numCC = length(channelBandwidths);

% CC configuration
referenceChannel = cell(1,numCC);
wavegen = cell(1,numCC);
for i = 1:numCC
    % Select a reference channel, FRC, based on the chosen modulation and
    % frequency range
    referenceChannel{i} = strcat('DL-FRC-',frequencyRange,'-',modulations{i});
    
    % Create a generator object for the above PDSCH FRC reference channel
    wavegen{i} = hNRReferenceWaveformGenerator(referenceChannel{i},...
        channelBandwidths(i),SCSs(i),duplexMode,nCellIDs(i),sv);
end

% Get the largest SCS configuration, |SCSConfig|, among the SCS
% configurations supported for each two consecutive channel bandwidths to
% obtain the minimum guardband, as described in TS 38.104 Section 5.4.1.2.
table = hGetGBTable(frequencyRange); % Minimum guardband table
if isequal(frequencyRange,'FR1')
    if any(channelBandwidths(:) == 5)
        SCSConfig = 1; 
    else
        SCSConfig = 2; 
    end
else
    SCSConfig = 3;
end

% Calculate the minimum guardband, as described in TS 38.104 Section 5.3.3
% Table 5.3.3-1 for FR1 and Table 5.3.3-1 for FR2, to obtain the minimum CC
% spacings.
minimumGuardBand = zeros(1,numCC);
NDLRB = zeros(1,numCC);
scs = strcat(num2str(2^SCSConfig*15),'kHz'); % Common SCS in kHz to extract 
                                             % minimum guardband 
for i = 1:numCC 
    NDLRB(i) = wavegen{i}.Config.SCSCarriers{1}.NSizeGrid;
    minimumGuardBand(i) = table{{scs},{strcat(num2str(channelBandwidths(i)),...
        'MHz');}}; % kHz
    minimumGuardBand(i) = minimumGuardBand(i)*1e-3; % MHz
end

% Compute the minimum CC spacings, as defined in TS 38.104 Section
% 5.4.1.2. Use these spacings to calculate the center frequencies and for
% filtering each CC.
minCCSpacings = zeros(1,numCC-1); % Minimum CC spacing
for k = 2:numCC
    minCCSpacings(k-1) = hNRCarrierAggregationChannelSpacing( ...
        channelBandwidths(k-1), channelBandwidths(k), minimumGuardBand(k-1), ...
        minimumGuardBand(k),channelRaster,SCSConfig); % MHz
end

% Determine the center frequency for each CC with respect to 0 Hz.
% Initially, the lower carrier frequency is at baseband (Fc_offset(1) = 0 Hz).
Fc_offset = zeros(1,numCC);
if isempty(ccSpacings)
    ccSpacings = minCCSpacings;
end
for k = 2:numCC
    Fc_offset(k) = Fc_offset(k-1) + ccSpacings(k-1); % MHz
end

F_offset_low = (NDLRB(1)*12+1)*(SCSs(1)*1e-3)/2 + minimumGuardBand(1); % MHz
F_offset_high = (NDLRB(end)*12-1)*(SCSs(end)*1e-3)/2 + minimumGuardBand(end); %MHz

F_edge_low = Fc_offset(1) - F_offset_low; % MHz
F_edge_high = Fc_offset(end) + F_offset_high; % MHz

BW_channel_CA = F_edge_high - F_edge_low; % MHz
fprintf('BW_channel_CA: %0.4f MHz\n',BW_channel_CA);

shiftToCenter = -1*(BW_channel_CA/2 + F_edge_low);

% Center aggregated bandwidth at baseband
Fc_offset = Fc_offset + shiftToCenter;
F_edge_low = Fc_offset(1) - F_offset_low;
F_edge_high = Fc_offset(end) + F_offset_high;

% Display frequency band edges
fprintf('F_edge_low:  %0.4f MHz\n',F_edge_low);

fprintf('F_edge_high: %0.4f MHz\n',F_edge_high);

fprintf('F_offset_low:  %0.4f MHz\n',F_offset_low);

fprintf('F_offset_high: %0.4f MHz\n',F_offset_high);

% Display carrier frequencies
fprintf('\n');
for i = 1:numCC
    fprintf('Component Carrier %d:\n',i);
    fprintf('   Fc: %0.4f MHz\n', Fc_offset(i));
end

% Obtain sample rates of the component carriers
CCSR = zeros(1,numCC);
carriers = cell(1,numCC);
for i = 1:numCC
    carriers{i} = nrCarrierConfig;
    carriers{i}.NCellID = nCellIDs(i);
    carriers{i}.NSizeGrid = NDLRB(i);
    carriers{i}.SubcarrierSpacing = SCSs(i);
    carriers{i}.CyclicPrefix = wavegen{i}.Config.BandwidthParts{1}.CyclicPrefix;
    info = nrOFDMInfo(carriers{i});
    CCSR(i) = info.SampleRate; % Hz
end

% Calculate the oversampling ratio for the largest BW CC to ensure the
% waveform occupies a maximum of 85% of the total bandwidth, |bwfraction|
bwfraction = 0.85; % Bandwidth utilization of 85%
OSR = (BW_channel_CA/bwfraction)/(max(CCSR)/1e6);

% To simplify the resampling operation, choose an oversampling ratio which
% is a power of 2: calculate the next power of two above OSR
OSR = 2^ceil(log2(OSR));

% Calculate the overall sample rate for the aggregated waveform
SR = OSR*max(CCSR); % Hz
fprintf('\nOutput sample rate: %0.4f Ms/s\n\n',SR/1e6);

% Calculate the individual oversampling factors for the component carriers
OSRs = SR./CCSR;

% Generate and aggregate the component carriers
tmwaveinfo = cell(1,numCC);
waveInfo = cell(1,numCC);
resourcesInfo = cell(1,numCC);
clear waveform
for i = numCC:-1:1
    % Generate each CC
    [wf,waveInfo{i},resourcesInfo{i}] = generateWaveform(wavegen{i},...
        numSubframes);
    
    % Resample the CCs so that they have the same sample rate  
    waveform(:,i) = resample(wf,OSRs(i),1)/OSRs(i);    
end

% Aggregate all CC. comm.MultibandCombiner applies the frequency offsets
% and combines the resulting signals together.
carrierAggregator = comm.MultibandCombiner(InputSampleRate = SR,...
    FrequencyOffsets = Fc_offset*1e6,...
    OutputSampleRateSource = 'Property',...
    OutputSampleRate = SR);

waveform = carrierAggregator(waveform);

specPlot = hNRCarrierAggregationPlotSpectrum(waveform,SR,...
    'Power Spectrum of Carrier Aggregation',{'Signal spectrum'});

% Verify carrier of interest ID
if CCofInterest > numCC || CCofInterest <= 0 || mod(CCofInterest,1) ~= 0
    error('nr5g:NRDownlinkCarrierAggregationExample:CCOutOfRange',...
        'Cannot demodulate CC number %d, choose an integer number that falls between 1 and %d\n',...
        CCofInterest, numCC) ;
end
fprintf(1,'Extracting CC number %d: \n', CCofInterest);

% Define downsampling filter order
filterOrder = 201;

% Precalculate the filter passband and stopband values for all CC
firPassbandVec = (NDLRB*12-1).*(SCSs*1e-3)/2 / (SR/1e6/2);
firStopbandVec = hNRCarrierAggregationStopband(minCCSpacings,NDLRB,SR,SCSs);

% Choose the passband and stopband values for the carrier of interest
firPassband = firPassbandVec(CCofInterest);
firStopband = firStopbandVec(CCofInterest);

% Pad signal with zeros to consider filter transient response length
waveform = [waveform; zeros(filterOrder*2,size(waveform,2))];

% Center the carrier of interest at 0 Hz
frequencyShifter = comm.PhaseFrequencyOffset(SampleRate = SR,...
    FrequencyOffset = -Fc_offset(CCofInterest)*1e6);
demodulatedCC = frequencyShifter(waveform);

% To ease the filter design requirements, apply the downsampling in two
% stages if necessary. Use a downsampling factor of 4 in the initial stage.
% If the quality of the resulting signal is not as required, consider a
% different filter design in this initial stage.
if (firStopband < 0.1)
    % Downsample by 4 in the initial stage
    SRC = 4;
    demodulatedCC = resample(demodulatedCC,1,SRC);
    % Update passband and stopband values
    firPassband = firPassband * SRC;
    firStopband = firStopband * SRC;
else
    % Do not apply an extra downsampling 
    SRC = 1;
end

% Design the lowpass filter to filter out the CC of your choice
frEdges = [0 firPassband firStopband 1];
firFilter = dsp.FIRFilter;
firFilter.Numerator = firpm(filterOrder,frEdges,[1 1 0 0]);

% Display the response of the designed filter
fvtool(firFilter,'Analysis','freq');

% Filter the signal to extract the component of interest
rxWaveform  = firFilter(demodulatedCC);

% Plot the demodulated and filtered spectra
filteredSpecPlot = ...
    hNRCarrierAggregationPlotSpectrum([demodulatedCC, rxWaveform],SR,...
            'Demodulated and Filtered Waveform Power Spectrum',...
            {'Carrier aggregated signal' 'Filtered signal'});

% Downsample the filtered carrier to its baseband rate
rxWaveform = downsample(rxWaveform,OSRs(CCofInterest)/SRC);

% Parameterize the channel estimator configuration using the structure |cfg|
cfg = struct();
cfg.Evm3GPP = true; % To measure EVM as defined in TS 38.104, Annex B(FR1) 
                    % / Annex C(FR2) set |Evm3GPP| to |true|. 
cfg.TargetRNTIs = [];
cfg.PlotEVM = true;
cfg.DisplayEVM = true;
cfg.Label = wavegen{CCofInterest}.ConfiguredModel{1};

% Perform EVM measurements and plot results
[evmInfo,eqSym,refSym] = hNRPDSCHEVM(wavegen{CCofInterest}.Config,...
    rxWaveform,cfg);

% Perform time synchronization on the input waveform
offset = nrTimingEstimate(carriers{CCofInterest},rxWaveform,...
    waveInfo{CCofInterest}.ResourceGridBWP);
rxWaveform = rxWaveform(1+offset:end,:);

% Perform OFDM demodulation
rxGrid = nrOFDMDemodulate(carriers{CCofInterest},rxWaveform);

% Get the allocated slots and OFDM symbols per slot
allocatedSlots = zeros(1,size(resourcesInfo{CCofInterest}.WaveformResources.PDSCH.Resources,2));
for i=1:length(allocatedSlots)
    allocatedSlots(i) = resourcesInfo{CCofInterest}.WaveformResources.PDSCH.Resources(i).NSlot;
end
L = carriers{CCofInterest}.SymbolsPerSlot; % OFDM symbols per slot

% Create a DLSCH decoder System object
decodeDLSCH = nrDLSCHDecoder;
decodeDLSCH.MultipleHARQProcesses = false;
decodeDLSCH.TargetCodeRate = wavegen{CCofInterest}.Config.PDSCH{1}.TargetCodeRate;
decodeDLSCH.LDPCDecodingAlgorithm = 'Normalized min-sum';
decodeDLSCH.MaximumLDPCIterationCount = 6;

for NSlot=1:length(allocatedSlots)
    % Extract slot
    SlotID = allocatedSlots(NSlot);
    rxSlot = rxGrid(:,(1:L)+(SlotID*L),:);
    refSlot = waveInfo{CCofInterest}.ResourceGridBWP(:,(1:L)+(SlotID*L),:);
    
    % Perform channel estimation
    [estChannelGrid,noiseEst] = nrChannelEstimate(carriers{CCofInterest},...
        rxSlot,refSlot);
            
    % Get PDSCH resource elements from the received grid and channel estimate
    pdschIndices = resourcesInfo{CCofInterest}.WaveformResources.PDSCH.Resources(NSlot).ChannelIndices;
    [pdschRx,pdschHest] = nrExtractResources(pdschIndices,rxSlot,estChannelGrid);
    
    % Perform equalization
    [pdschEq,csi] = nrEqualizeMMSE(pdschRx,pdschHest,noiseEst);
    
    % Perform layer demapping, symbol demodulation, and descrambling
    modulation = wavegen{CCofInterest}.Config.PDSCH{1}.Modulation;
    RNTI = wavegen{CCofInterest}.Config.PDSCH{1}.RNTI;
    [dlschLLRs,rxSymbols] = nrPDSCHDecode(pdschEq,modulation,...
        carriers{CCofInterest}.NCellID,RNTI,noiseEst);
    
    % Scale LLRs by CSI
    csi = nrLayerDemap(csi); % CSI layer demapping
    NumCW = size(resourcesInfo{CCofInterest}.WaveformResources.PDSCH.Resources(NSlot).Codeword,2);
    for cwIdx = 1:NumCW
        Qm = length(dlschLLRs{cwIdx})/length(rxSymbols{cwIdx}); % Bits per symbol
        csi{cwIdx} = repmat(csi{cwIdx}.',Qm,1);   % Expand by each bit per symbol
        dlschLLRs{cwIdx} = dlschLLRs{cwIdx} .* csi{cwIdx}(:);   % Scaled symbols
    end  
    
    % Calculate the transport block sizes for the codewords in the slot
    trBlkSize = resourcesInfo{CCofInterest}.WaveformResources.PDSCH.Resources(NSlot).TransportBlockSize;
    
    % Decode the DL-SCH transport channel
    decodeDLSCH.TransportBlockLength = trBlkSize;
    NLayers = wavegen{CCofInterest}.Config.PDSCH{1}.NumLayers;
    RVSequence = wavegen{CCofInterest}.Config.PDSCH{1}.RVSequence;
    [decbits,crc] = decodeDLSCH(dlschLLRs,modulation,NLayers,RVSequence);
     
    % Display the CRC status
    if crc
        disp(['Slot ' num2str(SlotID) ': CRC failed']);
    else
        disp(['Slot ' num2str(SlotID) ': CRC passed']);
    end
end