% Configure the set of NULRB values to describe the carriers to be
% aggregated
NULRB = [50 75 100];

% Establish the number of component carriers
numCC = numel(NULRB);

CCBWs = zeros(1,numCC);

if numCC<2
    error('Please specify more than one CC bandwidth value.');
end

% Specify the number of subframes to generate
numSubframes = 10;

% Configure CCs
frc = cell(1,numCC);
for i = 1:numCC
    switch NULRB(i)
        case 6
            frc{i} = lteRMCUL('A3-1');
            frc{i}.PUSCH.PRBSet = (0:2).';
        case 15
            frc{i} = lteRMCUL('A1-2');
            frc{i}.PUSCH.PRBSet = (0:8).';
        case 25
            frc{i} = lteRMCUL('A1-3');
            frc{i}.PUSCH.PRBSet = (0:15).';
        case 50
            frc{i} = lteRMCUL('A3-5');
            frc{i}.PUSCH.PRBSet = (0:39).';
        case 75
            frc{i} = lteRMCUL('A3-6');
            frc{i}.PUSCH.PRBSet = (0:59).';
        case 100
            frc{i} = lteRMCUL('A3-7');
            frc{i}.PUSCH.PRBSet = (0:80).';
        otherwise
            fprintf('Not a valid number of resource blocks: %d\n',...
                NULRB(i));
            return
    end
    CCBWs(i) = hNRBToBandwidth(NULRB(i));
    frc{i}.Bandwidth = CCBWs(i);
    frc{i}.TotSubframes = numSubframes;
    frc{i}.PUSCH.RVSeq = 0;
end

cec = struct;                        % Channel estimation config structure
cec.PilotAverage = 'UserDefined';    % Type of pilot symbol averaging
cec.FreqWindow = 15;                 % Frequency window size
cec.TimeWindow = 15;                 % Time window size
cec.InterpType = 'Cubic';            % 2D interpolation type
cec.InterpWindow = 'Centered';       % Interpolation window type
cec.InterpWinSize = 1;               % Interpolation window size
cec.Reference = 'Antennas';          % Specifies point of reference for
                                     % channel estimation.

% Define Delta_f1 parameter in MHz for nominal guard band and frequency
% offset calculations. In the uplink Delta_f1 is zero
% (TS 36.101, Table 5.6A-1)
Delta_f1 = 0; % in MHz
maxBW = max(CCBWs);

% Initialize center frequency and CC spacing vectors
F_c = zeros(1,numCC); % Center frequencies
ccSpacing = zeros(1,numCC-1); % CC spacing

% Calculate nominal guard band, TS 36.101 5.6A-1
nominalGuardBand = 0.05*maxBW-0.5*Delta_f1;

% Initially assume lower carrier frequency is at baseband
F_c(1) = 0;

% Calculate CC spacing and center frequencies
for k = 2:numCC
    ccSpacing(k-1) = hCarrierAggregationChannelSpacing( ...
        frc{k-1}.Bandwidth, frc{k}.Bandwidth);
    F_c(k) = F_c(k-1) + ccSpacing(k-1);
end

% Calculate lower and higher frequency offsets, TS 36.101 5.6A
F_offset_low = (0.18*NULRB(1)+Delta_f1)/2 + nominalGuardBand;
F_offset_high = (0.18*NULRB(end)+Delta_f1)/2 + nominalGuardBand;

% Calculate lower and higher frequency edges, TS 36.101 5.6A
F_edge_low = F_c(1) - F_offset_low;
F_edge_high = F_c(end) + F_offset_high;

% Calculate aggregated channel bandwidth, TS 36.101 5.6A
BW_channel_CA = F_edge_high - F_edge_low;
fprintf('BW_channel_CA: %0.4f MHz\n',BW_channel_CA);

% Calculate shift to center baseband
shiftToCenter = -1*(BW_channel_CA/2 + F_edge_low);

% Shift the center frequencies so that the aggregated bandwidth is centered
% at baseband
F_c = F_c + shiftToCenter;
F_edge_low = F_c(1) - F_offset_low;
F_edge_high = F_c(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', F_c(i));
end

% Considering the bandwidth utilization as 85%
bwfraction = 0.85;

% Calculate sampling rates of the component carriers
CCSR = zeros(1,numCC);
for i = 1:numCC
    info = lteSCFDMAInfo(frc{i});
    CCSR(i) = info.SamplingRate;
end

% Calculate overall sampling rate for the aggregated signal

% Calculate the oversampling ratio (for the largest BW CC) to make sure the
% signal occupies 85% (bwfraction) of the aggregated channel bandwidth
OSR = (BW_channel_CA/bwfraction)/(max(CCSR)/1e6);
% To simplify the resampling operation, choose an oversampling ratio as the
% smallest power of 2 greater than or equal to OSR
OSR = 2^ceil(log2(OSR));
SR = OSR*max(CCSR);
fprintf('\nOutput sample rate: %0.4f Ms/s\n\n',SR/1e6);

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

% Generate component carriers
tx = cell(1,numCC);
for i = 1:numCC
    tx{i} = lteRMCULTool(frc{i},randi([0 1],1000,1));
    tx{i} = resample(tx{i},OSRs(i),1)/OSRs(i);
end

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

waveform = sigAggregator([tx{:}]);

specPlot = hPlotSpectrum(waveform,SR,...
    'Power Spectrum of Carrier Aggregation Waveform',...
    {'Carrier aggregated signal'});

% Select CC to demodulate
CCofInterest = 1;
if CCofInterest > numCC || CCofInterest <= 0
    error('Cannot demodulate CC number %d, there are %d CCs\n',...
        CCofInterest, numCC) ;
end

% Define downsampling filter order
filterOrder = 301;

% Precalculate passband and stopband values for all CC
firPassbandVec = (0.18*NULRB+Delta_f1)/2 / (SR/1e6/2);
firStopbandVec = hCarrierAggregationStopband(ccSpacing,NULRB,SR);

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

% Extract and decode CC of interest
fprintf(1,'Extracting CC number %d: \n', CCofInterest);

% Pad signal with zeros to take into account filter transient response
% length
waveform = [waveform; zeros(filterOrder+1,size(waveform,2))];

% Demodulate carrier of interest
demodulatedCC = ...
    hCarrierAggregationModulate(waveform,SR,-F_c(CCofInterest)*1e6);

% Downsampling is done in two stages if the filter is too narrow. This
% eases the filter design requirements. If this is the case an initial
% downsampling factor of 4 is applied. You may want to consider a different
% filter design in this initial stage if the quality of the resulting
% signal is not acceptable.
if (firStopband < 0.1)
    % Down-sample by 4
    SRC = 4;
    demodulatedCC = resample(demodulatedCC,1,SRC);
    % Update pass and stopband to take first downsampling into account
    firPassband = firPassband * SRC;
    firStopband = firStopband * SRC;
else
    % No pre-filter
    SRC = 1;
end

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

fvtool(firFilter,'Analysis','freq');

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

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

rxWaveform = downsample(rxWaveform,OSRs(CCofInterest)/SRC);

% Get the parameters for CC of interest
FRC = frc{CCofInterest};
PUSCH = FRC.PUSCH;

% Synchronize received waveform
offset = lteULFrameOffset(FRC,PUSCH,rxWaveform);
rxWaveform = rxWaveform(1+offset:end, :);

% Configure channel estimation structure for EVM measurement
cecEVM = cec;
cecEVM.PilotAverage = 'TestEVM';
[evmpusch, evmdrs, emissions, plots] = hPUSCHEVM(frc{CCofInterest},cecEVM,rxWaveform);

% Plot the absolute in-band emissions
if (~isempty(emissions.DeltaRB))
    hPUSCHEVMEmissionsPlot(emissions);
end

% Perform SC-FDMA demodulation
rxGrid = lteSCFDMADemodulate(FRC,rxWaveform);

% Get the number of received subframes and OFDM symbols per subframe
dims = lteSCFDMAInfo(FRC);
samplesPerSubframe = dims.SamplingRate/1000;
nRxSubframes = floor(size(rxWaveform, 1)/samplesPerSubframe);
FRC.TotSubframes = 1;
resGridSize = lteResourceGridSize(FRC);
L = resGridSize(2);

disp('Decode transmitted subframes and check CRC.');

for n=0:nRxSubframes-1

    % Extract subframe
    rxSubframe = rxGrid(:,(1:L)+(n*L),:);

    % Get the transport block size for current subframe
    FRC.NSubframe = n;
    trBlkSize = PUSCH.TrBlkSizes(n+1);

    % Perform channel estimation
    [estChannelGrid, noiseEst] = lteULChannelEstimate( ...
        FRC, PUSCH, cec, rxSubframe);

    % Generate PUSCH indices and extract symbols from received and
    % channel estimate grids in preparation for PUSCH decoding
    ind = ltePUSCHIndices(FRC, PUSCH);
    [rxSym,hestSym] = lteExtractResources(ind, rxSubframe, estChannelGrid);

    % Update PUSCH to carry complete information of the UL-SCH coding
    PUSCH = lteULSCHInfo(FRC, PUSCH, trBlkSize, 'chsconcat');

    % Perform equalization, transform deprecoding, layer demapping,
    % demodulation, and descrambling on the received data using the channel
    % estimate
    [rxEncodedBits, puschsymbols] = ltePUSCHDecode(FRC, PUSCH, ...
        rxSym, hestSym, noiseEst);

    % Decode Uplink Shared Channel (UL-SCH)
    decState = [];
    [decbits,crc] = ...
        lteULSCHDecode(FRC, PUSCH, trBlkSize ,rxEncodedBits, decState);
    if crc
        disp(['Subframe ' num2str(n) ': CRC failed']);
    else
        disp(['Subframe ' num2str(n) ': CRC passed']);
    end
end