%% NR Synchronization Procedures
%
% This example demonstrates how to construct a waveform containing a
% synchronization signal burst (SS burst), pass the waveform through a
% fading channel with AWGN, and then blindly synchronize to the received
% waveform to decode the master information block (MIB). 

% Copyright 2018-2020 The MathWorks, Inc.

%% Introduction
%
% Before a User Equipment (UE) device can communicate with the network, it 
% must perform cell search and selection procedures and obtain initial
% system information. The first few steps in that process are acquiring
% frame synchronization, finding out the cell identity and decoding the
% Master Information Block (MIB). This example demonstrates how those steps
% can be performed with 5G Toolbox(TM).
%
% The figure below shows the main steps in the processing chain.
%
% <<../SynchProceduresExampleProcessingChain.png>>
%
% * _Burst generation_: A structure |txBurst| is created which configures
% an SS burst, and the function |hSSBurst| is used to create an OFDM
% resource grid containing the SS burst.
% * _Beam sweep_: The OFDM resource grid for the SS burst is beamformed
% onto a set of physical transmission antennas, with each SS/PBCH block in
% the burst having a different beamforming vector.
% * _Propagation channel_: The transmitted waveform is passed through a TDL
% propagation channel model.
% * _AWGN_: Additive White Gaussian Noise is applied to the receive
% antennas.
% * _Receiver_: Various synchronization and demodulation processes are
% applied to the received waveform in order to establish the cell identity
% and to decode the Master Information Block (MIB).
%
% The figure below shows the processing steps inside the receiver.
%
% <<../SynchProceduresExampleReceiver.png>>
%
% These processing steps are explained in detail below.

%% Burst Configuration
% A structure |txBurst| is created which configures an SS burst, including
% configuration of the pattern of SS/PBCH blocks within the burst and the
% content of the Master Information Block (MIB).

% 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;

%% Channel Configuration
% A Tapped Delay Line (TDL) propagation channel |channel| is configured, as
% well as the SNR for AWGN added at the receiver.

% 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;

%% Burst Generation
% The function |hSSBurst| is used to create an OFDM resource grid
% containing the SS burst. The resource grid is sized such that the
% corresponding OFDM modulated waveform has a sample rate equal to that
% specified by |txBurst.SampleRate|, which allows it to be easily added to
% a waveform carrying PDCCH and PDSCH. Time domain combining of the 
% waveforms is required in the case that the SS/PBCH block and PDCCH/PDSCH
% have different subcarrier spacings.

% Display burst 
disp(txBurst);

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

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

%% Beam Sweep
% The OFDM resource grid for the SS burst is beamformed onto a set of
% physical transmission antennas, with each SS/PBCH block in the burst
% having a different beamforming vector. The beamformed OFDM resource grid
% is then OFDM modulated to give a time domain waveform.

% 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);

%% Propagation Channel
% The transmitted waveform is passed through a TDL propagation channel
% model resulting in a received waveform for the configured number of
% antennas.

rxWaveform = channel(txWaveform);

%% AWGN
% Additive White Gaussian Noise is applied to the receive antennas.

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

%% Receiver Configuration
% In order to synchronize and demodulate the received waveform, the
% following information is needed:
%
% * The SS block pattern (Case A...E): the UE will know which block 
% patterns need to be searched based on the NR operating band, see 
% TS 38.104 Tables 5.4.3.3-1 and 5.4.3.3-2 [ <#18 1> ]. The SS block
% pattern determines the subcarrier spacing of the SS/PBCH blocks.
% * The waveform sample rate: required for OFDM demodulation of the 
% received waveform.
% * $L_{max}$, the number of SS/PBCH blocks in a burst: the UE will know
% the value of $L_{max}$ based on the SS block pattern and the NR operating
% band, TS 38.213 Section 4.1 [ <#18 2> ] describes the set of SS/PBCH
% blocks in a burst in each case. $L_{max}$ is used in calculation of
% parameters for the PBCH DM-RS sequences and PBCH descrambling. These
% parameters are a function of the SS/PBCH block index, see TS 38.211 
% Sections 7.3.3.1 and 7.4.1.4.1 [ <#18 2> ].

% 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);

%% PSS and Coarse Frequency Search
% PSS search is performed, which consists of correlating the received
% waveform (across all SS/PBCH blocks) with each of the three possible PSS
% sequences and extracting the strongest correlation peak. The SS/PBCH
% block with the strongest correlation peak indicates which beam in the
% beam sweep was most effective at directing the signal towards the
% receiver. The coarse frequency offset of the received waveform is also
% determined, by making various hypotheses about the frequency offset,
% correcting for each offset, and finding the offset that gives the
% strongest correlation.

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;

%% Fine Frequency Offset Estimation and Correction
% The timing offset for the strongest PSS sequence correlation can be used
% to synchronize the waveform in time. Then a fine frequency offset 
% estimate can be calculated by performing correlation between the cyclic 
% prefix of each OFDM symbol in the SSB, and the corresponding useful parts
% of the OFDM symbols. The phase of this correlation is proportional to the
% frequency offset in the waveform.

% 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,:);

%% SSS Search
% The timing of the PSS correlation peak is used to synchronize the
% waveform and OFDM demodulation is performed. The subcarriers associated
% with the SSS are extracted and correlated with each possible SSS
% sequence. The indices of the strongest PSS and SSS sequences are combined
% to give the physical layer cell identity, which is required for PBCH
% DM-RS and PBCH processing.

% 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;

%% PBCH DM-RS search
% In a process similar to SSS search, each possible PBCH DM-RS sequence is
% constructed and channel estimation and noise estimation is performed. The
% index of the PBCH DM-RS with the best SNR determines the LSBs of the
% SS/PBCH block index, required for PBCH scrambling initialization.

% 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)]);

%% Channel Estimation using PBCH DM-RS and SSS
% Now that the PBCH DM-RS sequence is known, a channel estimate for the
% SS/PBCH block can be created by estimating the channel in each PBCH DM-RS
% resource element location and interpolating across the SS/PBCH block. The
% SSS is also used to aid channel estimation. An estimate of the additive
% noise on the PBCH DM-RS / SSS is also performed.

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]);

%% PBCH Demodulation
% The subcarriers associated with the PBCH are extracted and the channel
% and noise estimates are used to perform MMSE equalization. The equalized
% PBCH symbols are then demodulated and descrambled to give bit estimates
% for the coded BCH block.

% 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') '%']);

%% BCH Decoding
% The BCH bit estimates are weighted with Channel State Information (CSI)
% from the MMSE equalizer then BCH decoding is performed, consisting of
% rate recovery, polar decoding, CRC decoding, descrambling and separating
% the 24 BCH transport block bits from the 8 additional timing-related
% payload bits.

% 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)]);

%% MIB Parsing
% Finally the 24 decoded BCH transport block bits are parsed into a
% structure which represents the MIB message fields. This includes
% reconstituting the 10-bit System Frame Number (SFN) |NFrame| from the 6
% MSBs in the MIB and the 4 LSBs in the PBCH payload bits. It also includes
% incorporating the MSB of the subcarrier offset |k_SSB| from the PBCH
% payload bits in the case of L=4 or 8 SS/PBCH blocks per burst.

% 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);

%% Appendix
% This example uses the following helper functions:
% 
% * <matlab:edit('hOFDMDemodulate.m') hOFDMDemodulate.m>
% * <matlab:edit('hOFDMInfo.m') hOFDMInfo.m>
% * <matlab:edit('hOFDMModulate.m') hOFDMModulate.m>
% * <matlab:edit('hSSBurst.m') hSSBurst.m>
% * <matlab:edit('hSSBurstInfo.m') hSSBurstInfo.m>

%% Selected Bibliography
% # 3GPP TS 38.104. "NR; Base Station (BS) radio transmission and
% reception (Release 15)." 3rd Generation Partnership Project; Technical
% Specification Group Radio Access Network.
% # 3GPP TS 38.213. "NR; Physical layer procedures for control (Release
% 15)." 3rd Generation Partnership Project; Technical Specification Group
% Radio Access Network.

%% Local functions

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
