% Burst configuration related to the burst structure itself:
txBurst.BlockPattern = 'Case B';
txBurst.SSBPeriodicity = 20;
txBurst.NFrame = 4;
txBurst.SSBTransmitted = [1 1 1 1 1 1 1 1];
txBurst.NCellID = 102;

% Burst configuration related to carrier (10 MHz, see TS 38.104 Table
% 5.3.2-1):
gnb.SubcarrierSpacing = 15;
gnb.NRB = 52;
gnb.CyclicPrefix = 'Normal';
carrierInfo = hOFDMInfo(gnb);
txBurst.SampleRate = carrierInfo.SamplingRate;
K = carrierInfo.NSubcarriers;
txBurst.FrequencyPointA = -K/2 * gnb.SubcarrierSpacing * 1e3;

% Burst configuration related to MIB content:
txBurst.DMRSTypeAPosition = 2;
txBurst.PDCCHConfigSIB1 = 17;
txBurst.CellBarred = 0;
txBurst.IntraFreqReselection = 0;
txBurst.SubcarrierSpacingCommon = carrierInfo.SubcarrierSpacing;
txBurst.DisplayBurst = true;

% Configure number of transmit and receive antennas
ntxants = 8;
nrxants = 2;

% Configure channel
velocity = 30.0;
fc = 4e9;
c = physconst('lightspeed');
fd = (velocity*1000/3600)/c*fc;
channel = nrTDLChannel;
channel.Seed = 24;
channel.DelayProfile = 'TDL-C';
channel.DelaySpread = 300e-9;
channel.MaximumDopplerShift = fd;
channel.MIMOCorrelation = 'Medium';
channel.Polarization = 'Cross-Polar';
channel.NumTransmitAntennas = ntxants;
channel.NumReceiveAntennas = nrxants;
channel.SampleRate = txBurst.SampleRate;

% Configure SNR for AWGN
SNRdB = 10;

% Display burst
disp(txBurst);

% Create and display burst information
txBurstInfo = hSSBurstInfo(txBurst);
disp(txBurstInfo);

% Create burst resource grid
[~,txBurstGrid] = hSSBurst(txBurst);

% Configure beamforming weights
W = fft(eye(ntxants)) / sqrt(ntxants);

% Beamform the OFDM symbols corresponding to each burst
beamformedGrid = zeros([size(txBurstGrid) ntxants]);
blockSubcarriers = txBurstInfo.OccupiedSubcarriers;
for ssb = 1:length(txBurstInfo.SSBIndex)

    blockSymbols = txBurstInfo.OccupiedSymbols(ssb,:);
    block = txBurstGrid(blockSubcarriers,blockSymbols);
    Wssb = W(mod(ssb-1,ntxants)+1,:);
    beamformedBlock = reshape(block(:) * Wssb,[size(block) ntxants]);
    beamformedGrid(blockSubcarriers,blockSymbols,:) = beamformedBlock;

end

% Perform OFDM modulation
beamformedGrid = beamformedGrid(:,1:max(txBurstInfo.OccupiedSymbols(:))+1,:);
txOfdmConfig.SubcarrierSpacing = txBurstInfo.SubcarrierSpacing;
txOfdmConfig.NRB = txBurstInfo.NRB;
txOfdmConfig.CyclicPrefix = txBurstInfo.CyclicPrefix;
txOfdmConfig.Windowing = 0;
[txWaveform,txOfdmInfo] = hOFDMModulate(txOfdmConfig,beamformedGrid);

rxWaveform = channel(txWaveform);

rng('default');
rxWaveform = awgn(rxWaveform,SNRdB,-10*log10(double(txOfdmInfo.Nfft)));

% Configure necessary burst parameters
burst.BlockPattern = txBurst.BlockPattern;
burst.SampleRate = txBurst.SampleRate;
L_max = numel(txBurst.SSBTransmitted);
burst.SSBTransmitted = ones(1,L_max);

% Create burst information structure and OFDM information from configured
% burst parameters
burstInfo = hSSBurstInfo(burst);
ofdmConfig.SubcarrierSpacing = burstInfo.SubcarrierSpacing;
ofdmConfig.NRB = burstInfo.NRB;
ofdmConfig.CyclicPrefix = burstInfo.CyclicPrefix;
ofdmInfo = hOFDMInfo(ofdmConfig);

pssIndices = nrPSSIndices;
pssGrid = zeros([240 4]);
refGrid = zeros([ofdmInfo.NSubcarriers ofdmInfo.SymbolsPerSlot]);
k = burstInfo.OccupiedSubcarriers;

fstep = burstInfo.SubcarrierSpacing * 1e3 / 2; % half subcarrier
fshifts = (-6:6) * fstep;
peak_value = zeros(numel(fshifts),3);
peak_index = zeros(numel(fshifts),3);
t = (0:size(rxWaveform,1)-1).' / burst.SampleRate;
for fIdx = 1:numel(fshifts)

    coarseFrequencyOffset = fshifts(fIdx);
    rxWaveformFreqCorrected = rxWaveform .* exp(-1i*2*pi*coarseFrequencyOffset*t);

    for NID2 = [0 1 2]

        pssRef = nrPSS(NID2);
        pssGrid(pssIndices) = pssRef;
        refGrid(k,2:5) = pssGrid;

        nSlot = 0;
        [~,corr] = nrTimingEstimate(rxWaveformFreqCorrected,ofdmConfig.NRB,ofdmConfig.SubcarrierSpacing,nSlot,refGrid,'CyclicPrefix',ofdmConfig.CyclicPrefix);
        corr = sum(abs(corr),2);
        [peak_value(fIdx,NID2+1),peak_index(fIdx,NID2+1)] = max(corr);
        peak_index(fIdx,NID2+1) = peak_index(fIdx,NID2+1) + ofdmInfo.SymbolLengths(1);

    end

end

% Plot PSS correlations
figure;
hold on;
plot(fshifts/1e3,peak_value);
title('PSS Correlations versus Frequency Offset');
ylabel('Magnitude');
xlabel('Frequency Offset (kHz)');

% Determine NID2 and coarse frequency offset by finding the strongest
% correlation
[fIdx,NID2] = find(peak_value==max(peak_value(:)));
coarseFrequencyOffset = fshifts(fIdx);
NID2 = NID2 - 1;

% Apply coarse frequency correction
rxWaveform = rxWaveform .* exp(-1i*2*pi*coarseFrequencyOffset*t);

% Plot selected NID2
plot(coarseFrequencyOffset/1e3,peak_value(fIdx,NID2+1),'kx','LineWidth',2,'MarkerSize',8);
lgd = legend;
lgd.Interpreter = 'latex';
legends = "$N_{ID}^{(2)}$ = " + num2cell(0:2);
legend([legends "coarse $\Delta_f$ = " + num2str(coarseFrequencyOffset) + ", $N_{ID}^{(2)}$ = " + num2str(NID2)],'Location','East');

% Determine timing offset
offset = peak_index(fIdx,NID2+1) - 1;

% Perform fine frequency offset estimation using CP correlation across
% the 4 OFDM symbols of the SSB
fineFrequencyOffset = frequencyOffsetSSB(rxWaveform(1+offset:end,:),ofdmInfo);

% Apply fine frequency correction
rxWaveform = rxWaveform .* exp(-1i*2*pi*fineFrequencyOffset*t);

% Extract strongest burst
offset = offset - ofdmInfo.SymbolLengths(1);
rxGrid = hOFDMDemodulate(ofdmConfig,rxWaveform(1+offset:end,:));
rxGrid = rxGrid(burstInfo.OccupiedSubcarriers,2:5,:);

% Extract the received SSS symbols from the SS/PBCH block
sssIndices = nrSSSIndices;
sssRx = nrExtractResources(sssIndices,rxGrid);

% Correlate received SSS symbols with each possible SSS sequence
sssEst = zeros(1,336);
for NID1 = 0:335

    ncellid = (3*NID1) + NID2;
    sssRef = nrSSS(ncellid);
    sssEst(NID1+1) = sum(abs(mean(sssRx .* conj(sssRef),1)).^2);

end

% Plot SSS correlations
figure;
stem(0:335,sssEst,'o');
title('SSS Correlations (frequency domain)');
xlabel('$N_{ID}^{(1)}$','Interpreter','latex');
ylabel('Magnitude');
axis([-1 336 0 max(sssEst)*1.1]);

% Determine NID1 by finding the strongest correlation
NID1 = find(sssEst==max(sssEst)) - 1;

% Plot selected NID1
hold on;
plot(NID1,max(sssEst),'kx','LineWidth',2,'MarkerSize',8);
lgd = legend;
lgd.Interpreter = 'latex';
legend(["correlations" "$N_{ID}^{(1)}$ = " + num2str(NID1)]);

% Form overall cell identity from NID1 and NID2
ncellid = (3*NID1) + NID2;

% Calculate PBCH DM-RS indices
dmrsIndices = nrPBCHDMRSIndices(ncellid);

% Perform channel estimation using DM-RS symbols for each possible DM-RS
% sequence and estimate the SNR
dmrsEst = zeros(1,8);
for ibar_SSB = 0:7

    refGrid = zeros([240 4]);
    refGrid(dmrsIndices) = nrPBCHDMRS(ncellid,ibar_SSB);
    [hest,nest] = nrChannelEstimate(rxGrid,refGrid,'CyclicPrefix',ofdmConfig.CyclicPrefix,'AveragingWindow',[0 1]);
    dmrsEst(ibar_SSB+1) = 10*log10(mean(abs(hest(:).^2)) / nest);

end

% Plot PBCH DM-RS SNRs
figure;
stem(0:7,dmrsEst,'o');
title('PBCH DM-RS SNR estimates');
xlabel('$\overline{i}_{SSB}$','Interpreter','latex');
xticks(0:7);
ylabel('Estimated SNR (dB)');
axis([-1 8 min(dmrsEst)-1 max(dmrsEst)+1]);

% Record ibar_SSB for the highest SNR
ibar_SSB = find(dmrsEst==max(dmrsEst)) - 1;

% Plot selected ibar_SSB
hold on;
plot(ibar_SSB,max(dmrsEst),'kx','LineWidth',2,'MarkerSize',8);
lgd = legend;
lgd.Interpreter = 'latex';
legend(["SNRs" "$\overline{i}_{SSB}$ = " + num2str(ibar_SSB)]);

refGrid = zeros([240 4]);
refGrid(dmrsIndices) = nrPBCHDMRS(ncellid,ibar_SSB);
refGrid(sssIndices) = nrSSS(ncellid);
[hest,nest] = nrChannelEstimate(rxGrid,refGrid,'CyclicPrefix',ofdmConfig.CyclicPrefix,'AveragingWindow',[0 1]);

% Extract the received PBCH symbols from the SS/PBCH block
[pbchIndices,pbchIndicesInfo] = nrPBCHIndices(ncellid);
pbchRx = nrExtractResources(pbchIndices,rxGrid);

% Plot received PBCH constellation before equalization
figure;
plot(pbchRx,'o');
title('Received PBCH constellation');
m = max(abs([real(pbchRx(:)); imag(pbchRx(:))])) * 1.1;
axis([-m m -m m]);

% Configure 'v' for PBCH scrambling according to TS 38.211 Section 7.3.3.1
% 'v' is also the 2 LSBs of the SS/PBCH block index for L=4, or the 3 LSBs
% for L=8 or 64
if (burstInfo.L==4)
    v = mod(ibar_SSB,4);
else
    v = ibar_SSB;
end
ssbIndex = v;

% PBCH equalization and CSI calculation
pbchHest = nrExtractResources(pbchIndices,hest);
[pbchEq,csi] = nrEqualizeMMSE(pbchRx,pbchHest,nest);
Qm = pbchIndicesInfo.G / pbchIndicesInfo.Gd;
csi = repmat(csi.',Qm,1);
csi = reshape(csi,[],1);

% Plot received PBCH constellation after equalization
figure;
plot(pbchEq,'o');
title('Equalized PBCH constellation');
m = max(abs([real(pbchEq(:)); imag(pbchEq(:))])) * 1.1;
axis([-m m -m m]);

% PBCH demodulation
pbchBits = nrPBCHDecode(pbchEq,ncellid,v,nest);

% Calculate RMS PBCH EVM
pbchRef = nrPBCH(pbchBits<0,ncellid,v);
evm = comm.EVM;
evm_rms = evm(pbchRef,pbchEq);

% Display calculated EVM
disp(['RMS PBCH EVM = ' num2str(evm_rms,'%0.3f') '%']);

% Apply CSI
pbchBits = pbchBits .* csi;

% Perform BCH decoding including rate recovery, polar decoding and CRC
% decoding. PBCH descrambling and separation of the BCH transport block
% bits 'trblk' from 8 additional payload bits A...A+7 is also performed:
%   A ... A+3: 4 LSBs of System Frame Number
%         A+4: half frame number
% A+5 ... A+7: for L=64, 3 MSBs of the SS/PBCH block index
%              for L=4 or 8, A+5 is the MSB of the subcarrier offset k_SSB
polarListLength = 8;
[~,err,trblk,sfn4lsb,nHalfFrame,msbidxoffset] = ...
    nrBCHDecode(pbchBits,polarListLength,burstInfo.L,ncellid);

% Use 'msbidxoffset' value to set bits of 'k_SSB' or 'ssbIndex', depending
% on the number of SS/PBCH blocks in the burst
if (burstInfo.L==64)
    ssbIndex = ssbIndex + (bi2de(msbidxoffset.','left-msb') * 8);
    k_SSB = 0;
else
    k_SSB = msbidxoffset * 16;
end

% Display the BCH CRC
disp(['BCH CRC = ' num2str(err)]);

% Displaying the SSB index
disp(['SSB index = ' num2str(ssbIndex)]);

% Create set of subcarrier spacings signaled by the 7th bit of the decoded
% MIB, the set is different for FR1 (L=4 or 8) and FR2 (L=64)
if (burstInfo.L==64)
    commonSCSs = [60 120];
else
    commonSCSs = [15 30];
end

% Create a structure of MIB fields from the decoded MIB bits. The BCH
% transport block 'trblk' is the RRC message BCCH-BCH-Message, consisting
% of a leading 0 bit then 23 bits corresponding to the MIB
mib.NFrame = bi2de([trblk(2:7); sfn4lsb] .','left-msb');
mib.SubcarrierSpacingCommon = commonSCSs(trblk(8) + 1);
mib.k_SSB = k_SSB + bi2de(trblk(9:12).','left-msb');
mib.DMRSTypeAPosition = 2 + trblk(13);
mib.PDCCHConfigSIB1 = bi2de(trblk(14:21).','left-msb');
mib.CellBarred = trblk(22);
mib.IntraFreqReselection = trblk(23);

% Display the MIB structure
disp(mib);

function frequencyOffset = frequencyOffsetSSB(waveform,ofdmInfo)

    % Get 'Lsym', the number of time domain samples in an OFDM symbol,
    % which is the sum of 'Lcp' and 'Lu', the number of cyclic prefix
    % samples and useful samples respectively
    Lcp = ofdmInfo.CyclicPrefixLengths(2);
    Lu = ofdmInfo.Nfft;
    Lsym = Lcp + Lu;

    % Multiply the waveform by itself delayed by Lu samples and conjugated
    delayed = [zeros(Lu,size(waveform,2)); waveform(1:end-Lu,:)];
    cpProduct = waveform .* conj(delayed);

    % Apply a moving sum filter with a window size equal to the CP length
    cpXCorr = filter(ones([Lcp 1]),1,cpProduct);

    % Moving sum over 4 OFDM symbols (i.e. the size of the SS block)
    y = cpXCorr;
    cpXCorrDelayed = cpXCorr;
    for k = 1:3
        cpXCorrDelayed = [zeros(Lsym,size(waveform,2)); cpXCorrDelayed(1:end-Lsym,:)];
        y = y + cpXCorrDelayed;
    end

    % Extract the correlation peak, average over the receive antennas,
    % then compute the phase and corresponding frequency offset
    cpCorrIndex = (Lsym * 4) + 1;
    frequencyOffset = ofdmInfo.SubcarrierSpacing * 1e3 * angle(mean(y(cpCorrIndex,:))) / (2*pi);

end