DVB-S2 Waveform Transmission and Reception Using Instrument Control Toolbox

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

% Set parameters for the spectrum analyzer
centerFrequency = 5.5e9;
measurementTime = 800e-3;
mechanicalAttenuation = 0;
startFrequency = 5.49e9;
stopFrequency = 5.51e9;
resolutionBandwidth = 220e3;
videoBandwidth = 220e3;

fecFrame = "normal";                            % DVB-S2 FECFrame type
modCod = 18;                                    % 32APSK 3/4
Fsamp = 4e6;                                    % Sampling rate in samples per second
Rsym = 1e6;                                     % Symbol rate in symbols per second
sps = Fsamp/Rsym;
numFrames = 50;
rxParams = getRxParams(fecFrame,modCod);
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.symbSyncLock              = 2;         % Number of frames required for symbol timing error convergence
rxParams.frameSyncLock             = 1;         % Number of frames required for frame synchronization
rxParams.sps = sps;
% Total frames taken for symbol timing and frame synchronization to happen
rxParams.initialTimeSync = rxParams.symbSyncLock + rxParams.frameSyncLock;
rxParams.initialTimeFreqSync = numFrames;
rxParams.totalSyncFrames = numFrames;

foundVISA = visadevlist;
resourceID = foundVISA(foundVISA.Model == "N9040B",:).ResourceName;
sigAnalyzerObj = visadev(resourceID);
sigAnalyzerObj.ByteOrder = "big-endian";
sigAnalyzerObj.Timeout = 20;

writeline(sigAnalyzerObj,"*RST");
instrumentInfo = writeread(sigAnalyzerObj,"*IDN?");
fprintf("Instrument identification information: %s",instrumentInfo);

% Set up signal analyzer mode to basic IQ mode
writeline(sigAnalyzerObj,":INSTrument:SELect BASIC");

% Set the center frequency
writeline(sigAnalyzerObj,strcat(":SENSe:FREQuency:CENTer ",num2str(centerFrequency)));

% Set the capture sample rate
writeline(sigAnalyzerObj,strcat(":SENSe:WAVeform:SRATe ",num2str(Fsamp)));

% Turn off averaging
writeline(sigAnalyzerObj,":SENSe:WAVeform:AVER OFF");

% Set the spectrum analyzer to take one single measurement after the trigger line goes high
writeline(sigAnalyzerObj,":INIT:CONT OFF");

% Set the trigger to external source 1 with positive slope triggering
writeline(sigAnalyzerObj,":TRIGger:WAVeform:SOURce IMMediate");
writeline(sigAnalyzerObj,":TRIGger:LINE:SLOPe POSitive");

% Set the time for which measurement needs to be made
writeline(sigAnalyzerObj,strcat(":WAVeform:SWE:TIME ",num2str(measurementTime)));

% Turn off electrical attenuation
writeline(sigAnalyzerObj,":SENSe:POWer:RF:EATTenuation:STATe OFF");

% Set the mechanical attenuation level
writeline(sigAnalyzerObj,strcat(":SENSe:POWer:RF:ATTenuation ",num2str(mechanicalAttenuation)));

% Turn IQ signal ranging to auto
writeline(sigAnalyzerObj,":SENSe:VOLTage:IQ:RANGe:AUTO ON");

% Set the endianness of returned data
writeline(sigAnalyzerObj,":FORMat:BORDer NORMal");

% Set the format of the returned data
writeline(sigAnalyzerObj,":FORMat:DATA REAL,64");

writeline(sigAnalyzerObj,"*TRG");
writeline(sigAnalyzerObj,":INITiate:WAVeform");
measureComplete = writeread(sigAnalyzerObj,"*OPC?");
writeline(sigAnalyzerObj,":READ:WAV0?");
data = readbinblock(sigAnalyzerObj,"double");

% Create the IQ vector
inphase = data(1:2:end);
quadrature = data(2:2:end);
rxWaveform = inphase+1i*quadrature;

writeline(sigAnalyzerObj,":FETCH:WAV1?");
signalSpec = readbinblock(sigAnalyzerObj,"double");
captureSampleRate = 1/signalSpec(1);
fprintf("Sample Rate (Hz) = %s",num2str(captureSampleRate));

fprintf("Number of points read = %s",num2str(signalSpec(4)));

fprintf("Max value of signal (dBm) = %s",num2str(signalSpec(6)));

fprintf("Min value of signal (dBm) = %s",num2str(signalSpec(7)));

rxWaveform = rxWaveform.';
specAn = spectrumAnalyzer(SampleRate=Fsamp);
specAn(rxWaveform)

timeFreqSync = HelperDVBS2TimeFreqSynchronizer( ...
    CarrSyncLoopBW=rxParams.carrSyncLoopBW, ...
    SymbSyncLoopBW=rxParams.symbSyncLoopBW, ...
    SamplesPerSymbol=rxParams.sps, ...
    DataFrameSize=rxParams.xFecFrameSize, ...
    SymbSyncTransitFrames=rxParams.symbSyncLock, ...
    FrameSyncAveragingFrames=rxParams.frameSyncLock);

agc = comm.AGC(AveragingLength=2000,AdaptationStepSize = 0.001);

dcBlock = dsp.DCBlocker(Algorithm="FIR",Length=100);

syncConst = comm.ConstellationDiagram(Title="Synchronized data", ...
XLimits=[-2 2], ...
YLimits=[-2 2], ...
ShowReferenceConstellation = false);
[numFramesLost,pktsErr,bitsErr,pktsRec,stIdx,numFreqEst] = deal(0);

% Initialize data indexing variables
plFrameSize = rxParams.plFrameSize;
symSyncOutLen = zeros(rxParams.initialTimeFreqSync,1);
freqCounter = 0;
isCoarseFreqLock = false;
cCFOEstMean = [];
isFineFreqLock = false;
cCFOEst = [];
rxParams.fineFreqEstTemp = zeros(numFrames-rxParams.initialTimeSync,1);
headerDataDisplayCounter = 0;
streamFormat = "packetized";

while stIdx < length(rxWaveform)

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

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

    % After coarse frequency offset loop is converged, the FLL works with a
    % reduced loop bandwidth.
    if rxParams.frameCount < rxParams.initialTimeFreqSync
        isCoarseFreqLock = false;
    else
        isCoarseFreqLock = true;
    end

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

    % Apply DC offset compensation, AGC, matched filtering, symbol timing synchronization, frame
    % synchronization, and coarse frequency offset compensation.
    syncIn = dcBlock(syncIn);
    syncIn = agc(syncIn)/sps;
    if rxParams.frameCount > 1
        [coarseFreqSyncOut,syncIndex,phEst] = timeFreqSync(syncIn,isCoarseFreqLock);
        coarseCFOEst = diff(phEst(1:sps:end)/(2*pi));
        if rxParams.frameCount <= rxParams.initialTimeSync
            symSyncOutLen = [symSyncOutLen;length(coarseFreqSyncOut)]; %#ok<*AGROW>
            if any(abs(diff(symSyncOutLen(1:rxParams.frameCount))) > 5)
                error(["Symbol timing synchronization failed. The loop will not " ...
                    "converge, resulting in frame recovery. Decrease the symbSyncLoopBW " ...
                    "parameter for proper loop convergence."]);
            end
        end
        if rxParams.frameCount > rxParams.initialTimeSync + 1
            cCFOEstMean = [cCFOEstMean;mean(coarseCFOEst)];
            cCFOEst = [cCFOEst;coarseCFOEst];
            if length(cCFOEstMean) > 2
                diffVal = diff(abs(cCFOEstMean));
                if all(diffVal(end-1:end) < 0.01) && ~isCoarseFreqLock
                    isCoarseFreqLock = true;
                    rxParams.initialTimeFreqSync = rxParams.frameCount;
                elseif isLastFrame && ~isCoarseFreqLock
                    fprintf("%s\n",["Coarse frequency error estimation failed. Try analyzing the cCFOEst " ...
                        "across frames and check carrier synchronization has converged. Try reducing the " ...
                        "carrSyncLoopBW parameters"]);
                end
            end
        end
        rxParams.syncIndex = syncIndex;
        % The PL frame start index lies somewhere in the middle of the frame being processed.
        % From fine frequency estimation, the processing happens as a PL frame.
        % A buffer stores the symbols required to fill one PL frame.
        if isLastFrame
            fineFreqIn = [rxParams.cfBuffer; coarseFreqSyncOut];
        elseif rxParams.frameCount > 1
            buffLen = rxParams.plFrameSize - length(rxParams.cfBuffer);
            if buffLen < length(coarseFreqSyncOut)
                fineFreqIn = [rxParams.cfBuffer; coarseFreqSyncOut(1:buffLen)];
            else
                fineFreqIn = [rxParams.cfBuffer; coarseFreqSyncOut];
            end
        end
        % 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) && ...
                ~isFineFreqLock
            rxParams.fineFreqCorrVal = HelperDVBS2FineFreqEst( ...
                fineFreqIn(rxParams.pilotInd),rxParams.numPilots, ...
                rxParams.refPilots,rxParams.fineFreqCorrVal);
            % Normalize the frequency estimate by the input symbol rate
            % freqEst = angle(R)/(pi*(N+1)), where N (18) is the number of elements
            % used to compute the mean of auto correlation (R) in
            % HelperDVBS2FineFreqEst.
            freqEst = angle(rxParams.fineFreqCorrVal)/(pi*(19));
            rxParams.fineFreqEstTemp(1:end-1) = rxParams.fineFreqEstTemp(2:end);
            rxParams.fineFreqEstTemp(end) = freqEst;
            numFreqEst = numFreqEst + 1;
            if numFreqEst >= 5
                cg = abs(diff(rxParams.fineFreqEstTemp));
                if all(cg(end-2:end) < 1e-4)
                    isFineFreqLock = true;
                    rxParams.totalSyncFrames = rxParams.frameCount;
fprintf("Estimated carrier frequency offset in Hz = %f\n",(coarseCFOEst(end) + freqEst).*Rsym);
                elseif isLastFrame && ~isFineFreqLock
                    fprintf("%s\n","Fine frequency error estimation failed. Try analyzing the pilot fields in the PL frame to debug the issue.")
                end
            end
        end
        if isFineFreqLock
            % Normalize the frequency estimate by the input symbol rate
            % freqEst = angle(R)/(pi*(N+1)), where N (18) is the number of elements
            % used to compute the mean of auto correlation (R) in
            % HelperDVBS2FineFreqEst.
            freqEst = angle(rxParams.fineFreqCorrVal)/(pi*(19));

            % Generate the symbol indices using frameCount and plFrameSize.
            % Subtract 2 from the rxParams.frameCount because the buffer used to get one
            % PL frame introduces a delay of one to the count.
            ind = (rxParams.frameCount-2)*plFrameSize:(rxParams.frameCount-1)*plFrameSize-1;
            phErr = exp(-1j*2*pi*freqEst*ind);
            fineFreqOut = fineFreqIn.*phErr(:);

            % Estimate the phase error estimation by using the HelperDVBS2PhaseEst
            % helper function.
            [phEstRes,rxParams.prevPhaseEst] = HelperDVBS2PhaseEst( ...
                fineFreqOut,rxParams.refPilots,rxParams.prevPhaseEst,rxParams.pilotInd);

            % Compensate for the residual frequency and phase offset by using
            % the HelperDVBS2PhaseCompensate helper function.
            if rxParams.frameCount >= rxParams.totalSyncFrames + 2
                phaseCompOut = HelperDVBS2PhaseCompensate(rxParams.ffBuffer, ...
                    rxParams.pilotEst,rxParams.pilotInd,phEstRes(2));
            end

            rxParams.ffBuffer = fineFreqOut;
            rxParams.pilotEst = phEstRes;

            % Perform phase compensation on the data portion by
            % interpolating the phase estimates computed on consecutive pilot
            % blocks. The second phase estimate is not available for the data
            % portion after the last pilot block in the last frame. Therefore,
            % the slope of phase estimates computed on all pilot blocks in the
            % last frame is extrapolated and used to compensate for the phase
            % error on the final data portion.
            if isLastFrame
                pilotBlkLen = 36;                                               % Symbols
                pilotBlkFreq = 1476;                                            % Symbols
                avgSlope = mean(diff(phEstRes(2:end)));
                chunkLen = rxParams.plFrameSize - rxParams.pilotInd(end) + ...
                    rxParams.pilotInd(pilotBlkLen);
                estEndPh = phEstRes(end) + avgSlope*chunkLen/pilotBlkFreq;
                phaseCompOut1 = HelperDVBS2PhaseCompensate(rxParams.ffBuffer, ...
                    rxParams.pilotEst,rxParams.pilotInd,estEndPh);
            end
        end

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

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

        if isValid  % Data valid signal

            % Decode the PL header by using the dvbsPLHeaderRecover
            % function
            syncOut = syncOut/sqrt(mean(abs(syncOut).^2));
            rxPLHeader = syncOut(1:90);
            phyParams = dvbsPLHeaderRecover(rxPLHeader,Mode="DVB-S2/S2X regular");
            M = phyParams.ModulationOrder;
            R = eval(phyParams.LDPCCodeIdentifier);
            fecFrame = phyParams.FECFrameLength;
            pilotStat = phyParams.HasPilots;
            xFECFrameLen = fecFrame/log2(M);
            % Validate the decoded PL header.
            if M ~= rxParams.modOrder || abs(R-rxParams.codeRate) > 1e-3 || ...
                    fecFrame ~= rxParams.cwLen || ~pilotStat
                fprintf("Pl frame number %d , %s\n",rxParams.frameCount, ...
                    "PL header decoding failed")
            else                                                                % Demodulation and decoding
                if ~headerDataDisplayCounter
fprintf("%s\n","Physical layer header detection summary")
detPLParams = generatePLParamsTable(M,R,fecFrame,pilotStat)
                    headerDataDisplayCounter = headerDataDisplayCounter + 1;
                end
                rxFrame = syncOut(1:plFrameSize);
                % Synchronized data constellation plot
syncConst(rxFrame)
                % Estimate noise variance by using
                % HelperDVBS2NoiseVarEstimate helper function.
                nVar = HelperDVBS2NoiseVarEstimate(rxFrame,rxParams.pilotInd,...
                    rxParams.refPilots,false);
                snr_dB = 10*log10(1/nVar);
fprintf("PL frame number = %d\n",rxParams.frameCount)
fprintf("Estimated SNR in dB = %f\n",snr_dB)
                % Recover the input bit stream by using
                % dvbs2BitRecover function
                [decBitsTemp,isFrameLost,pktCRC]  = dvbs2BitRecover(rxFrame,nVar);
                decBits = decBitsTemp{:};
                if ~isLastFrame
                    pktsErr = pktsErr + numel(pktCRC{:}) - sum(pktCRC{:});
                    pktsRec = pktsRec + numel(pktCRC{:});
                end
                if ~isFrameLost
fprintf("%s\n","BB header decoding passed")
                    if isempty(pktCRC{1})
                        streamFormat = "continuous";
                    else
                        streamFormat = "packetized";
                    end
                else
                    fprintf("%s\n","BB header decoding failed")
                end
                % Compute the number of frames lost. Consider CRC failure of the baseband header
                % as a frame loss.
                numFramesLost = isFrameLost + numFramesLost;
fprintf("Total number of frames lost = %1d\n",numFramesLost)
            end

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

% Error metrics display
% For GS and TS packetized streams
if strcmp(streamFormat,"packetized")
    if pktsRec == 0
        fprintf("All frames are lost. No packets are retrieved from BB frames.")
    else
        per = pktsErr/pktsRec;
fprintf("PER: %1.2e\n",per)
    end
end

function rxParams = getRxParams(fecFrame,modCod)
% Receiver parameters generation

[modOrder,codeRate,cwLen] = satcom.internal.dvbs.getS2PHYParams(modCod,fecFrame);

dataLen = cwLen/log2(modOrder);

% Pilot sequence and indices generation
slotLen = 90;
pilotBlkFreq = 16;                                                    % In slots
numPilotBlks = floor(dataLen/(slotLen*pilotBlkFreq));
if floor(dataLen/(slotLen*16)) == dataLen/(slotLen*pilotBlkFreq)
    numPilotBlks = numPilotBlks - 1;
end
pilotLen = numPilotBlks*36;                                           % One pilot block contains 36 pilot symbols
frameSize = dataLen + pilotLen + slotLen;
plScrambIntSeq = satcom.internal.dvbs.plScramblingIntegerSequence(0);
cMap = [1 1j -1 -1j].';
cSeq = cMap(plScrambIntSeq+1);
[~, pilotInd] = satcom.internal.dvbs.pilotBlock(numPilotBlks);

rxParams.plFrameSize = frameSize;
rxParams.xFecFrameSize = dataLen;
rxParams.modOrder = modOrder;
rxParams.codeRate = codeRate;
rxParams.cwLen = cwLen;
rxParams.frameCount = 1;
rxParams.numPilots = numPilotBlks;
rxParams.pilotInd = pilotInd + slotLen;
rxParams.refPilots = (1+1j)/sqrt(2).*cSeq(pilotInd);
rxParams.cfBuffer = [];
rxParams.ffBuffer = complex(zeros(frameSize, 1));
rxParams.pilotPhEst = zeros(numPilotBlks+1, 1);
[rxParams.prevPhaseEst, rxParams.fineFreqCorrVal] = deal(0);
rxParams.syncIndex = 1;
end
function tbl =  generatePLParamsTable(M,R,fecFrame,pilotStat)
modScheme = ["QPSK" "8PSK" "16APSK" "32APSK"];
Modulation = modScheme(log2(M)-1);
[n,d] = rat(R);
CodeRate = strcat(string(n),"/",string(d));
if fecFrame == 64800
    FECFrame = "normal";
else
    FECFrame = "short";
end
if pilotStat
    PilotStatus = "true";
else
    PilotStatus = "false";
end
tbl = table(Modulation,CodeRate,FECFrame,PilotStatus);
tbl = renamevars(tbl,["Modulation","CodeRate","FECFrame", "PilotStatus"],["Modulation","Code Rate","FEC Frame", "Pilot Status"]);
end