LDPC Link Simulation Using GPU Processing

fprintf(...
    'DVBS-2 Digital Video Broadcast Standard Bit Error Rate Simulation\n\n');

fprintf(...
    'Performance comparison of CPU- and GPU- accelerated decoders.\n');

try
    % Query the GPU
    dev = gpuDevice;
    
    % Print out information about the GPU that was found
    fprintf(...
        'GPU detected (%s, %d multiprocessors, Compute Capability %s)\n',...
        dev.Name,dev.MultiprocessorCount,dev.ComputeCapability);
    
    % Include a GPU-based simulation.
    doGPU = true;
    
catch % #ok<CTCH>
    
    % The GPU is not supported or not present, or the Parallel Computing
    %Toolbox was not present and licensed. Consider a CPU-only simulation.
    
    inp = input(['***NOTE: GPU not detected. ', ...
        'Continue with CPU-only simulation? [Y]/N '],'s');
    if strcmpi(inp, 'y') || isempty(inp)
        doGPU = false;
    else
        return;
    end
end

% DVBS-2 System Parameters
subsystemType = 'QPSK 1/2'; % Constellation and LDPC code rate
EsNodB = 0.75;              % Energy per symbol to noise PSD ratio in dB
SNR = convertSNR(EsNodB,'esno',SamplesPerSymbol=1);
numFramesPerCall = 20;      % Number of frames per call to ldpcDecode
numCalls = 10;              % Total number of frames = numFramesPerCall * numCalls
maxNumLDPCIterations = 50;  % LDPC Decoder iterations

dvb = getParamsDVBS2Demo(subsystemType,EsNodB,maxNumLDPCIterations);

% Create and configure the BCH Encoder and Decoder, Modulator, Demodulator.
createSimObjDVBS2Demo;

% Construct an LDPC Encoder configuration object
encoderCfg = ldpcEncoderConfig(dvb.LDPCParityCheckMatrix);

% Construct an LDPC Decoder configuration object
decoderCfg = ldpcDecoderConfig(dvb.LDPCParityCheckMatrix);

% Create an ErrorRate object to analyze the differences in bit error rate
% between the CPU and GPU.

BER = comm.ErrorRate;

if doGPU
    architectures = 2;
else
    architectures = 1;
end

% Initialize run time results vectors
runtime = zeros(architectures,numCalls);
avgtime = zeros(1,architectures);

% Seed the random number generator used for the channel and message
% creation.  This will allow a fair BER comparison between CPU and GPU.
% Cache the original random stream to restore later.

original_rs = RandStream.getGlobalStream;
rs = RandStream.create('mrg32k3a','seed',25);
RandStream.setGlobalStream(rs);

% Loop for each processing unit - CPU and GPU
for ii = 1:architectures
    
    % Do some initial setup for the execution loop
    if (ii == 1)
        arch = 'CPU'; % Use CPU
    else
        arch = 'GPU'; % Use GPU
    end
    
    % Reset the Error Rate object
    reset(BER);
    
    % Reset the random stream
    reset(rs);
    
    % Notice to the user that DVBS-2 simulation is beginning.
    fprintf(['\nUsing ' arch '-based LDPC Decoder:\n']);
    dels = repmat('\b',1,fprintf('  Initializing ...'));

    % Main simulation loop. Run numCalls+1 times and ignore the first
    % call (which has initialization overhead) for the run time
    % calculation. Use the first run for the BER calculation.
    for rr = 1:(numCalls+1)
        
        % Start timer
        ts = tic;
        
        % Create input messages
        msg = zeros(encbch.MessageLength, numFramesPerCall);
        msg(1:dvb.NumInfoBitsPerCodeword, :) = ...
            logical(randi([0 1],dvb.NumInfoBitsPerCodeword,numFramesPerCall));
        
        % Transmit
        bchencOut = encbch(msg(:));
        ldpcencOut = ldpcEncode(reshape(bchencOut,[],numFramesPerCall),encoderCfg);
        xlvrOut = intrlv(ldpcencOut,dvb.InterleaveOrder);
        modOut = pskmod(xlvrOut,dvb.ModulationOrder,dvb.PhaseOffset,'InputType','bit');
        
        % Corrupt with noise
        chanOut = awgn(modOut,SNR);
        
        % Receive
        demodOut = pskdemod(chanOut,dvb.ModulationOrder,dvb.PhaseOffset,'OutputType','approxllr','NoiseVariance',dvb.NoiseVar);
        
        if strcmpi(arch,'GPU')
            dexlvrOut = deintrlv(gpuArray(demodOut),dvb.InterleaveOrder);            
            ldpcdecOut = gather(logical(ldpcDecode(dexlvrOut,decoderCfg,dvb.LDPCNumIterations,'DecisionType','hard','Termination','max','OutputFormat','info')));
        else
            dexlvrOut = deintrlv(demodOut,dvb.InterleaveOrder);
            ldpcdecOut = logical(ldpcDecode(dexlvrOut,decoderCfg,dvb.LDPCNumIterations,'DecisionType','hard','Termination','max','OutputFormat','info'));
        end

        bchdecOut = decbch(ldpcdecOut(:));

        % Compute BER
        % Calculate BER at output of LDPC, not BCH.
        ber = BER(logical(bchencOut),ldpcdecOut(:));
        
        % Stop timer
        runtime(ii, rr) = toc(ts);
        
        % Don't report the first call with the initialization overhead.
        if (rr > 1)
            fprintf(dels);
            newCharsToDelete = fprintf('  Decode call %d : %.2f sec', ...
                rr-1, runtime(ii,rr));
            dels = repmat('\b',1,newCharsToDelete);
        end
    end % end of running frames through the DVBS-2 system.
    
    
    % Report the run time results to the Command Window.
    fprintf(dels); % Delete the last line printed out.
    
    % Calculate the average run time. Don't include call 1 because it
    % includes some System object initialization time.
    avgtime(ii) = mean(runtime(ii,2:end));
    
    fprintf('  %d frames decoded, %.2f sec/frame\n',numCalls*numFramesPerCall,avgtime(ii));
    fprintf('  Bit error rate: %g \n',ber(1) );
    
end % architecture loop

% Reset the random stream to the cached object
RandStream.setGlobalStream(original_rs);

if ~doGPU
    fprintf('\n*** GPU not present ***\n\n');
else
    % Calculate system-wide speedup
    fprintf(['\nFull system simulation runs %.2f times faster using ' ...
        'ldpcDecode with a GPU.\n\n'],avgtime(1) / avgtime(2));
end