%% End-to-End Bluetooth Low Energy PHY Simulation with RF Impairments and Corrections
% This example shows how the Communications Toolbox(TM) Library for the
% Bluetooth(R) Protocol can be used to measure the bit error rate (BER) and
% packet error rate (PER) for different modes of Bluetooth Low Energy (BLE)
% [ <#11 1> ] physical layer (PHY) packet transmissions that have radio
% front-end (RF) impairments and additive white gaussian noise (AWGN) added
% to them.

% Copyright 2018-2019 The MathWorks, Inc.

%% Introduction
% Bluetooth special interest group (SIG) introduced BLE for low power short
% range communications. BLE devices operate in the globally unlicensed
% industrial, scientific and medical (ISM) band in the frequency range 2.4
% GHz to 2.485 GHz. BLE specifies a channel spacing of 2 MHz, which results
% in 40 RF channels. The BLE standard specifies the *Link* layer which
% includes both *PHY* and *MAC* layers. BLE applications include image and
% video file transfers between mobile phones, home automation, and the
% internet of things (IoT).
%
% This end-to-end BLE PHY simulation determines BER and PER performance of
% the four BLE PHY transmission modes with RF impairments and AWGN added to
% the transmission packets. Nested for loops are used to compute error
% rates for each transmission mode at several bit energy to noise density
% ratio (Eb/No) settings. Inside the Eb/No loop, multiple transmission
% packets are generated using the <matlab:help('bleWaveformGenerator')
% bleWaveformGenerator> function and altered with RF impairments and AWGN
% to accumulate the error rate statistics. Each packet is distorted by
% these RF impairments:
%
% * DC offset
% * Carrier frequency offset
% * Carrier phase offset
% * Timing drift
%
% White gaussian noise is added to the transmitted BLE waveforms. The noisy
% packets are processed through a practical BLE receiver that performs the
% following operations:
%
% * Automatic gain control (AGC)
% * DC removal
% * Carrier frequency offset correction
% * Matched filtering
% * Packet detection
% * Timing error correction
% * Demodulation and decoding
% * Dewhitening
%
% The processing steps for each packet are summarized in the following
% diagram:
%
% <<../BLEPHYFlow.png>>
%
% The synchronized packets are then demodulated and decoded to recover the
% data bits. These recovered data bits are compared with the transmitted
% data bits to determine the BER and PER. BER and PER curves are generated
% for the following four PHY transmission throughput modes supported in
% BLE:
%
% * Uncoded PHY with data rate of 1 Mbps (LE1M)
% * Uncoded PHY with data rate of 2 Mbps (LE2M)
% * Coded PHY with data rate of 500 Kbps (LE500K)
% * Coded PHY with data rate of 125 Kbps (LE125K)
%

%% Check for Support Package Installation

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

%% Initialize the Simulation Parameters
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

%%
% The number of packets tested at each Eb/No point is controlled by two
% parameters:
%
% # |maxNumErrors| is the maximum number of bit errors simulated at each
% Eb/No point. When the number of bit errors reaches this limit, the
% simulation at this Eb/No point is complete.
% # |maxNumPackets| is the maximum number of packets simulated at each
% Eb/No point and limits the length of the simulation if the bit error
% limit is not reached.
%
% The numbers chosen for |maxNumErrors| and |maxNumPackets| in this example
% will lead to a very short simulation. For meaningful results we recommend
% increasing these numbers.

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

%% Simulating for Each Eb/No Point
% This example also demonstrates how a <docid:matlab_ref#f71-813245 parfor>
% loop can be used instead of the <docid:matlab_ref#buhafgy-1 for> loop
% when simulating each Eb/No point to speed up a simulation.
% <docid:matlab_ref#f71-813245 parfor>, as part of the
% <docid:distcomp_doccenter#index Parallel Computing Toolbox>, executes
% processing for each Eb/No point in parallel to reduce the total
% simulation time. To enable the use of parallel computing for increased
% speed, comment out the 'for' statement and uncomment the 'parfor'
% statement below. If Parallel Computing Toolbox (TM) is not installed,
% 'parfor' will default to the normal 'for' statement.

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

%% Simulation Results
% This section presents the BER and PER results w.r.t the input Eb/No range
% for the considered PHY modes
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');

%% Reference Results
% This section generates the reference BER, PER and Eb/No values for each
% PHY mode based on the receiver sensitivity and corresponding BER as
% specified in section 4.1 of the Bluetooth specification [ <#11 1> ].
[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

%% Further Exploration
% The number of packets tested at each Eb/No value is controlled by
% |maxNumErrors| and |maxNumPackets| parameters. For statistically
% meaningful results these values should be larger than those presented in
% this example. To generate the figure below, the simulation ran using a
% data length of 128 bytes, |maxNumErrors = 1e3|, and |maxNumPackets = 1e4|
% for all the four transmission modes.
%
% <<../BLEPracticalReceiverExample.png>>
%
% The figure shows that the reference BER and PER can be obtained at lower
% Eb/No points compared to the reference Eb/No value given in the Bluetooth
% specification. In this example, only the following impairments are added
% and passed through AWGN channel.
%
% * DC offset
% * Carrier frequency offset
% * Carrier phase offset
% * Timing drift
%
% The reference Eb/No values generated based on the BLE specification
% include margin for RF impairments and fading channel conditions that are
% not modeled in this simulation. As a result, the simulation results here
% outperform the standard reference results. If you modify this example to
% include additional impairments such as frequency drift, fading, and
% interference in the simulation the BER and PER curves will move right
% towards the reference Eb/No values generated based on the BLE standard
% receiver characteristics in [ <#11 1> ], Volume 6, Section 4.1.

%% Appendix
% The helpers used in this example are:
%
% * <matlab:edit('helperBLEImpairmentsAddition.m')
% helperBLEImpairmentsAddition.m>: Adds RF impairments to the BLE waveform
% * <matlab:edit('helperBLEPracticalReceiver.m')
% helperBLEPracticalReceiver.m>: Demodulate and decodes the received
% signal
% * <matlab:edit('helperBLEReceiverInit.m')
% helperBLEReceiverInit.m>: Initializes BLE receiver parameters
% * <matlab:edit('helperBLEImpairmentsInit.m')
% helperBLEImpairmentsInit.m>: Initializes RF impairment parameters
% * <matlab:edit('helperBLEReferenceResults.m')
% helperBLEReferenceResults.m>: Generates reference BER, PER and Eb/No
% values
%
%% Summary
% This example simulates a BLE PHY packet transmissions that have RF
% impairments and AWGN added to them. It shows how to generate BLE
% waveforms, demodulate and decode data bits using practical receiver and
% compute the BER and PER.

%% Selected Bibliography
% # Volume 6 of the Bluetooth Core Specification, Version 5.0 Core System 
% Package [Low Energy Controller Volume].
