% Check if the 'Communications Toolbox Library for the Bluetooth Protocol'
% support package is installed or not.
commSupportPackageCheck('BLUETOOTH');

EbNo = 2:4:10;                       % Eb/No in dB
sps = 4;                             % Samples per symbol, must be greater than 1
dataLen = 42;                        % Data length in bytes, includes header, payload and CRC
simMode = {'LE1M','LE2M','LE500K','LE125K'}; % PHY modes considered for the simulation

maxNumErrors = 100; % Maximum number of bit errors at an Eb/No point
maxNumPackets = 10; % Maximum number of packets at an Eb/No point

numMode = numel(simMode);          % Number of modes
ber = zeros(numMode,length(EbNo)); % Pre-allocate to store BER results
per = zeros(numMode,length(EbNo)); % Pre-allocate to store PER results
bitsPerByte = 8;                   % Number of bits per byte

% Fixed access address Ideally, this access address value should meet the
% requirements specified in Section 2.1.2 of the Bluetooth specification.
accessAdd = [0 1 1 0 1 0 1 1 0 1 1 1 1 1 0 1 1 0 0 ...
                            1 0 0 0 1 0 1 1 1 0 0 0 1]';
for iMode = 1:numMode

    phyMode = simMode{iMode};

    % Set signal to noise ratio (SNR) points based on mode.
    % For Coded PHYs (LE500K and LE125K), the code rate factor is included
    % in SNR calculation as 1/2 rate FEC encoder is used.
    if any(strcmp(phyMode,{'LE1M','LE2M'}))
        snrVec = EbNo - 10*log10(sps);
    else
        codeRate = 1/2;
        snrVec = EbNo + 10*log10(codeRate) - 10*log10(sps);
    end

%     parfor iSnr = 1:length(snrVec)  % Use 'parfor' to speed up the simulation
    for iSnr = 1:length(snrVec)   % Use 'for' to debug the simulation

        % Set random substream index per iteration to ensure that each
        % iteration uses a repeatable set of random numbers
        stream = RandStream('combRecursive','Seed',0);
        stream.Substream = iSnr;
        RandStream.setGlobalStream(stream);

        % Create an instance of error rate
        errorRate = comm.ErrorRate('Samples','Custom',...
                        'CustomSamples',1:(dataLen*bitsPerByte-1));

        % Create and configure the System objects for impairments
        initImp = helperBLEImpairmentsInit(phyMode,sps);

        % Create and configure the receiver System objects
        initRxParams = helperBLEReceiverInit(phyMode,sps,accessAdd);

        % Initialize error computation parameters
        [numErrs,perCnt] = deal(0);
        numPkt = 1;

        % Loop to simulate multiple packets
        while numErrs <= maxNumErrors && numPkt <= maxNumPackets

            % Generate BLE waveform
            txBits = randi([0 1],dataLen*bitsPerByte,1,'int8'); % Data bits generation
            chanIndex = randi([0 39],1,1); % Random channel index value for each packet
            txWaveform = bleWaveformGenerator(txBits,'Mode',phyMode,...
                                            'SamplesPerSymbol',sps,...
                                            'ChannelIndex',chanIndex,...
                                            'AccessAddress',accessAdd);

            % Define the RF impairment parameters
            initImp.pfo.FrequencyOffset = randsrc(1,1,-50e3:10:50e3); % In Hz, Max range is +/- 150 KHz
            initImp.pfo.PhaseOffset = randsrc(1,1,-10:5:10);          % In degrees
            initoff = 0.15*sps; % Static timing offset
            stepsize = 20*1e-6; % Timing drift in ppm, Max range is +/- 50 ppm
            initImp.vdelay = (initoff:stepsize:initoff+stepsize*(length(txWaveform)-1))'; % Variable timing offset
            initImp.dc = 20;     % Percentage w.r.t maximum amplitude value

            % Pass the generated waveform through RF impairments
            txImpairedWfm = helperBLEImpairmentsAddition(txWaveform,initImp);

            % Pass the transmitted waveform through AWGN channel
            rxWaveform = awgn(txImpairedWfm,snrVec(iSnr));

            % Recover data bits using practical receiver
            [rxBits,accessAddress] = helperBLEPracticalReceiver(rxWaveform,initRxParams,chanIndex);

            % Determine the BER and PER
            if(length(txBits) == length(rxBits))
                errors = errorRate(txBits,rxBits); % Outputs the accumulated errors
                ber(iMode,iSnr) = errors(1);       % Accumulated BER
                currentErrors = errors(2)-numErrs; % Number of errors in current packet
                if(currentErrors) % Check if current packet is in error or not
                    perCnt  = perCnt + 1;          % Increment the PER count
                end
                numErrs = errors(2);               % Accumulated errors
                numPkt = numPkt + 1;
            end
        end
        per(iMode,iSnr) = perCnt/(numPkt-1);

    disp(['Mode ' phyMode ', '...
        'Simulating for Eb/No = ', num2str(EbNo(iSnr)), ' dB' ', '...
        'BER:',num2str(ber(iMode,iSnr)), ', '...
        'PER:',num2str(per(iMode,iSnr))])
    end
end

markers = 'ox*s';
color = 'bmcr';
dataStr = {zeros(numMode,1)};
for iMode = 1:numMode
    subplot(2,1,1),semilogy(EbNo,ber(iMode,:).',['-' markers(iMode) color(iMode)]);
    hold on;
    dataStr(iMode) = simMode(iMode);

    subplot(2,1,2),semilogy(EbNo,per(iMode,:).',['-' markers(iMode) color(iMode)]);
    hold on;
    dataStr(iMode) = simMode(iMode);
end
subplot(2,1,1),
grid on;
xlabel('Eb/No (dB)');
ylabel('BER');
legend(dataStr);
title('BER of BLE under RF impairments');

subplot(2,1,2),
grid on;
xlabel('Eb/No (dB)');
ylabel('PER');
legend(dataStr);
title('PER of BLE under RF impairments');

[refBER,refPER,refEbNo] = deal(zeros(numMode,1));
headerLen = 2; % Header length in bytes
crcLen = 3; % CRC length in bytes
payloadLen = dataLen-headerLen-crcLen; % Payload length in bytes
for iMode = 1:numMode
    [refBER(iMode),refPER(iMode),refEbNo(iMode)] = ...
                        helperBLEReferenceResults(simMode(iMode),payloadLen);
    disp(['Mode ' simMode{iMode} ', '...
        'Reference Eb/No = ', num2str(refEbNo(iMode)), ' dB' ', '...
        'BER = ',num2str(refBER(iMode)), ', '...
        'PER = ',num2str(refPER(iMode)), ', '...
        'for payload length of ',num2str(payloadLen), ' bytes'])
end