if ~exist('dvbs2xLDPCParityMatrices.mat','file')
    if ~exist('s2xLDPCParityMatrices.zip','file')
        url = 'https://ssd.mathworks.com/supportfiles/spc/satcom/DVB/s2xLDPCParityMatrices.zip';
        websave('s2xLDPCParityMatrices.zip',url);
        unzip('s2xLDPCParityMatrices.zip');
    end
    addpath('s2xLDPCParityMatrices');
end

cfgDVBS2X.StreamFormat = "TS";
cfgDVBS2X.PLSDecimalCode = 129;                   
cfgDVBS2X.CanonicalMODCODName = "QPSK 2/9";
cfgDVBS2X.DFL = 14128;
cfgDVBS2X.RolloffFactor = 0.35;
cfgDVBS2X.SamplesPerSymbol = 2

simParams.sps = cfgDVBS2X.SamplesPerSymbol;    % Samples per symbol
simParams.numFrames = 20;                      % Number of frames to be processed
simParams.chanBW = 36e6;                       % Channel bandwidth in Hertz
simParams.cfo = 3e6;                           % Carrier frequency offset in Hertz due
                                               % to oscillator instabilities
simParams.dopplerRate = 3e3;                   % Doppler rate in Hertz/sec
simParams.peakDoppler = 20e3;                  % Peak doppler shift due to receiver motion
simParams.sco = 10;                            % Sampling clock offset in parts per
                                               % million
simParams.phNoiseLevel = "Level-1";            % Phase noise level
simParams.EsNodB = 5;                          % Energy per symbol to noise ratio

[data,txOut,rxIn,phyConfig,rxParams] = ...
    HelperDVBS2XVLSNRRxInputGenerate(cfgDVBS2X,simParams);
disp(phyConfig)

% Received signal constellation plot
rxConst = comm.ConstellationDiagram('Title','Received data', ...
    'XLimits',[-1 1],'YLimits',[-1 1], ...
    'ShowReferenceConstellation',false, ...
    'SamplesPerSymbol',simParams.sps);
rxConst(rxIn(1:rxParams.plFrameSize*simParams.sps))

% Transmitted and received signal spectrum visualization
Rsymb = simParams.chanBW/(1 + cfgDVBS2X.RolloffFactor);
Fsamp = Rsymb*simParams.sps;
specAn = dsp.SpectrumAnalyzer('SampleRate',Fsamp, ...
    'ChannelNames',{'Transmitted waveform','Received waveform'}, ...
    'ShowLegend',true);
specAn([txOut,rxIn(1:length(txOut))]);

rxParams.carrSyncLoopBW = 1e-4;          % Coarse frequency estimator loop bandwidth
                                         % normalized by symbol rate
rxParams.symbSyncLoopBW = 1e-4;          % Symbol timing synchronizer loop bandwidth
                                         % normalized by symbol rate  
rxParams.initialTimeFreqSync = 5;        % Number of frames required for coarse frequency
                                         % and symbol timing error convergence
rxParams.fineFreqLock = 4;               % Number of frames required for fine
                                         % frequency estimation
rxParams.NeedSmoothening = false;        % Smoothen the phase estimate

% Total frames used for overall synchronization
rxParams.totalSyncFrames = rxParams.initialTimeFreqSync + rxParams.fineFreqLock;

% Create coarse frequency synchronization System object by using
% HelperDVBS2XVLSNRCoarseFreqSynchronizer helper object
freqSync = HelperDVBS2XVLSNRCoarseFreqSynchronizer('SamplesPerSymbol',simParams.sps, ...
    'NormalizedLoopBandwidth',rxParams.carrSyncLoopBW);

% Create symbol timing synchronization System object by using
% comm.SymbolSynchronizer object
symSync = comm.SymbolSynchronizer('TimingErrorDetector','Gardner (non-data-aided)', ...
    'NormalizedLoopBandwidth',rxParams.symbSyncLoopBW);

% Create matched filter System object by using
% comm.RaisedCosineReceiveFilter object
if simParams.sps == 2
    decFac = 1;
else
    decFac = simParams.sps/(simParams.sps/2);
end
rxFilter = comm.RaisedCosineReceiveFilter( ...
    'RolloffFactor',0.35, ...
     'InputSamplesPerSymbol',simParams.sps,'DecimationFactor',decFac);
b = rxFilter.coeffs;
rxFilter.Gain = sum(b.Numerator);

% Initialize error computing parameters
[numFramesLost,pktsErr,bitsErr,pktsRec] = deal(0);

% Initialize data indexing variables
stIdx = 0;
dataSize = rxParams.inputFrameSize;
plFrameSize = rxParams.plFrameSize;
isLastFrame = false;
rxParams.fineFreqCorrVal = zeros(rxParams.fineFreqLock,1);
[formatIdx,formatIdxTemp] = deal(1);
vlsnrSyncStIdx = 93;
payloadStIdx = 899;
vlSNRFrameLen = plFrameSize - vlsnrSyncStIdx + 1;

while stIdx < length(rxIn)

    % Use one DVB-S2X PL frame for each iteration.
    endIdx = stIdx + rxParams.plFrameSize*simParams.sps;

    % In the last iteration, all the remaining samples in the received
    % waveform are considered.
    isLastFrame = endIdx > length(rxIn);
    endIdx(isLastFrame) = length(rxIn);
    rxData = rxIn(stIdx+1:endIdx);

    % After coarse frequency offset loop is converged, the FLL works with
    % previous frequency estimate.
    if rxParams.frameCount < rxParams.initialTimeFreqSync
        coarseFreqLock = false;
    else
        coarseFreqLock = true;
    end

    % Retrieve the last frame samples.
    if isLastFrame
        resSampCnt = plFrameSize*rxParams.sps - length(rxData);
        % Inadequate number of samples to fill last frame
        syncIn = [rxData; zeros(resSampCnt,1)];
    else
        syncIn = rxData;
    end

    % Apply coarse frequency offset compensation.
    [coarseFreqSyncOut,phEst] = freqSync(syncIn,coarseFreqLock);

    % Perform matched filtering and downsample the signal to 2 samples per
    % symbol.
    filtOut = rxFilter(coarseFreqSyncOut);

    % Apply symbol timing synchronization.
    symSyncOut = symSync(filtOut);

    % Apply frame synchronization and identify the MODCOD format. VL-SNR
    % sync frame is detected. PL header preceding the VL-SNR header is
    % ignored.
    if rxParams.frameCount > rxParams.initialTimeFreqSync && ~isLastFrame
        [~,rxParams.syncIndex,formatIdx] = ...
            HelperDVBS2XVLSNRFrameSync(symSyncOut,rxParams.SegLength);
        % MODCOD format identification failure verification
        formatFail = formatIdxTemp ~= rxParams.refFormat;
        if formatFail && ~isLastFrame
            % Update the counters, state variables, and buffers
            stIdx = endIdx;
            rxParams.frameCount = rxParams.frameCount + 1;
            formatIdxTemp = formatIdx;
            rxParams.cfBuffer = symSyncOut(rxParams.syncIndex:end);
            fprintf('%s\n','MODCOD format detection failed')
            continue;
        else
            fprintf('%s\n','MODCOD format detection passed')
            [setNum,phyParams] = getVLSNRParams(formatIdxTemp);
        end
    end

    % The PL frame start index lies somewhere in the middle of the data
    % being processed. From fine frequency estimation onwards, the
    % processing happens as a PL frame. A buffer is used to store symbols
    % required to fill one PL frame. PL frame is considered from
    % the start of VL-SNR header, precisely from the start of the 896 bit
    % length Walsh Hadamard (WH) sequence. 90 represents the PL header
    % length, and 2 accounts for the two zeros appended before the WH sequence.
    fineFreqIn = [rxParams.cfBuffer;...
        symSyncOut(1:vlSNRFrameLen-length(rxParams.cfBuffer))];

    % Estimate the fine frequency error by using the HelperDVBS2FineFreqEst
    % helper function.
    % Add 1 to the conditional check because the buffer used to get one PL
    % frame introduces a delay of one to the loop count.
    if (rxParams.frameCount > rxParams.initialTimeFreqSync + 1)
        % Extract the payload by removing header
        payload = fineFreqIn(payloadStIdx:end); 

        % Get the correlation estimate from the regular 36 length pilot blocks.
        est1 = HelperDVBS2FineFreqEst( ...
            payload(rxParams.regPilotInd),rxParams.regNumPilotBlks, ...
            rxParams.regPilotSeq,rxParams.fineFreqState,36,rxParams.NumLags);

        % Get the correlation estimate from the VL-SNR extra pilot blocks.
        % vlSNRPilotBlk1Params contains the VL-SNR type 1 pilot block
        % length and number of blocks in one frame.
        Lp = rxParams.vlSNRPilotBlk1Params(1);
        est2 = HelperDVBS2FineFreqEst( ...
            payload(rxParams.vlSNRPilotInd1),rxParams.vlSNRPilotBlk1Params(2), ...
            rxParams.vlSNRPilotSeq1,rxParams.fineFreqState,Lp,rxParams.NumLags);

        % vlSNRPilotBlk2Params contains the VL-SNR type 2 pilot block
        % length and number of blocks in one frame.
        Lp = rxParams.vlSNRPilotBlk2Params(1);
        est3 = HelperDVBS2FineFreqEst( ...
            payload(rxParams.vlSNRPilotInd2),rxParams.vlSNRPilotBlk2Params(2), ...
            rxParams.vlSNRPilotSeq2,rxParams.fineFreqState,Lp,rxParams.NumLags);

        estVal = est1 + est2 + est3;
        % Use the correlation values calculated over pilot fields spanning over
        % multiple frames to calculate the fine frequency error estimate.
        % The estimation uses a sliding window technique.
        rxParams.fineFreqCorrVal = [rxParams.fineFreqCorrVal(2:end);estVal];
    end

    if rxParams.frameCount >= rxParams.totalSyncFrames
        fineFreqLock = true;
    else
        fineFreqLock = false;
    end

    if fineFreqLock
        freqEst = angle(sum(rxParams.fineFreqCorrVal))/(pi*(rxParams.NumLags+1));
        ind = (rxParams.frameCount-2)*plFrameSize:(rxParams.frameCount-1)*plFrameSize-1;
        phErr = exp(-1j*2*pi*freqEst*ind).';
        % 92 accounts for the PL header (90), and 2 for zeros appended before the WH
        % sequence.
        fineFreqOut = fineFreqIn(1:vlSNRFrameLen).*phErr(vlsnrSyncStIdx:end);
        % 898 accounts for the 896 length WH sequence, and 2 zeros padded to the
        % header.
        rxPilots = fineFreqOut(rxParams.pilotInd+payloadStIdx-1);
        phErrEst = HelperDVBS2PhaseEst(rxPilots,rxParams.pilotSeq, ...
            rxParams.phErrState,rxParams.IsVLSNR,setNum,rxParams.Alpha);
        if rxParams.NeedSmoothening
            phErrEst =  smoothenEstimate(phErrEst);
        end
        phaseCompOut = HelperDVBS2PhaseCompensate(fineFreqOut(payloadStIdx:end), ...
            phErrEst,rxParams.pilotInd,setNum,rxParams.IsVLSNR);
    end

    % Recover the input bit stream.
    if rxParams.frameCount >= rxParams.totalSyncFrames
        isValid = true;
        syncOut = phaseCompOut;
    else
        isValid = false;
        syncOut = [];
    end

    % Update the buffers and counters.
    rxParams.cfBuffer = symSyncOut(rxParams.syncIndex:end);

    if isValid  % Data valid signal
        % Estimate noise variance by using
        % HelperDVBS2NoiseVarEstimate helper function.
        nVar = HelperDVBS2NoiseVarEstimate(syncOut,rxParams.pilotInd, ...
            rxParams.pilotSeq,false);
        % Recover the BB frame by using HelperDVBS2XBBFrameRecover
        % helper function.
        rxBBFrame = HelperDVBS2XBBFrameRecover(syncOut,phyParams, ...
            rxParams.plScramblingIndex,rxParams.regNumPilotBlks,nVar,true,setNum);
        % Recover the input bit stream by using
        % HelperDVBS2StreamRecover helper function.
        if strcmpi(cfgDVBS2X.StreamFormat,'GS') && ~rxParams.UPL
            [decBits,isFrameLost] = HelperDVBS2StreamRecover(rxBBFrame);
            if ~isFrameLost && length(decBits) ~= dataSize
                isFrameLost = true;
            end
        else
            [decBits,isFrameLost,pktCRC] = HelperDVBS2StreamRecover(rxBBFrame);
            if ~isFrameLost && length(decBits) ~= dataSize
                isFrameLost = true;
                pktCRC = zeros(0,1,'logical');
            end
            % Compute the PER for TS or GS packetized
            % mode.
            pktsErr = pktsErr + numel(pktCRC) - sum(pktCRC);
            pktsRec = pktsRec + numel(pktCRC);
        end
        if ~isFrameLost
            ts = sprintf('%s','BB header decoding passed.');
        else
            ts = sprintf('%s','BB header decoding failed.');
        end
        % Compute the number of frames lost. CRC failure of the
        % baseband header is considered a frame loss.
        numFramesLost = isFrameLost + numFramesLost;
        fprintf('%s(Number of frames lost = %1d)\n',ts,numFramesLost)
        % Compute the bits in error.
        if ~isFrameLost
            dataInd = (rxParams.frameCount-2)*dataSize+1:(rxParams.frameCount-1)*dataSize;
            errs = sum(data(dataInd) ~= decBits);
            bitsErr = bitsErr + errs;
        end
    end

    stIdx = endIdx;
    rxParams.frameCount = rxParams.frameCount + 1;
    formatIdxTemp = formatIdx;
end

% Synchronized data constellation plot
syncConst = comm.ConstellationDiagram('Title','Synchronized data', ...
    'XLimits',[-1.7 1.7],'YLimits',[-1.7 1.7], ...
    'ShowReferenceConstellation',false);
syncConst(syncOut)

pause(0.5)
% Error metrics display
% For GS continuous streams
if strcmpi(cfgDVBS2X.StreamFormat,'GS') && ~rxParams.UPL
    if (simParams.numFrames-rxParams.totalSyncFrames == numFramesLost)
        fprintf("All frames are lost. No bits are retrieved from BB frames.")
    else
        numFramesRec = simParams.numFrames - rxParams.totalSyncFrames - numFramesLost;
        ber = bitsErr/(numFramesRec*dataSize);
        fprintf('BER           : %1.2e\n',ber)
    end
else
    % For GS and TS packetized streams
    if pktsRec == 0
        fprintf("All frames are lost. No packets are retrieved from BB frames.")
    else
        if strcmpi(cfgDVBS2X.StreamFormat,'TS')
            pktLen = 1504;
        else
            pktLen = cfgDVBS2X.UPL;      % UP length including sync byte
        end
        ber = bitsErr/(pktsRec*pktLen);
        per = pktsErr/pktsRec;
        fprintf('PER: %1.2e\n',per)
        fprintf('BER: %1.2e\n',ber)
    end
end

% [coarseFreqSyncOut,phEst] = freqSync(syncIn,false);
% Frequency offset estimate normalized by symbol rate
% freqOffEst = diff(phEst(1:simParams.sps:end))/(2*pi);
% plot(freqOffEst)
% Rsym = simParams.chanBW/(1+cfgDVBS2X.RolloffFactor);
% actFreqOff = (simParams.cfo + simParams.peakDoppler)/Rsym;
% resCoarseCFO = abs(actFreqOff-freqOffEst(end));

% fineFreqEst = angle(sum(rxParams.fineFreqCorrVal))/(pi*(rxParams.NumLags+1));
% resFineCFO = abs(actFreqOff-freqOffEst(end)-fineFreqEst);

function [setNum,phyParams] = getVLSNRParams(formatIdx)

tableVLSNR = [4 64800 2/9 0; ...
    2 32400 1/5 0; ...
    2 32400 11/45 0; ...
    2 32400 1/3 0; ...
    2 16200 1/5 1; ...
    2 16200 11/45 1; ...
    2 16200 1/5 0; ...
    2 16200 4/15 0; ...
    2 16200 1/3 0];

params = tableVLSNR(formatIdx,:);
phyParams.ModOrder = params(1);
phyParams.FECFrameLen = params(2);
phyParams.CodeRate = params(3);
[n, d] = rat(phyParams.CodeRate);
phyParams.CodeIDF = [sprintf('%0.0f',n) '/' sprintf('%0.0f',d)];
phyParams.HasPilots = true;
phyParams.HasSpread = params(4);

if formatIdx > 6
    setNum = 2;
else
    setNum = 1;
end
end
% Smoothen the phase estimate using moving average filter
function newEst = smoothenEstimate(est)
errDiff = diff(est(2:end));
thres = -2*sign(mean(errDiff));
width = 5;
if thres > 0
   index = find(errDiff > thres); 
else
  index = find(errDiff < thres);  
end
if ~isempty(index)
    for k = 1:length(index)
        est(index(k)+2:end) = est(index(k)+2:end)-sign(thres)*2*pi;
    end
end
temp = est(2:end);
n = length(temp);
c = filter(ones(width,1)/width,1,temp);
cbegin = cumsum(temp(1:width-2));
cbegin = cbegin(1:2:end)./(1:2:(width-2))';
cend = cumsum(temp(n:-1:n-width+3));
cend = cend(end:-2:1)./(width-2:-2:1)';
c = [cbegin;c(width:end);cend];
newEst = [est(1);c];
end