End-to-End CCSDS High Photon Efficiency Telemetry Optical Link Simulation

cfgParams.ModOrder = 16;            % PPM modulation order
cfgParams.InfoSize = 7526;          % Information block size
cfgParams.NumRegisters = 18;        % Number of shift registers for Convolutional Interleaving
cfgParams.RegisterLengthStep = 420; % Additional registers in each step for Convolutional Interleaving
cfgParams.RepeatFactor = 2;         % Repetition factor
cfgParams.TFLength = 1024;          % Transfer frame length - 1024 bytes

numTF = 20;                    % Number of transfer frames
nsLen = 4;
ns = linspace(0.9,1.06,nsLen); % Average number of signal photons per pulsed slot
nb = 0.05;                     % Average number of noise photons per slot

% Initialize timing offset for the first iteration, to visualize the
% initialized offset with the calculated offset in the timing offset plot.
% From next iterations random offset is considered
slotOffset = 11;
ts = 0.45;                     % Normalized timing offset

% Initialize data for simulation
blockErrors = zeros(nsLen,1);
totalBlocks = zeros(nsLen,1);
maxNumBlocks = 60;
plotOffset = true;             % Flag to plot timing offset  
m = log2(cfgParams.ModOrder);
sps = 5*cfgParams.ModOrder/4;  % Slots per symbol

% Initialize waveform generator
hpeWavegen = ccsdsHPEWaveformGenerator(TransmissionType="telemetry", ...
    NumBytesInTransferFrame=cfgParams.TFLength, ...
    NumBitsInInformationBlock=cfgParams.InfoSize, ...
    PPMOrder=cfgParams.ModOrder, ...
    NumRows=cfgParams.NumRegisters, ...
    RegisterLengthIncrement=cfgParams.RegisterLengthStep, ...
    RepeatFactor=cfgParams.RepeatFactor);

% Initialize channel
dsocChan = dsocPoissonChannel(NumNoisePhotons=nb, ...
    TimingOffset=ts);

% Initialize PN Sequence
pnSequenceObj = comm.PNSequence(SamplesPerFrame=cfgParams.InfoSize, ...
    Polynomial=[1 1 0 1 0 1 0 0 1], ...
    InitialConditions=ones(8,1));
pnSequence = pnSequenceObj();
blockLen = 15120/m;
% Initialize convolutional deinterleaver
deintrlvrObj = comm.ConvolutionalDeinterleaver(NumRegisters=cfgParams.NumRegisters, ...
    RegisterLengthStep=cfgParams.RegisterLengthStep);

% Initialize frame synchronization object
ASM = [0 0 0 1 1 0 1 0 1 1 0 0 1 1 1, ...
    1 1 1 1 1 1 1 0 0 0 0 0 1 1 1 0 1]';
TFLengthInBits = cfgParams.TFLength*8;
fsObj = ccsdsTMFrameSynchronizer(SyncMarker=ASM, ...
      FrameLength=TFLengthInBits);

% Number of additional steps required for the last symbol to appear at the
% output in channel deinterleaving
numAdditionalSteps = cfgParams.RegisterLengthStep*cfgParams.NumRegisters*(cfgParams.NumRegisters-1);

for itr = 1:nsLen
    rng("default")
    dsocChan.NumSignalPhotons = ns(itr);
    disp("Running link simulation for ns = "+num2str(ns(itr))+"...")
    while totalBlocks(itr) < maxNumBlocks
        % Generate transfer frames with random data
        dataTF = randi([0 1],TFLengthInBits*numTF,1);

        % Generate CCSDS HPE TM Waveform
        txOut = hpeWavegen(dataTF);

        % Add slot offset and pass modulated data through poisson channel
        chanOutput = dsocChan([randi([0 1],slotOffset,1); txOut]);

        % Estimate timing offset, mean signal photons per pulsed slot and
        % background noise photons per slot. Plot maximum likelihood timing
        % offset estimate and log-likelihood function between each integer
        % interval for the first iteration. Highlight the final estimated
        % timing offset corresponding to the log-likelihood function peak
        % value.
        [estOffset,nsEst,nbEst] = HelperCCSDSHPETMTimingOffsetEstimator(chanOutput,m,plotOffset);
        plotOffset = false;                        % Clear the flag to prevent plotting in each iteration

        % Slot offset correction
        intOffset = floor(estOffset);              % Integer timing offset
        fracOffset = estOffset-intOffset;          % Fractional timing offset
        slotSyncOut = chanOutput(intOffset+1:end);

        % Compute LLRs
        zk0 = slotSyncOut*log(1 + (1-fracOffset)*(nsEst/nbEst));
        zk1 = [slotSyncOut(2:end); 0].*log(1 + fracOffset*(nsEst/nbEst));
        llr = zk0 + zk1;

        % Remove guard slots
        llrSym = buffer(llr,sps);
        ssIn = llrSym(1:cfgParams.ModOrder,:);

        % Perform repeat symbol synchronization and codeword synchronization
        ssOut = HelperCCSDSHPETMSymbolSynchronizer(ssIn,cfgParams.RepeatFactor);
        symLen = size(ssOut,2);

        % Perform channel deinterleaving
        chDeinterleaveOut = zeros(cfgParams.ModOrder,symLen-numAdditionalSteps);
        for idx = 1:cfgParams.ModOrder
            chOut = deintrlvrObj(ssOut(idx,:)');
            chDeinterleaveOut(idx,:) = chOut(1 + numAdditionalSteps:end);
            reset(deintrlvrObj);
        end
        release(deintrlvrObj);

        numInfoBlocks = size(chDeinterleaveOut,2)/blockLen;
        decData = zeros(cfgParams.InfoSize + 32 + 2,numInfoBlocks); % 32 CRC bits + 2 termination bits
        errBlock = zeros(1,numInfoBlocks);

        % SCPPM decoding
        for blockIdx = 1:numInfoBlocks
            sIdx = (blockIdx - 1)*blockLen + 1;
            eIdx = blockIdx*blockLen;
            code = chDeinterleaveOut(:,sIdx:eIdx);
            [decData(:,blockIdx),errBlock(blockIdx)] = ccsdsSCPPMDecode(code(:),hpeWavegen.CodeRate,m);
        end

        % Perform pseudo de-randomization on SMTFs
        psIn = repmat(pnSequence,1,numInfoBlocks);
        psOut = xor(psIn,decData(1:cfgParams.InfoSize,:));

        % Frame synchronization
        errBit = repmat(errBlock,size(psOut,1),1);
        [softFrames,sofPos] = fsObj(2*psOut(:)-1);
        outFrames = softFrames>0;

        % Obtain quality indicator and sequence indicator from the start of
        % frame position, sofPos. Set the sequence indicator, seqInd, to 0
        % if the transfer frame is the successor of the most recent
        % previous frame. Set the quality indicator, qualityInd, to 1 if
        % the recovered transfer frame is from the correctly decoded SCPPM
        % codewords.
        seqInd = [sofPos(1)~=length(ASM)+1, (sofPos(2:end)-sofPos(1:end-1)) ~= (TFLengthInBits+length(ASM))];
        qualityInd = zeros(length(sofPos),1);
        crcErr = errBit(:);
        for idx = 1:length(sofPos)
            qualityInd(idx) = sum(crcErr(sofPos(idx):sofPos(idx)+TFLengthInBits-1))==0;
        end

        % Update number of information block errors and total information
        % blocks for each iteration
        blockErrors(itr) = blockErrors(itr) + (numInfoBlocks - sum(~errBlock));
        totalBlocks(itr) = totalBlocks(itr) + numInfoBlocks;

        % Generate random timing offset for the next iteration
        slotOffset = randi(100,1,1);
        ts = rand;
        release(dsocChan)
        dsocChan.TimingOffset = ts;
        reset(hpeWavegen);
        reset(fsObj);
    end
end

figure
blockER = blockErrors./totalBlocks;
plot(ns, blockER,"*-")
xlim([0.9 2])
xlabel("Signal photons per slot, ns")
ylabel("Block error rate")
title("Performance for ModOrder = "+num2str(cfgParams.ModOrder)+", nb = "+num2str(nb))