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

simulationTime = 2*1e6; % Simulation time in microseconds
packetType     = 'DM1'; % Specify baseband packet type
mode           = 'BR'; % Specify PHY mode ('BR', 'EDR2M', 'EDR3M')
payloadLength  = 10; % Length of baseband packet in bytes

% Bluetooth frequency hopping
frequencyHop = bluetoothFrequencyHop;
frequencyHop.SequenceType = 'Connection Adaptive';

% Configure Bluetooth PHY transmission
phyTx = helperBluetoothPHY;
phyTx.Mode = mode;

% Configure Bluetooth PHY reception
phyRx = helperBluetoothPHY;
phyRx.Mode = mode;

% Specify as one of 'Generated' | 'BasebandFile' | 'None'
wlanInterferenceSource = 'BasebandFile';
wlanBBFilename = 'WLANNonHTDSSS.bb'; % Default baseband file
% Configure wireless channel
channel = helperBluetoothChannel;
channel.EbNo = 22; % Ratio of energy per bit (Eb) to the spectral noise density (No) in dB
channel.SIR = -20; % Signal to interference ratio in dB
if ~strcmpi(wlanInterferenceSource, 'None')
    % Generate the WLAN waveform
    wlanWaveform = helperBluetoothGenerateWLANWaveform(wlanInterferenceSource, wlanBBFilename);
    % Add the WLAN interference to Bluetooth channel
    addWLANWaveform(channel, wlanWaveform);
end

slotTime = 625; % Bluetooth slot duration in microseconds
% Simulation time in terms of slots
numSlots = floor(simulationTime/slotTime);
% Slot duration, including transmission and reception
slotValue = phyTx.slotsRequired(packetType)*2;
% Number of Master transmission slots
numMasterTxSlots = floor(numSlots/slotValue);
% Total number of Bluetooth physical channels
numBtChannels = 79;
% errorsBasic and errorsAdaptive store relevant bit and packet error
% information per channel. Each row stores the channel index, bit errors,
% packet errors, total bits, and BER per channel. errorsBasic and
% errorsAdaptive arrays store these values for basic frequency hopping
% and AFH, respectively.
[errorsBasic, errorsAdaptive] = deal(zeros(numBtChannels,5));
% Initialize first column with channel numbers
[errorsBasic(:,1), errorsAdaptive(:,1)] = deal(0:78);
% Initialize variables for calculating PER and BER
[berBasic, berAdaptive, bitErrors] = deal(0);
badChannels = zeros(1,0);
totalTransmittedPackets = numMasterTxSlots;
% Number of bits per octet
octetLength = 8;
% Sample rate and input clock used in PHY processing
samplePerSymbol = 88;
symbolRate = 1e6;
sampleRate = symbolRate*samplePerSymbol;
inputClock = 0;
% Store hop index
hopIndex = zeros(1, numMasterTxSlots);
% Index to hop index vector
hopIdx = 1;
% Baseband packet structure
basebandData = struct(...
    'LTAddr',1,             ... % Logical transport address
    'PacketType',packetType,... % Packet type
    'Payload',zeros(1,phyTx.MaxPayloadSize), ... % Payload
    'PayloadLength',0,  ... % Payload length
    'LLID',[0; 0],      ... % Logical link identifier
    'SEQN',0,           ... % Sequence number
    'ARQN',1,           ... % Acknowledgment flag
    'IsValid',false);   ... % Flag to identify the status of
    % cyclic redundancy check (CRC) and
% header error control (HEC)

% Bluetooth signal structure
bluetoothSignal = struct(...
    'PacketType',packetType, ... % Packet type
    'Waveform',[],           ... % Waveform
    'NumSamples',[],         ... % Number of samples
    'SampleRate',sampleRate, ... % Sample rate
    'SamplesPerSymbol',samplePerSymbol,      ... % Samples per symbol
    'Payload',zeros(1,phyTx.MaxPayloadSize), ... % Payload
    'PayloadLength',0, ... % Payload length
    'SourceID',0,      ... % Source identifier
    'Bandwidth',1,     ... % Bandwidth
    'NodePosition',[0 0 0], ... % Node position
    'CenterFrequency',centerFrequency(phyTx), ... % Center frequency
    'StartTime',0, ... % Waveform start time
    'EndTime',0,   ... % Waveform end time
    'Duration',0); ... % Waveform duration

% Clock ticks(one slot is 2 clock ticks)
clockTicks = slotValue*2;

% Spectrum analyzer for basic frequency hopping
spectrumAnalyzerBasic = dsp.SpectrumAnalyzer(...
    'Name','Bluetooth Basic Frequency Hopping', ...
    'ViewType','Spectrum and spectrogram', ...
    'TimeResolutionSource','Property', ...
    'TimeResolution',0.0005, ...
    'SampleRate',sampleRate, ...
    'TimeSpanSource','Property', ...
    'TimeSpan', 0.05, ...
    'FrequencyResolutionMethod', 'WindowLength', ...
    'WindowLength', 512, ...
    'AxesLayout', 'Horizontal', ...
    'FrequencyOffset',2441*1e6, ...
    'ColorLimits',[-20 15]);

% Spectrum analyzer for AFH
spectrumAnalyzerAdaptive = dsp.SpectrumAnalyzer(...
    'Name','Bluetooth Adaptive Frequency Hopping', ...
    'ViewType','Spectrum and spectrogram', ...
    'TimeResolutionSource','Property', ...
    'TimeResolution',0.0005, ...
    'SampleRate',sampleRate, ...
    'TimeSpanSource','Property', ...
    'TimeSpan',0.05, ...
    'FrequencyResolutionMethod','WindowLength', ...
    'WindowLength',512, ...
    'AxesLayout','Horizontal', ...
    'FrequencyOffset',2441*1e6, ...
    'ColorLimits',[-20 15]);

sprev = rng('default'); % Set random number generator seed
for slotIdx = 0:slotValue:numSlots-slotValue
    % Update clock
    inputClock = inputClock + clockTicks;

    % Frequency hopping
    [channelIndex,~] = nextHop(frequencyHop,inputClock);
    % PHY transmission
    stateTx = 1; % Transmission state
    TxBits = randi([0 1],payloadLength*octetLength,1);
    basebandData.Payload = TxBits;
    basebandData.PayloadLength = payloadLength;
    % Generate whiten initialization vector from clock
    clockBinary = comm.internal.utilities.de2biBase2RightMSB(inputClock,28);
    whitenInitialization = [clockBinary(2:7)'; 1];
    % Update the PHY with request from the baseband layer
    updatePHY(phyTx,stateTx,channelIndex,whitenInitialization,basebandData);
    % Initialize and pass elapsed time as zero
    elapsedTime = 0;
    [nextTxTime,btWaveform] = run(phyTx,elapsedTime); % Run PHY transmission
    run(phyTx, nextTxTime); % Update next invoked time

    % Channel
    bluetoothSignal.Waveform = btWaveform;
    bluetoothSignal.NumSamples = numel(btWaveform);
    bluetoothSignal.CenterFrequency = centerFrequency(phyTx);
    channel.ChannelIndex = channelIndex;
    bluetoothSignal = run(channel,bluetoothSignal,mode);
    distortedWaveform = bluetoothSignal.Waveform;

    % PHY reception
    stateRx = 2; % Reception state
    % Update the PHY with request from the baseband layer
    updatePHY(phyRx,stateRx,channelIndex,whitenInitialization);
    [nextRxTime,~] = run(phyRx,elapsedTime,bluetoothSignal);
    bluetoothSignal.NumSamples = 0;
    run(phyRx,nextRxTime,bluetoothSignal); % Run PHY reception
    chIdx = channelIndex + 1;

    % Calculate error rate upon successful decoding the packet
    if phyRx.Decoded
        rxBitsLength = phyRx.DecodedBasebandData.PayloadLength*octetLength;
        RxBits = phyRx.DecodedBasebandData.Payload(1:rxBitsLength);
        % BER calculation
        txSymLength = length(TxBits);
        rxSymLength = length(RxBits);
        minSymLength = min(txSymLength,rxSymLength);
        if minSymLength > 0
            bitErrors = sum(xor(TxBits(1:minSymLength),RxBits(1:minSymLength)));
            totalBits = minSymLength;
            % Bit errors found in channel
            errorsBasic(chIdx,2) = errorsBasic(chIdx,2) + bitErrors;
            % Total bits transmitted in channel
            errorsBasic(chIdx,4) = errorsBasic(chIdx,4) + totalBits;
        end
        if ~phyRx.DecodedBasebandData.IsValid || bitErrors
            % Packet errors found in channel
            errorsBasic(chIdx,3) = errorsBasic(chIdx,3) + 1;
        end
    else
        % Packet errors found in channel
        errorsBasic(chIdx,3) = errorsBasic(chIdx,3) + 1;
    end
    hopIndex(hopIdx) = channelIndex;
    hopIdx = hopIdx + 1;

    % Plot spectrum
    spectrumAnalyzerBasic(btWaveform + wlanWaveform(1:numel(btWaveform)));
    pause(0.01);
end

% Plot selected channel index per slot
figBasic = figure('Name','Basic frequency hopping');
axisBasic = axes(figBasic);
xlabel(axisBasic,'Slot');
ylabel(axisBasic,'Channel Index');
ylim(axisBasic,[0 numBtChannels+3]);
title(axisBasic,'Bluetooth Basic Frequency Hopping');
hold on;
plot(axisBasic,0:slotValue:numSlots-slotValue,hopIndex,'-o');

% Select 25 channels with highest packet errors as bad channels
[~,indexes] = sort(errorsBasic(:,3),'descend');
badChannelIdx = min(nnz(errorsBasic(:,3)),25);
if badChannelIdx ~= 0
    badChannels = indexes(1:badChannelIdx) - 1;
end
usedChannels = setdiff(0:numBtChannels-1,badChannels);

% BER per channel calculation
errorsBasic(:,5) = errorsBasic(:,2)./errorsBasic(:,4);
errorsBasic(:,5) = fillmissing(errorsBasic(:,5),'constant',0);
ber = nonzeros(errorsBasic(:,5));
if ~isempty(ber)
    berBasic = mean(ber);
end

% PER calculation
packetErrorsBasic = sum(errorsBasic(:,3));
perBasic = packetErrorsBasic/totalTransmittedPackets;

% Reset
hopIdx = 1;
[inputClock,bitErrors] = deal(0);
fprintf('PER of Bluetooth BR/EDR waveforms using basic frequency hopping: %.4f\n',perBasic);

fprintf('BER of Bluetooth BR/EDR waveforms using basic frequency hopping: %.4f\n',berBasic);

frequencyHop.SequenceType = 'Connection adaptive';
frequencyHop.UsedChannels = usedChannels;

%rng('default'); % Set random number generator seed
for slotIdx = 0:slotValue:numSlots-slotValue
    % Update clock
    inputClock = inputClock + clockTicks;

    % Frequency hopping
    [channelIndex,~] = nextHop(frequencyHop,inputClock);

    % PHY transmission
    stateTx = 1; % Transmission state
    TxBits = randi([0 1],payloadLength*octetLength,1);
    basebandData.Payload = TxBits;
    basebandData.PayloadLength = payloadLength;
    % Generate whiten initialization vector from clock
    clockBinary = comm.internal.utilities.de2biBase2RightMSB(inputClock,28);
    whitenInitialization = [clockBinary(2:7)'; 1];
    % Update the PHY with request from the baseband layer
    updatePHY(phyTx,stateTx,channelIndex,whitenInitialization,basebandData);
    % Initialize and pass elapsed time as zero
    elapsedTime = 0;
    [nextTxTime,btWaveform] = run(phyTx,elapsedTime); % Run PHY transmission
    run(phyTx,nextTxTime); % Update next invoked time

    % Channel
    bluetoothSignal.Waveform = btWaveform;
    bluetoothSignal.NumSamples = numel(btWaveform);
    bluetoothSignal.CenterFrequency = centerFrequency(phyTx);
    channel.ChannelIndex = channelIndex;
    bluetoothSignal = run(channel,bluetoothSignal,mode);
    distortedWaveform = bluetoothSignal.Waveform;

    % PHY reception
    stateRx = 2; % Reception state
    % Update the PHY with request from the baseband layer
    updatePHY(phyRx,stateRx,channelIndex,whitenInitialization);
    [nextRxTime,~] = run(phyRx,elapsedTime,bluetoothSignal);
    bluetoothSignal.NumSamples = 0;
    run(phyRx,nextRxTime,bluetoothSignal); % Run PHY reception
    chIdx = channelIndex + 1;
    
    % Calculate error rate upon successful decoding the packet
    if phyRx.Decoded
        rxBitsLength = phyRx.DecodedBasebandData.PayloadLength*octetLength;
        RxBits = phyRx.DecodedBasebandData.Payload(1:rxBitsLength);
        % BER calculation
        txSymLength = length(TxBits);
        rxSymLength = length(RxBits);
        minSymLength = min(txSymLength,rxSymLength);
        if minSymLength > 0
            bitErrors = sum(xor(TxBits(1:minSymLength),RxBits(1:minSymLength)));
            totalBits = minSymLength;
            % Bit errors found in channel
            errorsAdaptive(chIdx,2) = errorsAdaptive(chIdx,2) + bitErrors;
            % Total bits transmitted in channel
            errorsAdaptive(chIdx,4) = errorsAdaptive(chIdx,4) + totalBits;
        end

        if ~phyRx.DecodedBasebandData.IsValid || bitErrors
            % Packet errors found in channel
            errorsAdaptive(chIdx, 3) = errorsAdaptive(chIdx, 3) + 1;
        end
    else
        % Packet errors found in channel
        errorsAdaptive(chIdx, 3) = errorsAdaptive(chIdx, 3) + 1;
    end
    hopIndex(hopIdx) = channelIndex;
    hopIdx = hopIdx + 1;

    % Plot spectrum
    spectrumAnalyzerAdaptive(btWaveform + wlanWaveform(1:numel(btWaveform)));
    pause(0.01);
end

% Plot selected channel index per slot
figAdaptive = figure('Name','Adaptive Frequency Hopping');
axisAdaptive = axes(figAdaptive);
xlabel(axisAdaptive,'Slot');
ylabel(axisAdaptive,'Channel Index');
title(axisAdaptive,'Bluetooth Adaptive Frequency Hopping');
ylim(axisAdaptive,[0 numBtChannels+3]);
hold on;
plot(axisAdaptive,0:slotValue:numSlots-slotValue,hopIndex,'-o');

% BER per channel calculation
errorsAdaptive(:,5) = errorsAdaptive(:,2)./errorsAdaptive(:,4);
errorsAdaptive(:,5) = fillmissing(errorsAdaptive(:,5),'constant',0);
ber = nonzeros(errorsAdaptive(:,5));
if ~isempty(ber)
    berAdaptive = mean(ber);
end

% PER calculation
packetErrorsAdaptive = sum(errorsAdaptive(:,3));
perAdaptive = packetErrorsAdaptive/totalTransmittedPackets;
fprintf('PER of Bluetooth BR/EDR waveforms using adaptive frequency hopping: %.4f\n',perAdaptive);

fprintf('BER of Bluetooth BR/EDR waveforms using adaptive frequency hopping: %.4f\n',berAdaptive);

% Restore previous setting of random number generation
rng(sprev);